Superconducting qubits: improving the performance of single qubit gates

Increasing robustness against dephasing and control noise using Q-CTRL pulses

The Q-CTRL Python Package enables the user to create pulses that achieve a target operation on a given quantum hardware, and that are simultaneously robust against dephasing and control noise. The first section of the notebook presents the various steps necessary to create dephasing-robust pulses and characterize their performance using Q-CTRL tools. The second produces experimental validation of the robust pulses through their implementation on the IBM-Q hardware (note that executing hardware-related cells requires a valid IBM account). This section also contains all the hardware calibration required for proper pulse implementation. The third goes through the steps described in the previous two sections and presents the design, verification and implementation of pulses that are robust against control noise. Finally, we have a section where we discuss other challenges in superconducting qubit architectures and provide links to other notebooks where
we address them using Q-CTRL tools.

Some cells in this notebook require an account with IBM-Q to execute correctly. If you want to run them, please go to the IBM-Q experience to setup an account.

Table of Contents

  1. Creation and verification of dephasing-robust Q-CTRL pulses

    1. Creating robust pulses
    2. Time evolution of the driven qubit
    3. Robustness verification with filter functions and a quasi-static scan
  2. Experimental implementation on IBM-Q

    1. Pulse calibration: mapping $(I, Q)$ values onto hardware amplitude inputs
    2. Experimental time evolution of the driven qubit
    3. Robustness verification with a quasi-static scan
  3. Other noise sources: Q-CTRL pulses robust against control error

    1. Creating robust pulses
    2. Time evolution of the driven qubit: simulations and experiment
    3. Robustness verification with filter functions and a quasi-static scan
  4. Q-CTRL solutions for complex problems on superconducting qubits architecture

Imports and initialization

# Choose to run experiments or to use saved data
from pathlib import Path
import time

#disable warning from qiskit
import warnings
warnings.filterwarnings('ignore')

use_saved_data = True
if use_saved_data == False:
    timestr = time.strftime("%Y%m%d-%H%M%S")
    print("Time label for data saved throughout this experiment:" + timestr)
data_folder = Path("resources/")

# General imports
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from scipy.optimize import curve_fit
from scipy import interpolate
import pickle

# Parameters
giga = 1.0e9
mega = 1.0e6
micro = 1.0e-6
nano = 1.0e-9
scale_factor = 1e-14
colors = {"Unconstrained": "#BF04DC", "Bound slew": "#680CE9", "Square": "#000000", 'IBM default X-gate':'r'}

# Q-CTRL imports
from qctrl import Qctrl
from qctrlvisualizer import plot_controls

# Q-CTRL auxiliary functions
def get_detuning_scan(scan_array, control):

    durations = [segment["duration"] for segment in control["I"]]
    duration = np.sum(durations)

    I_values = [segment["value"] for segment in control["I"]]
    Q_values = [segment["value"] for segment in control["Q"]]

    # Define system object
    system = qctrl.factories.systems.create(name="qubit", hilbert_space_dimension=2)

    initial_state = np.array([1.0, 0.0])

    shift_I = qctrl.factories.shift_controls.create(
        name="I", system=system, operator=sigma_x / 2.0)
    shift_I_pulse = qctrl.factories.custom_pulses.create(
        control=shift_I,
        segments=[{"duration": d, "value": v} for d, v in zip(durations, I_values)])

    shift_Q = qctrl.factories.shift_controls.create(
        name="Q", system=system, operator=sigma_y / 2.0)
    shift_Q_pulse = qctrl.factories.custom_pulses.create(
        control=shift_Q,
        segments=[{"duration": d, "value": v} for d, v in zip(durations, Q_values)])

    drifts = qctrl.factories.additive_noises.create(
        name="drift", system=system, operator=sigma_z)

    quasi_static_function = qctrl.factories.quasi_static_functions.create(
        x_noise=drifts, x_coefficients=scan_array)

    # Define target. Note that setting target unitary and projection operators as below,
    # we ensure that the produced "infidelities" give the excited state population of
    # the final state after propagating a pure ground state. From this we may extract
    # the final ground state population, which is the quantity of interest.
    target = qctrl.factories.targets.create(
        system=system, unitary_operator=np.eye(2), projection_operator=np.array([[1, 0], [0, 0]]))

    # Compute quasi-static function
    result = qctrl.services.quasi_static_functions.calculate(quasi_static_function,)

    # Extract infidelities (excited state population)
    noises_and_infidelities = np.array(
        [[sampled_point["coefficients"][0], sampled_point["infidelity"]]
         for sampled_point in result.sampled_points])

    noises_and_infidelities_sorted = noises_and_infidelities[
        noises_and_infidelities.argsort(axis=0)[:, 0]]

    # return infidelity (ground state population)
    return 1 - noises_and_infidelities_sorted[:, 1]


def get_amplitude_scan(amplitudes, control):
    durations = [segment["duration"] for segment in control["I"]]
    duration = np.sum(durations)
    I_values = [segment["value"] for segment in control["I"]]
    Q_values = [segment["value"] for segment in control["Q"]]

    # Define system object
    system = qctrl.factories.systems.create(name="qubit", hilbert_space_dimension=2)

    initial_state = np.array([1.0, 0.0])

    drive = qctrl.factories.drive_controls.create(
        name="sigma_m", system=system, operator=sigma_m / 2)

    drive_pulse = qctrl.factories.custom_pulses.create(
        control=drive,
        segments=[
            {"duration": v, "value": I_values[idx] + Q_values[idx] * 1j}
            for idx, v in enumerate(durations)])

    amplitude_noise = qctrl.factories.control_noises.create(
        name="amplitude", control=drive)

    quasi_static_function = qctrl.factories.quasi_static_functions.create(
        x_noise=amplitude_noise, x_coefficients=amplitudes)

    # Define target. As above, this target allows us to calculate the ground
    # state populations.
    target = qctrl.factories.targets.create(
        system=system, unitary_operator=np.eye(2), projection_operator=np.array([[1, 0], [0, 0]]))

    # Compute quasi-static function
    result = qctrl.services.quasi_static_functions.calculate(quasi_static_function,)

    # Extract infidelities
    noises_and_infidelities = np.array(
        [[sampled_point["coefficients"][0], sampled_point["infidelity"]]
         for sampled_point in result.sampled_points])

    noises_and_infidelities_sorted = noises_and_infidelities[
        noises_and_infidelities.argsort(axis=0)[:, 0]]

    # return infidelities
    return 1 - noises_and_infidelities_sorted[:, 1]


def simulations(control):

    durations = [segment["duration"] for segment in control["I"]]
    I_values = [segment["value"] for segment in control["I"]]
    Q_values = [segment["value"] for segment in control["Q"]]
    duration = sum(durations)

    # Define system object
    system = qctrl.factories.systems.create(name="qubit", hilbert_space_dimension=2)

    shift_I = qctrl.factories.shift_controls.create(
        name="I", system=system, operator=sigma_x / 2.0)
    shift_I_pulse = qctrl.factories.custom_pulses.create(
        control=shift_I,
        segments=[{"duration": d, "value": v} for d, v in zip(durations, I_values)])

    shift_Q = qctrl.factories.shift_controls.create(
        name="Q", system=system, operator=sigma_y / 2.0)
    shift_Q_pulse = qctrl.factories.custom_pulses.create(
        control=shift_Q,
        segments=[{"duration": d, "value": v} for d, v in zip(durations, Q_values)])

    simulation = qctrl.factories.coherent_simulations.create(
        system=system,
        point_times=np.linspace(0, duration, 200),
        initial_state_vector=np.array([1.0, 0.0]))

    simulation_result = qctrl.services.coherent_simulations.run(system)

    gate_times = np.array(
        [frame["time"]
         for frame in simulation_result.simulations[0].trajectories[0].frames])
    state_vectors = np.array(
        [frame["state_vector"]
         for frame in simulation_result.simulations[0].trajectories[0].frames])

    return gate_times, state_vectors

