How to create leakage-robust single-qubit gates

Design pulses that minimize leakage to unwanted states

Boulder Opal exposes a highly-flexible optimization engine for general-purpose gradient-based optimization. It allows one to define the physical problem in a higher dimensional quantum system while retaining the ability to optimize for a target operation in a particular subspace. In this notebook we demonstrate the optimization of a single-qubit gate in a superconducting transmon system with three levels.

Summary workflow

1. Define optimization subspace in the graph

The flexible Q-CTRL optimization engine expresses all optimization problems as data flow graphs, which describe how optimization variables (variables that can be tuned by the optimizer) are transformed into the cost function (the objective that the optimizer attempts to minimize). For an optimal control problem, the cost is typically given by the gate infidelity with respect to a user-defined target operation. You can restrict this target operation to a subspace of your full quantum system by multiplying it by the appropriate subspace projector. This ensures that the projected target treats population outside the desired subspace as leakage error. For example, the qubit projector in a three-dimensional space is

$$ P_{\rm{qubit}} = \begin{pmatrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix}.$$

The target operator is defined using

    target_operator =, filter_function_projector=qubit_projector

where full_target is your original target operator in the full space. The filter_function_projector parameter is optional and should be used when you also want to add robustness to noise in that subspace.

Once this is defined, we pass this operator to the infidelity using qctrl.operations.infidelity_pwc

        infidelity = qctrl.operations.infidelity_pwc(

2. Execute graph-based optimization

With the graph object created, an optimization can be run using the qctrl.functions.calculate_optimization function. The cost, the outputs, and the graph must be provided. The function returns the results of the optimization.

Worked example: Optimize single-qubit Hadamard in a qutrit system

In this example, we consider a qutrit system in which we effect a single-qubit gate robust against control noise while simultaneously minimizing leakage out of the computational subspace. The system is described by the following Hamiltonian:

\begin{align*} H(t) = & \frac{\chi}{2} (a^\dagger)^2 a^2 + (1+\beta(t))\left(\gamma(t) a + \gamma^*(t) a^\dagger \right) + \frac{\alpha(t)}{2} a^\dagger a \,, \end{align*}

where $\chi$ is the anharmonicity, $\gamma(t)$ and $\alpha(t)$ are, respectively, complex and real time-dependent pulses, $\beta(t)$ is a small, slowly-varying stochastic amplitude noise process, and $a = |0 \rangle \langle 1 | + \sqrt{2} |1 \rangle \langle 2 |$.

import matplotlib.pyplot as plt
import numpy as np
from qctrlvisualizer import plot_controls

from qctrl import Qctrl

# Starting a session with the API
qctrl = Qctrl()

We start by defining the operators and physical parameters:

# Define standard matrices
transmon_levels = 3
a = np.diag(np.sqrt(np.arange(1, transmon_levels)), k=1)
ad = a.T.conj()
ada = ad @ a
ad2a2 = ad @ ada @ a
hadamard = np.array(
    [[1.0, 1.0, 0], [1.0, -1.0, 0], [0, 0, np.sqrt(2)]], dtype=complex
) / np.sqrt(2)
qubit_projector = np.pad(np.eye(2), ((0, 1), (0, 1)), mode="constant")

# Define physical constants and parameters
chi = 2 * np.pi * -300.0 * 1e6  # Hz
gamma_max = 2 * np.pi * 30e6  # Hz
alpha_max = 2 * np.pi * 30e6  # Hz
segment_count = 50
duration = 100e-9  # s
sinc_cutoff_frequency = 300e6

Below we show how to create a data flow graph for optimizing the system described above. Note that we are using a filter to produce smooth pulses as explained in our How to add smoothing and band-limits to optimized controls user guide.

# Create graph object
with qctrl.create_graph() as graph:
    # Create the complex optimizable gamma(t) signal
    gamma = qctrl.operations.complex_pwc_signal(
            count=segment_count, lower_bound=0, upper_bound=gamma_max
            upper_bound=2 * np.pi,

    # Create the optimizable alpha(t) signal
    alpha = qctrl.operations.pwc_signal(
            count=segment_count, lower_bound=-alpha_max, upper_bound=alpha_max

    # Create filtered signals
    sinc_kernel = qctrl.operations.sinc_convolution_kernel(sinc_cutoff_frequency)
    gamma_filtered = qctrl.operations.convolve_pwc(pwc=gamma, kernel=sinc_kernel)
    alpha_filtered = qctrl.operations.convolve_pwc(pwc=alpha, kernel=sinc_kernel)
    rediscretized_gamma = qctrl.operations.discretize_stf(
    rediscretized_alpha = qctrl.operations.discretize_stf(

    # Create Hamiltonian terms
    anharmonicity = ad2a2 * chi / 2
    drive = qctrl.operations.pwc_operator_hermitian_part(2 * rediscretized_gamma * a)
    shift = rediscretized_alpha * ada / 2
    hamiltonian = anharmonicity + drive + shift

    # Create the target operator in the qubit subspace
    target_operator =, filter_function_projector=qubit_projector

    infidelity = qctrl.operations.infidelity_pwc(

We can now run the optimization and plot the resulting control pulses:

# Run the optimization
optimization_result = qctrl.functions.calculate_optimization(
    cost_node_name="infidelity", output_node_names=["alpha", "gamma"], graph=graph

print(f"\nOptimized cost:\t{optimization_result.cost:.3e}")

# Plot the optimized controls
        "$\\alpha$": optimization_result.output["alpha"],
        "$\\gamma$": optimization_result.output["gamma"],
Your task calculate_optimization (action_id="534473") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_optimization (action_id="534473") has completed.

Optimized cost:	6.078e-09