Custom Hamiltonians

This tutorial explains how to build custom Hamiltonians in QubitOS for advanced quantum control scenarios beyond the standard single-qubit gates.

Prerequisites

• Completed Pulse Generation Tutorial
• Completed GRAPE Optimizer Deep Dive
• Understanding of quantum Hamiltonians

---

Hamiltonian Basics

In quantum mechanics, the Hamiltonian \(H\) governs time evolution:

\[ i\hbar\frac{d}{dt}|\psi(t)\rangle = H(t)|\psi(t)\rangle \]

In quantum control, we decompose the Hamiltonian into:

\[ H(t) = H_0 + \sum_k u_k(t) H_k \]

Where:

• \(H_0\) is the drift Hamiltonian (always-on terms like qubit frequencies)
• \(H_k\) are control Hamiltonians (how pulses couple to the system)
• \(u_k(t)\) are the control amplitudes (what GRAPE optimizes)

---

Pauli String Notation

QubitOS uses a convenient string notation for building Hamiltonians from Pauli operators.

Basic Syntax

from qubitos.pulsegen.hamiltonians import parse_pauli_string

# Single Pauli operator on qubit 0
H = parse_pauli_string("X0", num_qubits=1)

# Pauli operator on qubit 1 in a 2-qubit system
H = parse_pauli_string("Z1", num_qubits=2)

# Tensor product of operators
H = parse_pauli_string("Z0 Z1", num_qubits=2)  # Z ⊗ Z

Coefficients and Sums

# Single term with coefficient
H = parse_pauli_string("0.5 * X0", num_qubits=1)

# Sum of terms
H = parse_pauli_string("0.5 * X0 + 0.3 * Z0", num_qubits=1)

# Subtraction
H = parse_pauli_string("1.0 * Z0 - 0.5 * X0", num_qubits=1)

Multi-Qubit Examples

# Two-qubit Ising coupling
H = parse_pauli_string("0.01 * Z0 Z1", num_qubits=2)

# Heisenberg interaction
H = parse_pauli_string("0.01 * X0 X1 + 0.01 * Y0 Y1 + 0.01 * Z0 Z1", num_qubits=2)

# Transverse field Ising model
H = parse_pauli_string(
    "1.0 * Z0 + 1.0 * Z1 + 0.5 * Z0 Z1 + 0.1 * X0 + 0.1 * X1",
    num_qubits=2
)

---

Building Hamiltonians Programmatically

For more complex systems, use the build_hamiltonian function:

from qubitos.pulsegen.hamiltonians import build_hamiltonian

# Build with Pauli strings
H0, Hc = build_hamiltonian(
    drift="5.0 * Z0",                    # Qubit frequency term
    controls=["X0", "Y0"],               # I and Q control
    num_qubits=1,
)

Using NumPy Arrays Directly

import numpy as np
from qubitos.pulsegen.hamiltonians import PAULI_X, PAULI_Y, PAULI_Z

# Define drift Hamiltonian as matrix
omega_qubit = 5.0  # GHz
H0 = omega_qubit * PAULI_Z / 2

# Control Hamiltonians
Hc = [PAULI_X / 2, PAULI_Y / 2]

# Use in optimizer
result = optimizer.optimize(
    target_unitary=target,
    num_qubits=1,
    drift_hamiltonian=H0,
    control_hamiltonians=Hc,
)

---

Standard Gate Unitaries

QubitOS provides target unitaries for common gates:

from qubitos.pulsegen.hamiltonians import get_target_unitary
import numpy as np

# Single-qubit gates
X_gate = get_target_unitary("X")       # Pauli-X
Y_gate = get_target_unitary("Y")       # Pauli-Y
Z_gate = get_target_unitary("Z")       # Pauli-Z
H_gate = get_target_unitary("H")       # Hadamard
S_gate = get_target_unitary("S")       # Phase gate
T_gate = get_target_unitary("T")       # T gate
SX_gate = get_target_unitary("SX")     # √X gate