# Q-CTRL login
qctrl = Qctrl()


# IBM-Q imports
from qiskit import IBMQ

from qiskit.tools.jupyter import *
from qiskit.tools.monitor import job_monitor

from qiskit.compiler import assemble
from qiskit.pulse.commands import SamplePulse

import qiskit.pulse as pulse
import qiskit.pulse.pulse_lib as pulse_lib

from qiskit.pulse import (DriveChannel, MeasureChannel,
                          AcquireChannel, Acquire, MemorySlot)

# IBM-Q auxiliary functions
def get_closest_multiple_of_16(num):
    return int(num + 8) - (int(num + 8) % 16)

def fit_function(x_values, y_values, function, init_params):
    fitparams, conv = curve_fit(function, x_values, y_values, init_params)
    y_fit = function(x_values, *fitparams)

    return fitparams, y_fit

def get_amplitude(vec):
    i_signal = np.imag(vec)
    r_signal = np.real(vec)

    mvec = [np.mean(r_signal), np.mean(i_signal)]
    src_mat = np.vstack((r_signal - mvec[0], i_signal - mvec[1])).T
    (_, _, v_mat) = np.linalg.svd(src_mat)
    dvec = v_mat[0, 0:2]
    if dvec.dot(mvec) < 0:
        dvec = -dvec
    return src_mat.dot(dvec)

# IBM credentials and backend selection
IBMQ.enable_account('YOUR TOKEN HERE')
provider = IBMQ.get_provider(hub="ibm-q", group="open", project="main")
backend = provider.get_backend("ibmq_armonk")
backend_defaults = backend.defaults()
backend_config = backend.configuration()

# Backend properties
qubit = 0
qubit_freq_est = backend_defaults.qubit_freq_est[qubit]
dt = backend_config.dt
print(f"Qubit: {qubit}")
print(f"Hardware sampling time: {dt/nano} ns")
print(f"Qubit frequency estimate: {qubit_freq_est/giga} GHz")
Qubit: 0
Hardware sampling time: 0.2222222222222222 ns
Qubit frequency estimate: 4.974295283079766 GHz

Return to Table of Contents

1. Creation and verification of dephasing-robust Q-CTRL pulses

In this section, we use BOULDER OPAL to perform optimizations to achieve an X-gate operation on the qubit. The total Hamiltonian of the driven quantum system is:

$$ H(t) = \left(1 + \beta_\gamma (t)\right)H_c(t) + \eta(t) \sigma_z, $$

where $\beta_\gamma(t)$ is a fractional time-dependent amplitude fluctuation process, $\eta(t)$ is a small stochastic slowly-varying dephasing noise process, and $H_c(t)$ is the control term given by

\begin{align*} H_c(t) = & \frac{1}{2}\left(\gamma(t) \sigma_- + \gamma^*(t) \sigma_+ \right) \\ = & \frac{1}{2}\left(I(t) \sigma_x + Q(t) \sigma_y \right). \end{align*}

Here, $\gamma(t) = I(t) + i Q(t)$ is the time-dependent complex-valued control pulse waveform and $\sigma_k$, $k=x,y,z$, are the Pauli matrices. The technical documentation explains the conventions Q-CTRL uses when splitting the total Hamiltonian of a system into a set of controls and noises.

In this notebook, we create optimal pulses using two different control scenarios:

  • Unconstrained - no constraints in the pulse shape.
  • Constrained slew rate - with a bounded rate of change in the pulse amplitude.

These pulses are constructed with the hardware backend characteristics in mind. The duration of each segment of the piecewise constant control pulse, for example, is required to be an integer multiple of the backend sampling time, dt. This guarantees that Q-CTRL pulses can be properly implemented given the time resolution of the hardware. In the example below, we create two robust pulses, one with 128 and the other with 256.

Return to Table of Contents

Creating robust pulses

We begin setting up the optimization by defining basic operators and constants. The $I$ and $Q$ terms of the control Hamiltonian are also defined in the cell below, together with the specific optimization contraints.

sigma_z = np.array([[1.0, 0.0], [0.0, -1.0]], dtype=np.complex)
sigma_x = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex)
sigma_y = np.array([[0.0, -1.0j], [1.0j, 0.0]], dtype=np.complex)
sigma_m = np.array([[0.0, 0.0], [1.0, 0.0]], dtype=np.complex)
X_gate = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex)

number_of_segments = {}
segment_scale = {}
duration_int = {}
duration = {}
rabi_rotation = np.pi
number_of_runs = 20
omega_max = 2 * np.pi * 8.5e6
I_max = omega_max / np.sqrt(2)
Q_max = omega_max / np.sqrt(2)

# Check if hardware time resolution has been defined 
# otherwise use default dt for armonk backend   
try:
    dt
except NameError:
    dt = 2/9 * 1e-9
    
number_of_segments['Unconstrained'] = 64
segment_scale['Unconstrained'] = 20
number_of_segments['Bound slew'] = 256
segment_scale['Bound slew'] = 5
scheme_names = ["Unconstrained", "Bound slew"]
shifts = {}

for shift_name in scheme_names:
    duration_int[shift_name] = number_of_segments[
        shift_name] * segment_scale[shift_name]
    duration[shift_name] = duration_int[shift_name] * dt
    shifts[shift_name] = {'I': {'operator': sigma_x/2.,
                          'segments_count': number_of_segments[shift_name],
                          'amplitude_upper_bound': I_max},
                          'Q': {'operator': sigma_y/2.,
                                'segments_count': number_of_segments[shift_name],
                                'amplitude_upper_bound': Q_max}}
shifts['Bound slew']['I']['maximum_slew_rate'] = I_max/10
shifts['Bound slew']['Q']['maximum_slew_rate'] = Q_max/10

dephasing_term = {'operator': sigma_z,'noise': True}

The optimization can now be performed by a single call. Note that the execution cells in this notebook contain an option (set at the initialization cell) to run all the commands from scratch or to use previously saved data.

if use_saved_data == False:
    robust_dephasing_controls = {}
    for shift_name, shift in shifts.items(): 
        # Create system  
        system = qctrl.factories.systems.create(name="system " + shift_name, 
                                                hilbert_space_dimension=2) 
         
        # Create shifts 
        shift_I = qctrl.factories.shift_controls.create( 
            name="I", 
            system=system, 
            operator=shift['I']['operator']) 
        
        shift_Q = qctrl.factories.shift_controls.create(
            name="Q", 
            system=system, 
            operator=shift['Q']['operator']) 
         
        # Create pulses 
        pulse_I = qctrl.factories.optimum_pulses.create( 
            control=shift_I, 
            upper_bound=shift['I']['amplitude_upper_bound'], 
            fixed_modulus=False, 
            segment_count=shift['I']['segments_count'], 
            duration=duration[shift_name], 
            maximum_slew_rate=shift['I'].get('maximum_slew_rate', None)) 
         
        pulse_Q = qctrl.factories.optimum_pulses.create( 
            control=shift_Q, 
            upper_bound=shift['Q']['amplitude_upper_bound'], 
            fixed_modulus=False, 
            segment_count=shift['Q']['segments_count'], 
            duration=duration[shift_name], 
            maximum_slew_rate=shift['Q'].get('maximum_slew_rate', None)) 
         
        # Create noises 
        dephasing = qctrl.factories.additive_noises.create( 
            name="Dephasing", 
            system=system, 
            operator=sigma_z) 
         
        # Create target 
        target = qctrl.factories.targets.create( 
            unitary_operator=X_gate, 
            projection_operator=np.identity(2), 
            system=system) 
         
        # Run optimization 
        system = qctrl.services.robust_optimization.run(system=system) 
         
        # Store results 
        robust_dephasing_controls[shift_name] = {control.name: control.pulse.segments  
                                                     for control in system.controls if control.pulse} 
        print(shift_name + ' cost = ', system.cost, '\n') 
    filename = 'robust_dephasing_controls_'+ timestr 
    outfile = open(filename,'wb') 
    pickle.dump(robust_dephasing_controls,outfile) 
    outfile.close() 
