Engineering Low-Level Control: Implementing Custom Gates via Qiskit Pulse in Qiskit 1.0
In noisy intermediate-scale quantum (NISQ) computing, abstract quantum logic gates ($X$, $Y$, $Z$, $\text{CNOT}$) must eventually be translated into physical, continuous-time control signals. On superconducting transmon architectures, these signals manifest as microwave voltage pulses applied to drive lines routed directly into dilution refrigerators.
While standardized compilers like Qiskit Transpiler map abstract circuits into standard calibrated gate sets (e.g., $ECR, R_z, X, \sqrt{X}$), algorithmic execution often suffers from pulse-level inefficiencies. Standard pulse sequences are calibrated for general-purpose execution, leaving significant performance on the table when executing specific subroutines, suppressing crosstalk, or implementing non-standard rotations.
By bypassing standard compiler abstractions and utilizing Qiskit Pulse, developers can directly design, modulate, and calibrate pulse envelopes at the physical hardware layer. This guide explores the physical and mathematical foundations of pulse-level quantum control, derives the Derivative Removal by Adiabatic Gate (DRAG) technique, and demonstrates how to attach custom pulse schedules to quantum circuits in Qiskit 1.0.
The Physics of Transmon Control: Why Pulse Level Matters
A transmon qubit is a weakly anharmonic Josephson junction oscillator. Unlike an ideal two-level spin system, a transmon possesses an infinite hierarchy of energy levels $|0\rangle, |1\rangle, |2\rangle, \dots, |n\rangle$.
Energy
^
| |3> (Non-computational)
| ---
| w12 + delta
| --- |2> (Non-computational leakage state)
|
| w01 (Drive frequency)
| --- |1>
|
|
| --- |0>
+-------------------------------------> State
The simplified Hamiltonian of an uncoupled transmon modeled as a Duffing oscillator is:
$$H_0 = \hbar \omega_{01} a^\dagger a + \frac{\hbar \delta}{2} a^\dagger a^\dagger a a$$
where $a^\dagger$ and $a$ are the bosonic creation and annihilation operators, $\omega_{01}$ is the transition frequency between the computational ground state $|0\rangle$ and first excited state $|1\rangle$, and $\delta = \omega_{12} - \omega_{01}$ is the transmon anharmonicity (typically negative, on the order of $-2\pi \times 300\text{ MHz}$ to $-2\pi \times 350\text{ MHz}$).
When a microwave control field $V(t) = \text{Re}[\Omega(t) e^{-i \omega_d t}]$ at drive frequency $\omega_d$ is coupled to the qubit via a drive line, the drive Hamiltonian in the dipole approximation is:
$$H_d(t) = \hbar \left( \Omega(t) e^{-i \omega_d t} + \Omega^*(t) e^{i \omega_d t} \right) (a + a^\dagger)$$
Transforming into the rotating frame of the drive ($\omega_d = \omega_{01}$) and applying the Rotating Wave Approximation (RWA) yields:
$$H_{\text{rot}} = \hbar \frac{\delta}{2} a^\dagger a^\dagger a a + \frac{\hbar}{2} \left( \Omega(t) a^\dagger + \Omega^*(t) a \right)$$
The Problem: Leakage to non-computational state $|2\rangle$
If we drive the transition $|0\rangle \to |1\rangle$ with a fast pulse of short duration $T$, its spectral width in the frequency domain broadens proportionally to $\sim 1/T$. Because the $|1\rangle \to |2\rangle$ transition frequency lies at $\omega_{12} = \omega_{01} + \delta$, the tail of a short Gaussian pulse spectrum overlaps with $\omega_{12}$.
This causes unwanted population transfer (leakage) into the non-computational state $|2\rangle$ and introduces phase errors ($\sigma_z$ errors) on the qubit due to dynamic AC Stark shifts.
Mathematical Formulation of DRAG (Derivative Removal by Adiabatic Gate)
To mitigate off-resonant leakage while maintaining fast gate durations, Motzoi et al. (2009) introduced the DRAG pulse strategy.
Instead of a purely real Gaussian pulse envelope $\Omega(t) = \Omega_x(t)$, we introduce a complex quadrature component $\Omega_y(t)$:
$$\Omega(t) = \Omega_x(t) + i \Omega_y(t)$$
We perform a time-dependent unitary transformation (Schrieffer-Wolff transformation) $U_S(t) = \exp(-i S(t))$ into a frame that eliminates driving terms to state $|2\rangle$ to first order in $1/\delta$:
$$S(t) = -\frac{\dot{\Omega}_x(t)}{2 \delta} (a^\dagger a^\dagger a - a^\dagger a a)$$
Under this frame transformation, the transformed Hamiltonian $H_{\text{eff}} = U_S H_{\text{rot}} U_S^\dagger - i U_S \dot{U}_S^\dagger$ cancels the matrix element driving $|1\rangle \to |2\rangle$ when the quadrature drive component is set precisely to the time derivative of the in-phase envelope:
$$\Omega_y(t) = -\frac{\beta}{\delta} \dot{\Omega}_x(t)$$
where $\beta$ is a dimensionless scaling factor (typically tuned near $\beta \approx 1$ in theoretical models and calibrated experimentally to account for higher-order shifts and transfer function distortions).
For a standard Gaussian envelope defined by amplitude $A$, standard deviation $\sigma$, and center $t_0$:
$$\Omega_x(t) = A \exp\left( -\frac{(t - t_0)^2}{2\sigma^2} \right)$$
$$\Omega_y(t) = - \beta \cdot \frac{t - t_0}{\sigma^2} \cdot A \exp\left( -\frac{(t - t_0)^2}{2\sigma^2} \right)$$
Practical Implementation in Qiskit 1.0
Below is a complete script demonstrating how to define custom gate objects, construct a low-level pulse schedule using Qiskit Pulse, attach calibrations to a circuit, and simulate the quantum state evolution under 3-level transmon Hamiltonian dynamics using NumPy and SciPy.
import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit
from qiskit.circuit import Gate
from qiskit import pulse
from qiskit.pulse import DriveChannel, Drag, build
# =====================================================================
# 1. Quantum Circuit & Custom Gate Calibration Construction
# =====================================================================
class CustomSXGate(Gate):
"""Custom pi/2 rotation gate around X-axis (sqrt(X)) using DRAG."""
def __init__(self, label=None):
super().__init__('custom_sx', 1, [], label=label)
def build_custom_calibrated_circuit():
"""Builds a QuantumCircuit and attaches custom pulse calibration."""
qc = QuantumCircuit(1, 1)
custom_sx = CustomSXGate()
# Append custom gate to circuit
qc.append(custom_sx, [0])
qc.measure(0, 0)
# Hardware pulse parameters
duration = 160 # pulse duration in samples (dt units)
amp = 0.18 # peak pulse amplitude [0, 1]
sigma = 40 # Gaussian standard deviation
beta = 1.85 # DRAG parameter
# Construct low-level pulse schedule using pulse builder
with build(name="custom_sx_schedule") as custom_sched:
d0 = DriveChannel(0)
# Create DRAG pulse envelope
drag_pulse = Drag(duration=duration, amp=amp, sigma=sigma, beta=beta)
pulse.play(drag_pulse, d0)
# Attach pulse calibration to circuit for qubit 0
qc.add_calibration(custom_sx, (0,), custom_sched)
return qc, custom_sched
# Generate circuit and schedule
circuit, schedule = build_custom_calibrated_circuit()
print("--- Quantum Circuit Representation ---")
print(circuit)
print("\n--- Attached Pulse Schedule ---")
print(schedule)
# =====================================================================
# 2. Physics Simulation: 3-Level Transmon Dynamics (DRAG vs Gaussian)
# =====================================================================
def simulate_transmon_dynamics(use_drag=True):
"""
Simulates time-dependent Schrödinger equation for a 3-level transmon
system subjected to a Gaussian or DRAG pulse envelope.
"""
# System Constants
dt = 0.222e-9 # Qiskit sampling interval (222 ps)
duration_samples = 160 # 160 dt = ~35.5 ns
t_total = duration_samples * dt
delta = -2 * np.pi * 330e6 # Transmon anharmonicity (-330 MHz in rad/s)
rabi_scale = 2 * np.pi * 45e6 # Coupling scaling (45 MHz peak Rabi rate)
# Time vector
t = np.linspace(0, t_total, duration_samples)
t0 = t_total / 2
sigma = 40 * dt
amp = 0.18
beta = 1.85 / (2 * np.pi * 330e6) # Dimensional DRAG parameter
# Envelope components
gauss = amp * np.exp(-0.5 * ((t - t0) / sigma) ** 2)
dgauss = -((t - t0) / (sigma ** 2)) * gauss
Omega_x = gauss
Omega_y = -beta * dgauss if use_drag else np.zeros_like(gauss)
# 3-level Transmon Operators
# |0> = [1,0,0]^T, |1> = [0,1,0]^T, |2> = [0,0,1]^T
H_0 = np.diag([0.0, 0.0, delta]) # Static Hamiltonian in rotating frame
a = np.array([[0, 1, 0],
[0, 0, np.sqrt(2)],
[0, 0, 0]], dtype=complex)
a_dag = a.T.conj()
# Time evolution starting from ground state |0>
psi = np.array([1.0, 0.0, 0.0], dtype=complex)
step_dt = dt
pop_0, pop_1, pop_2 = [], [], []
for i in range(duration_samples):
# Time-dependent drive Hamiltonian: H_d = 0.5 * (Omega * a^\dagger + Omega* * a)
Omega_t = rabi_scale * (Omega_x[i] + 1j * Omega_y[i])
H_drive = 0.5 * (Omega_t * a_dag + np.conj(Omega_t) * a)
H_total = H_0 + H_drive
# Propagation via matrix exponential over timestep dt
U_step = la.expm(-1j * H_total * step_dt)
psi = U_step @ psi
# Record populations
pop_0.append(np.abs(psi[0])**2)
pop_1.append(np.abs(psi[1])**2)
pop_2.append(np.abs(psi[2])**2)
return t * 1e9, np.array(pop_0), np.array(pop_1), np.array(pop_2)
# Run simulations
t_ns, p0_g, p1_g, p2_g = simulate_transmon_dynamics(use_drag=False)
_, p0_d, p1_d, p2_d = simulate_transmon_dynamics(use_drag=True)
print("\n--- Simulation Results ---")
print(f"Standard Gaussian Peak Leakage (|2> Population): {np.max(p2_g):.6e}")
print(f"DRAG Corrected Peak Leakage (|2> Population) : {np.max(p2_d):.6e}")
print(f"Leakage Suppression Ratio : {np.max(p2_g) / np.max(p2_d):.2f}x")
Hardware Limitations & Real-World Considerations
While pulse-level calibration unlocks fine-grained quantum control, physical execution on quantum hardware introduces strict system constraints:
1. Arbitrary Waveform Generator (AWG) Sample Granularity
Control electronics map digital pulse definitions into analog RF waveforms via high-speed DACs. Hardware backends enforce strict constraints on pulse waveforms: - Sample Time Step ($dt$): On current IBM Quantum backends, $dt \approx 0.222\text{ ns}$ ($4.5\text{ GSamples/s}$). Pulse durations must be defined as integer multiples of $dt$. - Alignment Constraints: Pulse durations and start times must satisfy alignment constraints (e.g., duration must be a multiple of 16 samples, $16 \cdot dt \approx 3.55\text{ ns}$) due to FPGA buffer alignment requirements.
2. Environmental Noise & Calibration Drift
Physical pulse shape parameters ($\text{amp}, \beta, \sigma$) are sensitive to environmental variations: - $T_1$ Relaxation & $T_2^*$ Dephasing: Shorter pulses reduce $T_1$ decay during gate execution but require higher peak amplitude $\Omega_{\max}$, increasing power broadening and state leakage. - Thermal Fluctuation & Phase Drift: System drift over hours causes drive frequency offsets ($\Delta \omega = \omega_d - \omega_{01}$), requiring periodic automated recalibration via automated tune-up pipelines (e.g., ORBIT or Randomized Benchmarking).
3. Architecture Evolution Beyond Qiskit 1.x
In Qiskit 1.0+, qiskit.pulse remains supported for low-level schedule creation and circuit calibration attachment. However, modern physics-based pulse simulation and optimal control workflows (such as GRAPE and CRAB algorithms) are increasingly transitioned to specialized high-performance simulation frameworks such as Qiskit Dynamics.