How to create dephasing and amplitude robust single-qubit gates

Incorporate robustness into the design of optimal pulses

Boulder Opal exposes a highly-flexible optimization engine for general-purpose gradient-based optimization. It can be directly applied to model-based control optimization for arbitrary-dimensional quantum systems. Boulder Opal's optimization engine also allows the user to design pulses that are robust against certain types of noise. In this notebook, we demonstrate how to produce single qubit gates that are simultaneously robust to dephasing and amplitude noises.

Summary optimization workflow

1. Define robustness condition in computational 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 using e.g. qctrl.operations.infidelity_pwc. To enforce robustness, we use the same function and just add a list with noise operators to the parameter noise_operators:

    cost = qctrl.operations.infidelity_pwc(

You would typically choose this list to be made of the operators characterizing the dominant noise in your system dynamics, such as the dephasing and control noise operators demonstrated in this guide. The addition of this parameter in qctrl.operations.infidelity_pwc ensures that the optimization cost will take into account both the infidelity with respect to the defined target operation and the robustness term.

2. Execute graph-based optimization

With the graph object created, an optimization can be run using the qctrl.functions.calculate_optimization function. The function returns the results of the optimization. Note this example code block uses naming that should be replaced with the naming used in your graph.

optimization_result = qctrl.functions.calculate_optimization(
    output_node_names=["alpha", "gamma"],

Worked example: single qubit gate robust to amplitude and dephasing noises

We present a detailed worked example of robust optimization in a single-qubit system. Specifically, we consider a single-qubit system represented by the following Hamiltonian:

\begin{align*} H(t) &= \frac{(1+\beta_{\gamma}(t))}{2}\left(\gamma(t)\sigma_{-} + \gamma^*(t)\sigma_{+}\right) + \frac{\alpha(t)}{2} \sigma_{z} + \eta(t) \sigma_{z} \,, \end{align*}

where $\gamma(t)$ and $\alpha(t)$ are, respectively, complex and real time-dependent pulses, $\sigma_{\pm}$ are the qubit ladder operators and $\sigma_{z}$ is the Pauli-Z operator. $\eta(t)$ and $\beta_{\gamma}(t)$ are small, slowly-varying stochastic processes corresponding to dephasing and amplitude noise, respectively.

The functions of time $\gamma(t)$ and $\alpha(t)$ are not predetermined, and instead are optimized by the Q-CTRL optimization engine in order to achieve some target operation.

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 operators and parameters:

# Define standard matrices
sigma_y = np.array([[0, -1j], [1j, 0]])
sigma_z = np.array([[1, 0], [0, -1]])
sigma_m = np.array([[0, 1], [0, 0]])

# Define physical constants
gamma_max = 2 * np.pi * 0.5e6  # Hz
alpha_max = 2 * np.pi * 0.5e6  # Hz
segment_count = 50
duration = 10e-6  # s
sinc_cutoff_frequency = 5e6

Below we show how to create a data flow graph for optimizing the single-qubit system described above. Comments in the code explain the details of each step. 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.

# Define the data flow graph describing the system
with qctrl.create_graph() as graph:

    # Create a complex piecewise-constant (PWC) signal, with optimizable modulus
    # and phase, representing gamma(t)
    gamma = qctrl.operations.complex_pwc_signal(
            count=segment_count, lower_bound=0, upper_bound=gamma_max
            upper_bound=2 * np.pi,
    # Create a real PWC signal, with optimizable amplitude, representing
    # alpha(t)
    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 PWC operators representing the Hamiltonian terms
    drive = qctrl.operations.pwc_operator_hermitian_part(rediscretized_gamma * sigma_m)
    shift = rediscretized_alpha * sigma_z / 2

    # Create a constant PWC operator representing the dephasing noise
    # (note that we scale by 1/duration to ensure consistent units between
    # the noise Hamiltonian and the control Hamiltonian)
    dephasing = qctrl.operations.constant_pwc_operator(
        duration=duration, operator=sigma_z / duration

    # Create the noise list with the dephasing and drive operators
    noise_list = [dephasing, drive]

    # Create the target operator
    target_operator =

    # Create infidelity
    cost = qctrl.operations.infidelity_pwc(
        hamiltonian=drive + shift,

Execute graph-based optimization

optimization_result = qctrl.functions.calculate_optimization(
    cost_node_name="robust_cost", output_node_names=["alpha", "gamma"], graph=graph
print("Optimized cost:\t", optimization_result.cost)
Your task calculate_optimization (action_id="534472") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_optimization (action_id="534472") has completed.
Optimized cost:	 1.29497772818305e-10

We may use the qctrlvisualizer package to plot the optimized pulses, which are available in the result object. By setting polar=False, we plot $\gamma(t) = I(t) + i Q(t)$ in terms of its real and imaginary parts.

        "$\\alpha$": optimization_result.output["alpha"],
        "$\\gamma$": optimization_result.output["gamma"],