else:
    filename = data_folder / "robust_dephasing_controls_demo"
    infile = open(filename, "rb")
    robust_dephasing_controls = pickle.load(infile)
    infile.close()

We also create a 'primitive' benchmark pulse: a square pulse with the same duration as the optimized controls to serve as a reference.

number_of_segments['Square'] = 1
segment_scale['Square'] = 1280
duration_int['Square'] = number_of_segments['Square'] * segment_scale['Square']
duration['Square'] = duration_int['Square'] * dt
square_pulse_durations = np.array([duration['Square']])
square_pulse_values_control = np.array([rabi_rotation / duration['Square']])

square_sequence = {"I": [{"duration": d, "value": v}
                         for d, v in zip(square_pulse_durations, square_pulse_values_control)],
                   "Q": [{"duration": d, "value": v}
                         for d, v in zip(square_pulse_durations, np.zeros([1]))]}
scheme_names.append('Square')

Now, robust and primitive controls are combined into a single dictionary and plotted.

gate_schemes = {**robust_dephasing_controls, **{"Square": square_sequence}}

for scheme_name, gate in gate_schemes.items():
    print(scheme_name)
    plot_controls(plt.figure(), gate)
    plt.show()
Unconstrained
Bound slew
Square

The control plots represent the $I(t)$ and $Q(t)$ components of the control pulses for: unconstrained (top) and bound slew (middle) optimizations, as well as for the primitive control (bottom). The unconstrained solution may present a highly jagged control profile while the bound slew pulse provides a smoother control by constraining how much the control values can change from one segment to the next.

Return to Table of Contents

Time evolution of the driven qubit

Another useful way to characterize the pulses is to use the Q-CTRL simulation tool to obtain the time evolution of the system. This provides a good way of observing the effect of arbitrary pulses on the qubit dynamics and to check if the pulses are performing the desired operation. The next two cells run the simulations, extract the results, and plot them.

probabilities = {}
gate_times = {}
for scheme_name, control in gate_schemes.items():
    gate_durations, gate_states = simulations(control)
    probabilities[scheme_name] = np.abs(np.array(gate_states))**2
    gate_times[scheme_name] = gate_durations
100%|██████████| 2/2 [00:09<00:00,  4.99s/it, running=0]
100%|██████████| 2/2 [00:09<00:00,  4.63s/it, running=0]
100%|██████████| 2/2 [00:04<00:00,  2.41s/it, running=0]
gs = gridspec.GridSpec(1, 3, wspace=0.2)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(16)
fig.suptitle("Time evolution under control pulses", fontsize=16, y=1.15)
idx = 0
for scheme_name, gate in gate_schemes.items():
    ax = fig.add_subplot(gs[idx])
    ax.plot(gate_times[scheme_name],
            probabilities[scheme_name][:, 0], color="#BF04DC",
            label="$P_0$: Q-CTRL simulation")
    ax.plot(gate_times[scheme_name],
            probabilities[scheme_name][:, 1], color="#680CE9",
            label="$P_1$: Q-CTRL simulation")
    ax.set_title(scheme_name, fontsize=14)
    ax.set_xlabel("Time [s]", fontsize=14)
    ax.set_ylabel("Populations", fontsize=14)
    ax.tick_params(labelsize=14)
    idx += 1
ax.legend(loc="best", bbox_to_anchor=(-0.4, 1.3))
plt.show()

Ground ($P_0$) and excited ($P_1$) state population dynamics for the duration of the robust and primitive pulses. The large variation of the unconstrained pulse shape causes the rugged evolution of the populations in the left plot.

Return to Table of Contents

Robustness verification with filter functions and quasi-static scan

Filter functions provide a useful way to determine the sensitivity of pulses to different time varying noise mechanisms. In the next cells, we generate the filter functions for the Q-CTRL robust pulses and the primitive square pulse under dephasing noise.

sample_count = 4096
frequencies = np.logspace(-2, np.log10(omega_max), 2000)
schemes = {scheme: {} for scheme in gate_schemes.keys()}

for scheme_name, control in gate_schemes.items():
    scheme_objects = schemes[scheme_name]
    durations = [segment["duration"] for segment in control["I"]]
    I_values = [segment["value"] for segment in control["I"]]
    Q_values = [segment["value"] for segment in control["Q"]]

    system = qctrl.factories.systems.create(
        name=scheme_name, hilbert_space_dimension=2)

    shift_I = qctrl.factories.shift_controls.create(
        name="I", system=system, operator=sigma_x / 2.0)
    shift_I_pulse = qctrl.factories.custom_pulses.create(
        control=shift_I,
        segments=[{"duration": d, "value": v} for d, v in zip(durations, I_values)])

    shift_Q = qctrl.factories.shift_controls.create(
        name="Q", system=system, operator=sigma_y / 2.0)
    shift_Q_pulse = qctrl.factories.custom_pulses.create(
        control=shift_Q,
        segments=[{"duration": d, "value": v} for d, v in zip(durations, Q_values)])

    dephasing_noise = qctrl.factories.additive_noises.create(
        name="Dephasing", system=system, operator=sigma_z)

    scheme_objects["system"] = system
    scheme_objects["dephasing_noise"] = dephasing_noise


# Create filter function objects and calculate
for scheme_objects in schemes.values():
    scheme_objects[
        "dephasing_filter_function"] = qctrl.factories.filter_functions.create(
        noise=scheme_objects["dephasing_noise"],
        sample_count=sample_count,
        interpolated_frequencies=frequencies)
    scheme_objects[
        "calculated_dephasing_filter_function"] = qctrl.services.filter_functions.calculate(
        scheme_objects["dephasing_filter_function"])

filter_functions = {
    scheme: scheme_objects["calculated_dephasing_filter_function"].interpolated_points
    for scheme, scheme_objects in schemes.items()}
100%|██████████| 4/4 [00:24<00:00,  6.02s/it, running=0]
100%|██████████| 4/4 [00:22<00:00,  5.59s/it, running=0]
100%|██████████| 4/4 [00:26<00:00,  6.55s/it, running=0]

We also verify the pulse robustness against quasi-static noises by comparing the optimized $\pi$-pulses with the primitive pulse while detuning the drive frequency from the qubit frequency: $\Delta = \omega_\text{drive}-\omega_\text{qubit}$. This is tantamount to an effective dephasing term in the Hamiltonian, the strength and effect of which may be measured by scanning the detuning $\Delta$. The scan is calculated in the cell below.

gate_infidelity = {}
scan_array = np.linspace(-np.pi * 1.0 * mega, np.pi * 1.0 * mega, 21)
gate_infidelity["dephasing"] = {
    scheme_name: get_detuning_scan(scan_array, control)
    for scheme_name, control in gate_schemes.items()}
100%|██████████| 2/2 [00:07<00:00,  3.62s/it, running=0]
100%|██████████| 2/2 [00:09<00:00,  4.84s/it, running=0]
100%|██████████| 2/2 [00:07<00:00,  3.76s/it, running=0]

Here we plot the filter functions and the detuning scan

gs = gridspec.GridSpec(1, 2, wspace=0.2)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(16)
fig.suptitle("Robustness verification", fontsize=16, y=1.15)

ax = fig.add_subplot(gs[0])
ax.set_xlabel("Frequency (Hz)", fontsize=14)
ax.set_ylabel("Inverse power", fontsize=14)
ax.set_title("Filter functions", fontsize=14)
ax.tick_params(labelsize=14)
ax.set_xlim(np.min(frequencies), np.max(frequencies))
for scheme_name, _data in filter_functions.items():
    inverse_powers = np.array([_d["inverse_power"] for _d in _data])
    inverse_power_precisions = np.array([_d["inverse_power_precision"] for _d in _data])

    plt.loglog(frequencies, inverse_powers, "-", linewidth=2, 
               label=scheme_name, color=colors[scheme_name])
    y_upper = inverse_powers + inverse_power_precisions
    y_lower = inverse_powers - inverse_power_precisions
    ax.fill_between(frequencies, y_lower, y_upper,
                    hatch="||", alpha=0.35, facecolor="none", 
                    edgecolor=colors[scheme_name], linewidth=0.0)