# Rotation gates (parameterized)
RX_90 = get_target_unitary("RX", angle=np.pi/2)
RY_45 = get_target_unitary("RY", angle=np.pi/4)
RZ_180 = get_target_unitary("RZ", angle=np.pi)

# Two-qubit gates
CZ_gate = get_target_unitary("CZ", num_qubits=2)
CNOT_gate = get_target_unitary("CNOT", num_qubits=2)
ISWAP_gate = get_target_unitary("ISWAP", num_qubits=2)
SWAP_gate = get_target_unitary("SWAP", num_qubits=2)

Embedding Gates in Larger Systems

# X gate on qubit 1 in a 3-qubit system
X_on_q1 = get_target_unitary("X", num_qubits=3, qubit_indices=[1])

# CZ between qubits 0 and 2 in a 3-qubit system
CZ_02 = get_target_unitary("CZ", num_qubits=3, qubit_indices=[0, 2])

---

Custom Target Unitaries

For gates not in the standard library, define your own:

import numpy as np
from qubitos.pulsegen import GrapeOptimizer, GrapeConfig

# Custom gate: √Y gate
sqrtY = np.array([
    [1+1j, -1-1j],
    [1+1j,  1+1j],
], dtype=np.complex128) / 2

config = GrapeConfig(
    num_time_steps=100,
    duration_ns=50,
    target_fidelity=0.999,
)

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

print(f"√Y gate fidelity: {result.fidelity:.4f}")

Verify Your Unitary

Always check that custom unitaries are valid:

def is_unitary(U: np.ndarray, tol: float = 1e-10) -> bool:
    """Check if a matrix is unitary: U^dag @ U = I."""
    d = U.shape[0]
    identity = np.eye(d)
    product = U.conj().T @ U
    return np.allclose(product, identity, atol=tol)

# Verify
print(f"√Y is unitary: {is_unitary(sqrtY)}")

---

Two-Qubit Systems

Two-qubit optimization requires more care with Hamiltonians.

Drift Hamiltonian

A typical superconducting two-qubit system has:

\[ H_0 = \frac{\omega_0}{2}Z_0 + \frac{\omega_1}{2}Z_1 + J \cdot Z_0 Z_1 \]

import numpy as np
from qubitos.pulsegen.hamiltonians import (
    parse_pauli_string,
    build_hamiltonian,
    get_target_unitary,
)
from qubitos.pulsegen import GrapeOptimizer, GrapeConfig

# System parameters
omega_0 = 5.0   # Qubit 0 frequency (GHz)
omega_1 = 5.2   # Qubit 1 frequency (GHz)
J = 0.01        # Coupling strength (GHz)

# Build drift Hamiltonian
H0 = parse_pauli_string(
    f"{omega_0/2} * Z0 + {omega_1/2} * Z1 + {J} * Z0 Z1",
    num_qubits=2
)

Control Hamiltonians

For independent control of each qubit:

# Control on qubit 0: I and Q channels
Hc_0_I = parse_pauli_string("X0", num_qubits=2)
Hc_0_Q = parse_pauli_string("Y0", num_qubits=2)

# Control on qubit 1: I and Q channels
Hc_1_I = parse_pauli_string("X1", num_qubits=2)
Hc_1_Q = parse_pauli_string("Y1", num_qubits=2)

control_hamiltonians = [Hc_0_I, Hc_0_Q, Hc_1_I, Hc_1_Q]

Optimizing a CZ Gate

# Target: CZ gate
target = get_target_unitary("CZ", num_qubits=2)

# Configuration for 2-qubit gates
config = GrapeConfig(
    num_time_steps=200,
    duration_ns=100,
    max_iterations=500,
    target_fidelity=0.999,
    learning_rate=0.5,
)

optimizer = GrapeOptimizer(config)
result = optimizer.optimize(
    target_unitary=target,
    num_qubits=2,
    drift_hamiltonian=H0,
    control_hamiltonians=control_hamiltonians,
)

print(f"CZ fidelity: {result.fidelity:.4f}")
print(f"Converged: {result.converged}")

---

Advanced: XY Coupling

For systems with XY-type coupling (common in transmon qubits):

\[ H_{XY} = \frac{g}{2}(X_0 X_1 + Y_0 Y_1) \]

This enables iSWAP-family gates:

# XY coupling Hamiltonian
g = 0.02  # Coupling strength (GHz)
H_coupling = parse_pauli_string(
    f"{g/2} * X0 X1 + {g/2} * Y0 Y1",
    num_qubits=2
)

# Add to drift
H0 = parse_pauli_string(
    f"{omega_0/2} * Z0 + {omega_1/2} * Z1",
    num_qubits=2
)
H0 = H0 + H_coupling

# Optimize iSWAP
target = get_target_unitary("ISWAP", num_qubits=2)
result = optimizer.optimize(
    target_unitary=target,
    num_qubits=2,
    drift_hamiltonian=H0,
)

---

Transmon Anharmonicity

Real transmon qubits have anharmonicity that affects gate quality. While the current implementation uses a two-level approximation, you can account for leakage to the \(|2\rangle\) state by using a 3-level model:

import numpy as np

def transmon_drift_3level(omega: float, alpha: float) -> np.ndarray:
    """Build 3-level transmon drift Hamiltonian.

    Args:
        omega: Qubit frequency (GHz)
        alpha: Anharmonicity (GHz), typically negative

    Returns:
        3x3 drift Hamiltonian
    """
    # Energy levels: |0⟩ = 0, |1⟩ = ω, |2⟩ = 2ω + α
    H0 = np.diag([0, omega, 2*omega + alpha]).astype(np.complex128)
    return H0

def transmon_drive_3level() -> tuple[np.ndarray, np.ndarray]:
    """Build 3-level transmon control Hamiltonians.

    Returns:
        Tuple of (H_I, H_Q) control Hamiltonians
    """
    # X-like coupling: |0⟩↔|1⟩ and |1⟩↔|2⟩ with √2 factor
    H_I = np.array([
        [0, 1, 0],
        [1, 0, np.sqrt(2)],
        [0, np.sqrt(2), 0],
    ], dtype=np.complex128)

    # Y-like coupling
    H_Q = np.array([
        [0, -1j, 0],
        [1j, 0, -1j*np.sqrt(2)],
        [0, 1j*np.sqrt(2), 0],
    ], dtype=np.complex128)

    return H_I, H_Q

# Usage (note: requires 3x3 target unitary)
H0 = transmon_drift_3level(omega=5.0, alpha=-0.3)
H_I, H_Q = transmon_drive_3level()

# Embed X gate in 3-level space
X_3level = np.array([
    [0, 1, 0],
    [1, 0, 0],
    [0, 0, 1],  # |2⟩ unchanged (no leakage)
], dtype=np.complex128)

3-Level Optimization

The current GRAPE implementation assumes 2-level qubits. For 3-level optimization, you'll need to modify the optimizer or use a num_qubits value that gives the correct Hilbert space dimension.

---

Tensor Product Utilities

For building multi-qubit operators manually:

from qubitos.pulsegen.hamiltonians import (
    tensor_product,
    PAULI_I,
    PAULI_X,
    PAULI_Y,
    PAULI_Z,
)

# X ⊗ I (X on qubit 0, I on qubit 1)
X_I = tensor_product([PAULI_X, PAULI_I])

# I ⊗ Z (I on qubit 0, Z on qubit 1)  
I_Z = tensor_product([PAULI_I, PAULI_Z])

# Z ⊗ Z (ZZ coupling)
Z_Z = tensor_product([PAULI_Z, PAULI_Z])

# Three-qubit: X ⊗ I ⊗ Z
X_I_Z = tensor_product([PAULI_X, PAULI_I, PAULI_Z])

---

Rotation Gates

Generate arbitrary rotation gates:

from qubitos.pulsegen.hamiltonians import rotation_gate
import numpy as np