ax = fig.add_subplot(gs[1])
ax.set_title("Robustness scan", fontsize=14)
ax.set_ylabel("Ground state population", fontsize=14)
ax.set_xlabel("Detuning (MHz)", fontsize=14)
for name, infidelity in gate_infidelity["dephasing"].items():
    ax.plot(scan_array * 1e-6 / np.pi, infidelity, label=name, color=colors[name])
ax.legend(loc="best", bbox_to_anchor=(0.0, 1.25))
plt.show()

Left plot: Filter functions for the three control pulses under dephasing noise. The hallmark of robustness is the significant drop of the filter function values at low frequencies. This feature is clearly observed for the optimized Q-CTRL pulses. Right plot: theoretical ground state population as a function of the detuning. The two Q-CTRL solutions show a flat behavior for a large range of detunings as compared to the primitive pulse, indicating robustness against dephasing.

Return to Table of Contents

2. Experimental implementation on IBM-Q

Before implementing the Q-CTRL pulses on IBM-Q, it is necessary to perform a calibration of the pulse amplitudes on the hardware.

Pulse calibration: mapping $(I, Q)$ values onto hardware amplitude inputs

The Q-CTRL optimization module returns a piecewise drive pulse, given in terms of a complex-valued waveform

\begin{align*} \gamma(t) & = I(t) + i Q(t)\\ &= \Omega(t) e^{i\phi(t)} \end{align*}

as a function of time $t$. The IQ variables relate to the (on-resonance) Rabi rate $\Omega$ and drive phase $\phi$ as

\begin{align*} I = \Omega \cos(\phi)\\ Q = \Omega \sin(\phi)\\ |\gamma|\le\Omega_\text{max} \end{align*}

On the IBM hardware, pulses are defined in terms of dimensionless variables $A = A_I + i A_Q$, with real and imaginary components proportional to physical input voltages, and bounded as

\begin{align*} A_I \in [-1,1] \\ A_Q \in [-1,1] \\ |A|\le 1. \end{align*}

The goal of this calibration routine is establish the map

\begin{align*} (A_I, A_Q) \longleftrightarrow (I, Q) \end{align*}

from hardware amplitudes to control pulse amplitudes. The calibration steps below assume that $I$ and $Q$ can be independently calibrated. The calibration steps are:

  1. Pick a given hardware input amplitude $A_I$ and create a constant rectangular pulse to act on the system.

  2. With the pulse from (1), perform a Rabi experiment by varying the pulse duration and measuring the qubit population. This will show Rabi oscillations over time, which we use to fit a cosine function and find the corresponding Rabi rate. Varying the pulse duration has the benefit of avoiding possible nonlinearities that can distort fitting of Rabi oscillations, encountered when scanning the amplitudes to observe Rabi oscillations.

  3. Repeat steps (1) and (2) with different values of the amplitude $A_I$ and make a table with the correspondence between $A_I$ and $I$. This correspondence will not necessarily be linear. Using this data set, we create a interpolating function that connects the $I$ values to the $A_I$ amplitudes. Performing the same procedure driving the $Q$ component allows us to obtain the $Q \leftrightarrow A_Q$ correspondence.

  4. The result from (3) can be used to translate the $(I, Q)$ values in an arbitrary Q-CTRL pulse into the appropriate $(A_I, A_Q)$ hardware input amplitudes.

In this setup cell we create a dictionary with the Rabi oscillation programs to run on the machine for each pulse amplitude.

# Get qubit frequency
qubit = 0
backend.properties(refresh=True)
qubit_frequency_updated = backend.properties().qubit_property(qubit, 'frequency')[0]

# Drive and measurement channels
drive_chan = DriveChannel(qubit)
meas_chan = MeasureChannel(qubit)
inst_sched_map = backend_defaults.instruction_schedule_map
measure_schedule = inst_sched_map.get('measure', qubits=backend_config.meas_map[qubit])

# Pulse times are picked so that roughly the same number of Rabi oscillations is observed for the
# different pulse amplitudes
# The minimum duration of a pulse in the armonk backend is 64*dt
num_time_points = 40
pulse_amp_array = np.linspace(0.1, 1, 10)
pulse_times = np.array([
    64 + np.arange(0,get_closest_multiple_of_16(32 * 10 * num_time_points / (idx + 1)),
                   get_closest_multiple_of_16(32 * 10 / (idx + 1)))
    for idx in np.arange(pulse_amp_array.shape[0])])

pulse_amp_array = np.concatenate((-pulse_amp_array[::-1], pulse_amp_array))
pulse_times = np.concatenate((pulse_times[::-1], pulse_times))

num_shots_per_point = 1024
rabi_programs_dic_I = {}

for idx, pulse_amp in enumerate(pulse_amp_array):
    rabi_schedules_I = []
    for integer_duration in pulse_times[idx]:
        waveform_I = np.ones([integer_duration]) * pulse_amp
        this_schedule = pulse.Schedule(name=f"I pulse: duration = {integer_duration}")
        this_schedule += pulse.SamplePulse(waveform_I)(drive_chan)
        this_schedule += measure_schedule << this_schedule.duration
        rabi_schedules_I.append(this_schedule)
    rabi_experiment_program_I = assemble(rabi_schedules_I, backend=backend,
                                         meas_level=1, meas_return="avg",
                                         shots=num_shots_per_point,
                                         schedule_los=[{drive_chan: qubit_frequency_updated}]
                                         * len(pulse_times[idx]))

    rabi_programs_dic_I[str(pulse_amp)] = rabi_experiment_program_I

Now we run the $I$ calibration with input amplitudes varying in the range [-1,1]. For this backend, we've run this calibration multiple times for both $I$ and $Q$ and they always coincide. For this reason, here we only run the calibration for $I$ and use it to calibrate both $I$ and $Q$ values. Note that in our fit we use the function $\cos(2\pi \, \Omega \, t +\phi)$ and the fitted rates will be given in Hz, not rad/s. It is important to keep this in mind when using the Q-CTRL pulse solutions.

if use_saved_data == False:
    rabi_calibration_exp_I = []
    rabi_oscillations_results = []
    for idx, pulse_amp in enumerate(pulse_amp_array):
        job = backend.run(rabi_programs_dic_I[str(pulse_amp)])
        job_monitor(job)
        rabi_results = job.result(timeout=120)
        rabi_values = []
        time_array = pulse_times[idx] * dt

        for i in range(len(pulse_times[idx])):
            rabi_values.append(rabi_results.get_memory(i)[qubit])
        rabi_values = get_amplitude(rabi_values)
        rabi_oscillations_results.append(rabi_values)
        A0 = np.max(rabi_values)
        fit_params, y_fit = fit_function(
            time_array, rabi_values,
            lambda x, A, rabi_freq, phi: A * np.cos(2 * np.pi * rabi_freq * x + phi),
            [A0, pulse_amp * 1.1e7, 0])  # initial parameters for curve_fit
        rabi_calibration_exp_I.append(fit_params[1])
        print("Drive amplitude:", pulse_amp)
        print("Fitted Rabi frequency [Hz]:", fit_params[1])
    
    plt.title("Examplary Rabi oscillation data with fitting")
    plt.xlabel("Time [s]")
    plt.ylabel("Measured signal [a.u.]")
    plt.scatter(time_array, rabi_values, color="black")
    plt.xlim(0, time_array[-1])
    plt.plot(time_array, y_fit, color="red")
    plt.show()
        
    filename = "rabi_calibration_exp_I_" + timestr
    outfile = open(filename, "wb")
    pickle.dump(rabi_calibration_exp_I, outfile)
    outfile.close()
else:
    filename = data_folder / "rabi_calibration_exp_I_demo"
    infile = open(filename, "rb")
    rabi_calibration_exp_I = pickle.load(infile)
    infile.close()
amplitude_interpolated_list = np.linspace(-1, 1, 100)

f_amp_to_rabi = interpolate.interp1d(pulse_amp_array, rabi_calibration_exp_I)
rabi_interpolated_exp_I = f_amp_to_rabi(amplitude_interpolated_list)

plt.title("IBM Armonk Rabi rates calibration", fontsize=16, y=1.05)
plt.xlabel("Hardware input amplitude", fontsize=14)
plt.ylabel("Rabi rate [Hz]", fontsize=14)
plt.scatter(pulse_amp_array, rabi_calibration_exp_I, label="$I$-calibration")
plt.plot(amplitude_interpolated_list, rabi_interpolated_exp_I)
plt.axvline(0, color="black", linestyle="dashed")
plt.axhline(0, color="black", linestyle="dashed")
plt.legend()
plt.show()

The plot shows the interpolation results from the Rabi oscillation fitting data. The curve provides the relationship between the hardware dimensionless input amplitudes and the Rabi rate ($I$ values) in Hz. Note that the relationship is non-linear for larger values of the amplitude.

Return to Table of Contents

Experimental time evolution of the driven qubit

In this section we will compare the implementation and performance of the pulses discussed previously. A good test to check that the Q-CTRL pulses are being implemented properly in the IBM hardware is to compare the theoretical with the experimental time evolution of the qubit state. The next cell executes the different steps to obtain such a comparison:

  • Perform pulse calibration on both the optimized and primitive pulses using the interpolation obtained in the previous section to establish the $\{I(t), Q(t)\}\rightarrow \{A_I (t), A_Q(t)\}$ conversion.
  • Setup the schedules and run them on IBM-Q
  • Extract and plot the resulting dynamics
# Calibration of pulses
f_rabi_to_amp = interpolate.interp1d(
    rabi_interpolated_exp_I, amplitude_interpolated_list)
waveform = {}
for scheme_name, control in gate_schemes.items():
    I_values = [segment["value"] / 2 / np.pi for segment in control["I"]]
    Q_values = [segment["value"] / 2 / np.pi for segment in control["Q"]]
    if scheme_name == "Square":
        A_I_values = np.repeat(f_rabi_to_amp(I_values), 
                               segment_scale[scheme_name] * number_of_segments[scheme_name])
        A_Q_values = np.repeat(f_rabi_to_amp(Q_values), 
                               segment_scale[scheme_name] * number_of_segments[scheme_name])
    else:
        A_I_values = np.repeat(f_rabi_to_amp(I_values), segment_scale[scheme_name])
        A_Q_values = np.repeat(f_rabi_to_amp(Q_values), segment_scale[scheme_name])
    waveform[scheme_name] = A_I_values + 1j * A_Q_values

pulse_evolution_program = {}
for scheme_name in gate_schemes.keys():
    time_min = 64
    time_max = duration_int[scheme_name]
    time_step = 32
    num_shots_per_point = 2048
    times_int = np.arange(time_min, time_max + time_step, time_step)
    ibm_evolution_times = times_int * dt
    num_time_points = len(times_int)

    evolution_schedules = []
    for time_idx in times_int:
        schedule = pulse.Schedule(name="duration_%d" % time_idx)
        schedule = schedule | pulse.SamplePulse(waveform[scheme_name][:time_idx])(
            drive_chan)
        schedule = schedule | measure_schedule << schedule.duration
        evolution_schedules.append(schedule)

    pulse_evolution_program[scheme_name] = assemble(
        evolution_schedules,
        backend=backend,
        meas_level=2,
        meas_return="single",
        shots=num_shots_per_point,
        schedule_los=[{drive_chan: qubit_frequency_updated}] * num_time_points)
    
# Running the jobs
if use_saved_data == False:
    evolution_exp_results_dephasing = {}
    for scheme_name in gate_schemes.keys():
        job = backend.run(pulse_evolution_program[scheme_name])
        job_monitor(job)
        evolution_exp_results_dephasing[scheme_name] = job.result(timeout=120)
    filename = "evolution_exp_results_dephasing_" + timestr
    outfile = open(filename, "wb")
    pickle.dump(evolution_exp_results_dephasing, outfile)
    outfile.close()
else:
    filename = data_folder / "evolution_exp_results_demo_dephasing"
    infile = open(filename, "rb")
    evolution_exp_results_dephasing = pickle.load(infile)
    infile.close()

# Extract results and plot
evolution_ibm = {}
for scheme_name in gate_schemes.keys():
    evolution_exp_data = np.zeros(num_time_points)
    for idx, time_idx in enumerate(times_int):
        counts = evolution_exp_results_dephasing[scheme_name].get_counts(
            "duration_%d" % time_idx)
        excited_pop = 0
        for bits, count in counts.items():
            excited_pop += count if bits[::-1][qubit] == "1" else 0
        evolution_exp_data[idx] = excited_pop / num_shots_per_point
    evolution_ibm[scheme_name] = evolution_exp_data

gs = gridspec.GridSpec(1, 3, wspace=0.2)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(16)
fig.suptitle("Time evolution under control pulses", fontsize=16, y=1.05)
idx = 0
for scheme_name in gate_schemes.keys():
    ax = fig.add_subplot(gs[idx])
    ax.plot(gate_times[scheme_name], probabilities[scheme_name][:, 0],
            color=colors["Unconstrained"], label="$P_0$: Q-CTRL simulation")
    ax.plot(gate_times[scheme_name], probabilities[scheme_name][:, 1],
            color=colors["Bound slew"], label="$P_1$: Q-CTRL simulation")
    ax.scatter(ibm_evolution_times, 1 - evolution_ibm[scheme_name],
               color=colors["Unconstrained"], label="$P_0$: IBM experiments")
    ax.scatter(ibm_evolution_times, evolution_ibm[scheme_name],
               color=colors["Bound slew"], label="$P_1$: IBM experiments")
    ax.set_title(scheme_name, fontsize=14)
    ax.set_xlabel("Time [s]", fontsize=14)
    ax.set_ylabel("Populations", fontsize=14)
    ax.tick_params(labelsize=14)
    idx += 1
ax.legend(loc="best", bbox_to_anchor=(1, 1.0))
plt.show()

When it comes to implementing a time-varying pulse, any hardware will have bandwidth limitations. The abrupt variations on the waveform of our Unconstrained control pulse seems to be beyond what can be reasonably implemented in the hardware, causing the deviations between IBM experiments and Q-CTRL simulations observed on the left plot. This highlights the importance of developing hardware-tailored solutions, as in the Bound slew pulse. Constraining how much the pulse amplitude can vary from segment to segment delivers a smoother control pulse that is successfully implemented on the IBM hardware, as shown in the center plot.

Return to Table of Contents

Robustness verification with a quasi-static scan

To demonstrate robustness, here we perform the experimental detuning scan, similar to what was done for the theoretical controls. As an extra comparison case, we have also added to this scan the default X-gate for this IBM backend. The steps to produce achieve this, all combined in the next execution cell, are the following:

  • Define detuning sweep interval
  • Create schedules and running jobs on IBM-Q
  • Extracting results and produce robustness plot
scheme_names.append('IBM default X-gate')
# Define detuning sweep interval
detuning_array = np.linspace(
    qubit_frequency_updated - 1e6, qubit_frequency_updated + 1e6, 21)
drive_frequencies = [{drive_chan: freq} for freq in detuning_array]
number_of_runs = 20
num_shots_per_point = 512