# Rotations around Pauli axes
RX_pi2 = rotation_gate("X", np.pi/2)   # R_X(π/2)
RY_pi4 = rotation_gate("Y", np.pi/4)   # R_Y(π/4)
RZ_pi = rotation_gate("Z", np.pi)      # R_Z(π)

# Verify: R_X(π) = -iX
RX_pi = rotation_gate("X", np.pi)
print(np.allclose(RX_pi, -1j * PAULI_X))  # True

Arbitrary Axis Rotation

For rotation around an arbitrary axis \(\hat{n} = (n_x, n_y, n_z)\):

\[ R_{\hat{n}}(\theta) = \cos(\theta/2)I - i\sin(\theta/2)(n_x X + n_y Y + n_z Z) \]

def rotation_arbitrary(n: tuple[float, float, float], theta: float) -> np.ndarray:
    """Rotation around arbitrary axis.

    Args:
        n: Unit vector (nx, ny, nz) defining rotation axis
        theta: Rotation angle in radians

    Returns:
        2x2 rotation matrix
    """
    nx, ny, nz = n
    c = np.cos(theta / 2)
    s = np.sin(theta / 2)

    return c * PAULI_I - 1j * s * (nx * PAULI_X + ny * PAULI_Y + nz * PAULI_Z)

# Example: rotation around (1,1,0)/√2 axis
axis = (1/np.sqrt(2), 1/np.sqrt(2), 0)
R_custom = rotation_arbitrary(axis, np.pi/2)

---

Complete Example: Optimizing a Custom Gate

Putting it all together:

import numpy as np
from qubitos.pulsegen import GrapeOptimizer, GrapeConfig
from qubitos.pulsegen.hamiltonians import (
    parse_pauli_string,
    get_target_unitary,
    rotation_gate,
)

# Define the system
omega = 5.0  # Qubit frequency (GHz)

H0 = parse_pauli_string(f"{omega/2} * Z0", num_qubits=1)
Hc = [
    parse_pauli_string("X0", num_qubits=1),
    parse_pauli_string("Y0", num_qubits=1),
]

# Define a custom gate: R_X(π/3) - an irrational angle
target = rotation_gate("X", np.pi/3)

# Configure optimizer
config = GrapeConfig(
    num_time_steps=100,
    duration_ns=50,
    max_iterations=300,
    target_fidelity=0.9999,
    learning_rate=0.5,
)

# Run optimization
optimizer = GrapeOptimizer(config)
result = optimizer.optimize(
    target_unitary=target,
    num_qubits=1,
    drift_hamiltonian=H0,
    control_hamiltonians=Hc,
)

# Report results
print(f"Gate: R_X(π/3)")
print(f"Fidelity: {result.fidelity:.6f}")
print(f"Converged: {result.converged}")
print(f"Iterations: {result.iterations}")
print(f"Pulse duration: {config.duration_ns} ns")
print(f"Pulse samples: {len(result.i_envelope)}")

---

Troubleshooting

Low Fidelity for Custom Gates

Issue	Solution
Fidelity stuck at ~0.5	Check Hamiltonian normalization and units
Fidelity oscillates	Reduce learning rate
Very slow convergence	Increase pulse duration or time steps
Never converges	Verify target unitary is achievable with your controls

Common Mistakes

1. Wrong units: QubitOS uses GHz for frequencies and MHz for pulse amplitudes
2. Missing identity: For multi-qubit systems, operators must be tensored with identity on other qubits
3. Unitary check: Custom targets must be unitary (\(U^\dagger U = I\))
4. Control mismatch: Ensure control Hamiltonians match your physical system

---

Summary

Task	Function/Class
Parse Pauli string	parse_pauli_string()
Build Hamiltonian	build_hamiltonian()
Get standard gate	get_target_unitary()
Rotation gate	rotation_gate()
Tensor product	tensor_product()
Embed in larger space	embed_gate()

---

Next Steps

• API Reference - Complete Hamiltonian API
• Notebooks - Interactive examples
• Architecture - System design overview