# Create schedules and run jobs
if use_saved_data == False:
    pi_experiment_results = {}
    for run in range(number_of_runs):
        for scheme_name in scheme_names:
            # Default IBM X-gate
            if scheme_name=='IBM default X-gate':
                pi_schedule = pulse.Schedule(name="IBM default X-gate")
                pi_schedule |= inst_sched_map.get('u3', [qubit], P0=np.pi, P1=0.0, P2=np.pi)
            else:
                pi_schedule = pulse.Schedule(name="Q-CTRL robust pulse")
                pi_schedule += pulse.SamplePulse(waveform[scheme_name])(drive_chan)
            pi_schedule += measure_schedule << pi_schedule.duration

            pi_experiment_ibm = assemble(
                pi_schedule,
                backend=backend,
                meas_level=2,
                meas_return="single",
                shots=num_shots_per_point,
                schedule_los=drive_frequencies)
            job = backend.run(pi_experiment_ibm)
            print(job.job_id())
            job_monitor(job)
            pi_experiment_results[scheme_name + str(run)] = job.result(timeout=120)

# Extract results
if use_saved_data == False:
    detuning_sweep_results = {}
    for scheme_name in scheme_names:
        qctrl_sweep_values = np.zeros((number_of_runs, len(detuning_array)))
        for run in range(number_of_runs):
            i = 0
            for result in pi_experiment_results[scheme_name + str(run)].results:
                counts = result.data.counts.to_dict()["0x1"]
                excited_pop = counts / num_shots_per_point
                qctrl_sweep_values[run][i] = 1 - excited_pop
                i = i + 1
        detuning_sweep_results[scheme_name] = qctrl_sweep_values
    filename = "detuning_sweep_results_" + timestr
    outfile = open(filename, "wb")
    pickle.dump(detuning_sweep_results, outfile)
    outfile.close()
else:
    filename = data_folder / "detuning_sweep_results_demo"
    infile = open(filename, "rb")
    detuning_sweep_results = pickle.load(infile)
    infile.close()
# Robustness plot
fig, ax = plt.subplots()
idx = 0
for scheme_name, probability in detuning_sweep_results.items():
    if scheme_name == "Bound slew" or scheme_name == "Square":
        lw = 2
        ms = 5
        el = 1
    else:
        lw = 1
        ms = 2.5
        el = 0.5
    ax.errorbar(detuning_array / mega - qubit_frequency_updated / mega,
                np.mean(probability, axis=0), yerr=np.std(probability, axis=0)
                / np.sqrt(len(detuning_sweep_results[scheme_name])), fmt="-o",
                linewidth=lw, markersize=ms, elinewidth=el, capsize=3,
                label=scheme_name + " run on IBM-Q", color=colors[scheme_name])
    if scheme_name != "IBM default X-gate":
        plt.plot(detuning_array / mega - qubit_frequency_updated / mega,
                 gate_infidelity["dephasing"][scheme_name] + 0.11, linestyle="dashed",
                 label=scheme_name + str(": simulations"), color=colors[scheme_name])
    #idx += 1
plt.xlabel("Detuning [MHz]", fontsize=14)
plt.ylabel("Ground state population ", fontsize=14)
plt.title("Detuning scan", fontsize=16, y=1.05)
plt.legend(loc="best", bbox_to_anchor=(1, 1.0))
plt.show()

Experimental (symbols) and simulations (lines) ground state populations as a function of the detuning. The Bound slew control solution is mostly flat across the detuning scan, showing its robustness against dephasing as compared to the primitive pulse. The irregular behavior for the Unconstrained solution is a consequence of the imperfect implementation of this pulse, as already seen in the simulation plots. The theoretical curves have been shifted up by $0.11$ to account for measurement errors in the experimental data. Note the drift in the square pulse curve: the calibration has gone off from the time it has been performed and the time this robustness scan was generated, indicating the benefit of the robust controls.

Return to Table of Contents

3. Other noise sources: Q-CTRL pulses robust against control error

Besides robustness to dephasing, the Q-CTRL Python package can also be used to design solutions that are robust against errors in the control amplitudes. The workflow is exactly as in the dephasing case:

  • Initialization cell defining Q-CTRL optimization parameters
  • Run optimization and extract results
  • Characterize robustness using filter functions
  • Pulse calibration and implementation
  • Robustness verification with amplitude scan

Creating robust pulses

The next cell contains all Q-CTRL commands to generate pulses for the unconstrained and bound slew control scenarios defined previously. The only difference to the dephasing case is the setting of noise=True for the $I$ and $Q$ controls (shift terms) and removing the dephasing noise (drift term) when calculating the optimized controls.

number_of_segments = 128
segment_scale = 20
duration_int = number_of_segments * segment_scale
duration = duration_int * dt

shift_unconstrained = {'I': {'operator': sigma_x/2.,
                             'segments_count': number_of_segments,
                             'amplitude_upper_bound': I_max,
                             'noise': True},
                       'Q': {'operator': sigma_y/2.,
                             'segments_count': number_of_segments,
                             'amplitude_upper_bound': Q_max,
                             'noise': True}}

shift_bound_slew = {'I': {'operator': sigma_x/2.,
                          'segments_count': number_of_segments,
                          'amplitude_upper_bound': I_max,
                          'maximum_slew_rate':I_max/5,
                          'noise': True},
                    'Q': {'operator': sigma_y/2.,
                          'segments_count': number_of_segments,
                          'amplitude_upper_bound': Q_max,
                          'maximum_slew_rate':Q_max/5,
                          'noise': True}}

shift_list = [shift_unconstrained, shift_bound_slew]
scheme_names = ["Unconstrained", "Bound slew"]

# Run optimization
if use_saved_data == False:
    robust_amplitude_controls = {}
    for idx, shift in enumerate(shift_list): 
        start_time = time.time() 
         
        # Create system  
        system = qctrl.factories.systems.create( 
            name="system " + str(idx), 
            hilbert_space_dimension=2) 
         
        # Create shifts 
        shift_I = qctrl.factories.shift_controls.create( 
            name="I", 
            system=system, 
            operator=shift['I']['operator']) 
         
        shift_Q = qctrl.factories.shift_controls.create(
            name="Q", 
            system=system, 
            operator=shift['Q']['operator']) 
         
        # Create pulses 
        pulse_I = qctrl.factories.optimum_pulses.create( 
            control=shift_I, 
            upper_bound=shift['I']['amplitude_upper_bound'], 
            fixed_modulus=False, 
            segment_count=shift['I']['segments_count'], 
            duration=duration, 
            maximum_slew_rate=shift['I'].get('maximum_slew_rate', None)) 
         
        pulse_Q = qctrl.factories.optimum_pulses.create( 
            control=shift_Q, 
            upper_bound=shift['Q']['amplitude_upper_bound'], 
            fixed_modulus=False, 
            segment_count=shift['Q']['segments_count'], 
            duration=duration, 
            maximum_slew_rate=shift['Q'].get('maximum_slew_rate', None)) 
         
        # Create noises 
        amplitude_noise = qctrl.factories.control_noises.create(
            name="Amplitude noise", control=shift_I)
         
        # Create target 
        target = qctrl.factories.targets.create(
            unitary_operator=X_gate, projection_operator=np.identity(2), system=system) 
         
        # Run optimization 
        system = qctrl.services.robust_optimization.run(system=system) 
         
        # Store results 
        robust_amplitude_controls[scheme_names[idx]] = {control.name: control.pulse.segments  
                                                     for control in system.controls if control.pulse} 
        print('Execution time for '+ str(number_of_runs)+ ' runs (in sec):', time.time()-start_time) 
        print(scheme_names[idx] + ' cost = ', system.cost, '\n') 
        

    filename = "robust_amplitude_controls_" + timestr
    outfile = open(filename, "wb")
    pickle.dump(robust_amplitude_controls, outfile)
    outfile.close()
else:
    filename = data_folder / "robust_amplitude_controls_demo"
    infile = open(filename, "rb")
    robust_amplitude_controls = pickle.load(infile)
    infile.close()

# Create square pulse
square_pulse_durations = np.array([duration])
square_pulse_values_control = np.array([rabi_rotation / duration])
square_sequence = { "I": [{"duration": d, "value": v}
                          for d, v in zip(square_pulse_durations, square_pulse_values_control)],
                   "Q": [{"duration": d, "value": v}
                         for d, v in zip(square_pulse_durations, np.zeros([number_of_segments]))]}

# Combine all pulses into a single dictionary
gate_schemes = {**robust_amplitude_controls, **{"Square": square_sequence}}

Here we plot the controls

for scheme_name, gate in gate_schemes.items():
    print(scheme_name)
    plot_controls(plt.figure(), gate)
    plt.show()
Unconstrained
Bound slew
Square

$I(t)$ and $Q(t)$ components of the control pulses for: unconstrained (top) and bound slew (middle) optimizations, and primitive control (bottom).

Return to Table of Contents

Time evolution of the driven qubit: simulations and experiment

As for the dephasing robust pulses, we calibrate the control error robust pulses using the interpolation obtained previously and run them on the IBM machine to compare the theoretical with the experimental time evolution of the qubit state under no added noise. The cell below combines the following steps:

  • Pulse calibration
  • Setup IBM schedules
  • Run jobs on IBM-Q and extract results
  • Run Q-CTRL simulations
  • Plot results
# Pulse calibration
f_rabi_to_amp = interpolate.interp1d(
    rabi_interpolated_exp_I, amplitude_interpolated_list)

waveform = {}
for scheme_name, control in gate_schemes.items():
    I_values = [segment["value"] / 2 / np.pi for segment in control["I"]]
    Q_values = [segment["value"] / 2 / np.pi for segment in control["Q"]]
    if scheme_name == "Square":
        A_I_values = np.repeat(
            f_rabi_to_amp(I_values), segment_scale * number_of_segments
        )
        A_Q_values = np.repeat(
            f_rabi_to_amp(Q_values), segment_scale * number_of_segments
        )
    else:
        A_I_values = np.repeat(f_rabi_to_amp(I_values), segment_scale)
        A_Q_values = np.repeat(f_rabi_to_amp(Q_values), segment_scale)
    waveform[scheme_name] = A_I_values + 1j * A_Q_values

time_min = 64
time_max = duration_int
time_step = 48
num_shots_per_point = 2048
times_int = np.arange(time_min, time_max + time_step, time_step)
ibm_evolution_times = times_int * dt
num_time_points = len(times_int)

pulse_evolution_program = {}
for scheme_name in gate_schemes.keys():
    evolution_schedules = []
    for time_idx in times_int:
        schedule = pulse.Schedule(name="duration_%d" % time_idx)
        schedule = schedule | pulse.SamplePulse(waveform[scheme_name][:time_idx])(drive_chan)
        schedule = schedule | measure_schedule << schedule.duration
        evolution_schedules.append(schedule)

    pulse_evolution_program[scheme_name] = assemble(
        evolution_schedules, backend=backend, meas_level=2, meas_return="single",
        shots=num_shots_per_point, schedule_los=[{drive_chan: qubit_frequency_updated}] * num_time_points)

# Run jobs
if use_saved_data == False:
    evolution_exp_results = {}
    for scheme_name in gate_schemes.keys():
        job = backend.run(pulse_evolution_program[scheme_name])
        job_monitor(job)
        evolution_exp_results[scheme_name] = job.result(timeout=120)
    filename = "evolution_exp_results_" + timestr
    outfile = open(filename, "wb")
    pickle.dump(evolution_exp_results, outfile)
    outfile.close()
else:
    filename = data_folder / "evolution_exp_results_demo_amplitude"
    infile = open(filename, "rb")
    evolution_exp_results = pickle.load(infile)
    infile.close()

# Extract results
evolution_ibm = {}
for scheme_name in gate_schemes.keys():
    evolution_exp_data = np.zeros(num_time_points)
    for idx, time_idx in enumerate(times_int):
        counts = evolution_exp_results[scheme_name].get_counts("duration_%d" % time_idx)
        excited_pop = 0
        for bits, count in counts.items():
            excited_pop += count if bits[::-1][qubit] == "1" else 0
        evolution_exp_data[idx] = excited_pop / num_shots_per_point
    evolution_ibm[scheme_name] = evolution_exp_data

# Q-CTRL simulation
probabilities = {}
for scheme_name, control in gate_schemes.items():
    gate_times, gate_states = simulations(control)
    probabilities[scheme_name] = np.abs(np.array(gate_states)) ** 2

# Plotting
gs = gridspec.GridSpec(1, 3, wspace=0.2)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(16)
fig.suptitle("Time evolution under control pulses", fontsize=16, y=1.1)
idx = 0
for scheme_name in gate_schemes.keys():
    ax = fig.add_subplot(gs[idx])
    ax.plot(gate_times, probabilities[scheme_name][:, 0],
            color=colors["Unconstrained"], label="$P_0$: Q-CTRL simulations")
    ax.plot(gate_times, probabilities[scheme_name][:, 1],
            color=colors["Bound slew"], label="$P_1$: Q-CTRL simulations")
    ax.scatter(ibm_evolution_times, 1 - evolution_ibm[scheme_name],
               color=colors["Unconstrained"], label="$P_0$: IBM experiments")
    ax.scatter(ibm_evolution_times, evolution_ibm[scheme_name],
               color=colors["Bound slew"], label="$P_1$: IBM experiments")
    ax.set_title(scheme_name, fontsize=14)
    ax.set_xlabel("Time [s]", fontsize=14)
    ax.set_ylabel("Populations", fontsize=14)
    ax.tick_params(labelsize=14)
    idx += 1
ax.legend(loc="best", bbox_to_anchor=(1, 1.0))
plt.show()
100%|██████████| 2/2 [00:08<00:00,  4.34s/it, running=0]
100%|██████████| 2/2 [00:07<00:00,  3.93s/it, running=0]
100%|██████████| 2/2 [00:08<00:00,  4.42s/it, running=0]

Time evolution of ground ($P_0$) and excited ($P_1$) state populations for simulations (lines) and experimental (symbols) runs under the control-robust pulses.

Return to Table of Contents

Robustness verification with filter functions and a quasi-static scan

We again determine the sensitivity of the pulses to the noise and compare them to the primitive using filter functions and static scans. The only difference here is the change in the noise process, which is achieved by omitting the dephasing term and setting the noise on the controls to be True.

For the filter functions calculated in the cell below, we use noise in the $I$ component of the control but similar filter functions can be obtained with noise in the $Q$ component.

sample_count = 4096
frequencies = np.logspace(-2, np.log10(omega_max), 2000)

schemes = {scheme: {} for scheme in gate_schemes.keys()}
for scheme_name, control in gate_schemes.items():
    scheme_objects = schemes[scheme_name]
    durations = [segment["duration"] for segment in control["I"]]
    I_values = [segment["value"] for segment in control["I"]]
    Q_values = [segment["value"] for segment in control["Q"]]
    system = qctrl.factories.systems.create(
        name=scheme_name, hilbert_space_dimension=2)
    shift_I = qctrl.factories.shift_controls.create(
        name="I", system=system, operator=sigma_x / 2.0)
    shift_pulse_I = qctrl.factories.custom_pulses.create(
        control=shift_I,
        segments=[{"duration": d, "value": v} for d, v in zip(durations, I_values)])
    shift_Q = qctrl.factories.shift_controls.create(
        name="Q", system=system, operator=sigma_y / 2.0)
    shift_pulse_Q = qctrl.factories.custom_pulses.create(
        control=shift_Q,
        segments=[{"duration": d, "value": v} for d, v in zip(durations, Q_values)])
    
    # Define noises
    amplitude_noise = qctrl.factories.control_noises.create(
        name="Amplitude noise", control=shift_I,
    )
    scheme_objects["system"] = system
    scheme_objects["amplitude_noise"] = amplitude_noise


# Create filter function objects and calculate
for scheme_objects in schemes.values():
    scheme_objects[
        "amplitude_filter_function"] = qctrl.factories.filter_functions.create(
        noise=scheme_objects["amplitude_noise"],
        sample_count=sample_count,
        interpolated_frequencies=frequencies)
    scheme_objects[
        "calculated_amplitude_filter_function"] = qctrl.services.filter_functions.calculate(
        scheme_objects["amplitude_filter_function"])

filter_functions = {
    scheme: scheme_objects["calculated_amplitude_filter_function"].interpolated_points
    for scheme, scheme_objects in schemes.items()}
100%|██████████| 4/4 [00:26<00:00,  6.58s/it, running=0]
100%|██████████| 4/4 [00:25<00:00,  6.50s/it, running=0]
100%|██████████| 4/4 [00:26<00:00,  6.71s/it, running=0]

We also demonstrate robustness by performing an control error scan, both numerically and experimentally. The steps follow the same structure as in the dephasing case and are combined in the next cell:

  • Define control error range
  • Create schedules and run them on IBM-Q
  • Create Q-CTRL theoretical control error scan
# Define control amplitude error range
amp_error = 0.25
amp_scalings = np.linspace(1 - amp_error, amp_error + 1, 21)

# Running jobs
if use_saved_data == False:
    pi_experiment_results = {}
    for run in range(number_of_runs):
        for scheme_name in gate_schemes.keys():
            pi_schedule = []
            for amp_scaling in amp_scalings:
                this_schedule = pulse.Schedule(name="Q-CTRL robust pulse")
                this_schedule += pulse.SamplePulse(waveform[scheme_name] 
                                                   * amp_scaling)(drive_chan) << 1
                this_schedule += measure << duration_int
                pi_schedule.append(this_schedule)
            num_shots_per_point = 512

            pi_experiment_ibm = assemble(
                pi_schedule,
                backend=backend,
                meas_level=2,
                meas_return="single",
                shots=num_shots_per_point,
                schedule_los=[{drive_chan: qubit_frequency_updated}]
                * len(amp_scalings))
            job = backend.run(pi_experiment_ibm)
            print(job.job_id())
            job_monitor(job)
            pi_experiment_results[scheme_name + str(run)] = job.result(timeout=120)

else:
    filename = data_folder / "amplitude_sweep_results_demo"
    infile = open(filename, "rb")
    amplitude_sweep_results = pickle.load(infile)
    infile.close()

# Extracting results
if use_saved_data == False:
    amplitude_sweep_results = {}
    for scheme_name in gate_schemes.keys():
        qctrl_sweep_values = np.zeros((number_of_runs, len(amp_scalings)))
        for run in range(number_of_runs):
            i = 0
            for result in pi_experiment_results[scheme_name + str(run)].results:
                counts = result.data.counts.to_dict()["0x1"]
                excited_pop = counts / num_shots_per_point
                qctrl_sweep_values[run][i] = 1 - excited_pop
                i = i + 1
        amplitude_sweep_results[scheme_name] = qctrl_sweep_values
    filename = "amplitude_sweep_results_" + timestr
    outfile = open(filename, "wb")
    pickle.dump(amplitude_sweep_results, outfile)
    outfile.close()
else:
    filename = data_folder / "amplitude_sweep_results_demo"
    infile = open(filename, "rb")
    amplitude_sweep_results = pickle.load(infile)
    infile.close()

# Q-CTRL control amplitude scan
gate_infidelity = {}
amplitude_array = np.linspace(-0.25, 0.25, 21)
gate_infidelity["amplitude"] = {
    scheme_name: get_amplitude_scan(amplitude_array, control)
    for scheme_name, control in gate_schemes.items()}
100%|██████████| 2/2 [00:07<00:00,  3.96s/it, running=0]
100%|██████████| 2/2 [00:08<00:00,  4.15s/it, running=0]
100%|██████████| 2/2 [00:04<00:00,  2.30s/it, running=0]

Here we plot the numerical filter functions and the control error scan

gs = gridspec.GridSpec(1, 2, wspace=0.2)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(16)
fig.suptitle("Robustness verification", fontsize=16, y=1.1)

ax = fig.add_subplot(gs[0])
ax.set_xlabel("Frequency (Hz)", fontsize=14)
ax.set_ylabel("Inverse power", fontsize=14)
ax.set_title("Filter functions", fontsize=14)
ax.tick_params(labelsize=14)
ax.set_xlim(np.min(frequencies), np.max(frequencies))

for scheme_name, _data in filter_functions.items():
    inverse_powers = np.array([_d["inverse_power"] for _d in _data])
    inverse_power_precisions = np.array([_d["inverse_power_precision"] for _d in _data])

    plt.loglog(frequencies, inverse_powers, "-", linewidth=2, 
               label=scheme_name, color=colors[scheme_name])
    y_upper = inverse_powers + inverse_power_precisions
    y_lower = inverse_powers - inverse_power_precisions
    ax.fill_between(frequencies, y_lower, y_upper, hatch="||", alpha=0.35,
                    facecolor="none", edgecolor=colors[scheme_name], linewidth=0.0)
ax.legend(loc="best")

ax = fig.add_subplot(gs[1])
ax.set_title("Robustness scan", fontsize=14)
ax.set_ylabel("Ground state population", fontsize=14)
ax.set_xlabel("Control amplitude error", fontsize=14)
for scheme_name, probability in amplitude_sweep_results.items():
    if scheme_name == "Bound slew" or scheme_name == "Square":
        lw = 2
        ms = 5
        el = 1
    else:
        lw = 1
        ms = 2.5
        el = 0.5
    ax.errorbar(amp_scalings - 1, np.mean(probability, axis=0),
                yerr=np.std(probability, axis=0)
                / np.sqrt(len(amplitude_sweep_results[scheme_name])),
                fmt="-o", linewidth=lw, markersize=ms, elinewidth=el,
                capsize=3, label=scheme_name + " run on IBM-Q", color=colors[scheme_name])
    plt.plot(amp_scalings - 1, gate_infidelity["amplitude"][scheme_name] + 0.13,
             linestyle="dashed", label=scheme_name + str(": simulations"), color=colors[scheme_name])
plt.xlabel("Control amplitude Error", fontsize=14)
plt.ylabel("Ground state population ", fontsize=14)
plt.title("Control error scan", fontsize=16, y=1.05)
ax.legend(loc="best")
plt.show()

Left plot: Filter functions for the three control pulses under control noise in the $I$ component. Again, robustness is clearly observed for the optimized Q-CTRL pulses. Right plot: experimental (symbols) and simulation (lines) ground state populations as a function of the error in the control amplitude (note that the control amplitudes vary between -1 and 1). The Bound slew control solution is robust, with a flat response for a large range of control errors as compared to the primitive pulse. The theoretical curves have been shifted up by $0.13$ to account for measurement errors in the experimental data.

4. Q-CTRL solutions to complex problems on the superconducting qubit architecture

The Q-CTRL tools illustrated above can also be used to solve more complex problems in superconducting qubit hardware. Below we provide links to other notebooks where some of these issues are addressed:

  1. Leakage suppression
    In the linked notebook, we deal with the fact that superconducting qubits are actually oscillators and their description as a two-level is an approximation. This means that, during operations, leakage to other levels can occur. Considering the system as an oscillator truncated to three levels (qutrit), we can use the optimization tools in the Q-CTRL Python package to perform a desired gate with high fidelity while, at the same time, minimizing leakage. The notebook illustrates this for the implementation of the iSWAP gate between two qutrits.

  2. Cross-resonance gates
    Performing high-fidelity two-qubit gates is an essential step in any quantum computing architecture. In this notebook we present robust solutions to the implementation of cross-resonance gates with superconducting qubits.

  3. $T_1$ effects on optimized pulses
    Through independent simulations of a decaying proccess, this notebook analyses the effect of $T_1$ on Q-CTRL pulses.

Wiki

Comprehensive knowledge base of quantum control theory

Explore