How to calculate and use filter functions for arbitrary controls

Calculate the frequency-domain noise sensitivity of driven controls

The Q-CTRL Python package allows the calculation of filter functions as a way of estimating the sensitivity of an arbitrary control on a D-dimensional Hilbert space to time-dependent noise expressed in the frequency domain. In this guide we show how to define, calculate, and visualize filter functions obtained using the Q-CTRL Python package.

For more general information on the mathematical underpinning of the filter function framework see Section III.B. of Software tools for quantum control: Improving quantum computer performance through noise and error suppression

Overview of filter function workflows

Usually, workflows involving filter function calculations include these four steps. A concrete example of these is in the next section.

1. Calculate filter functions

The function qctrl.functions.calculate_filter_function takes the following parameters:

  • duration, the duration of the control pulses,
  • frequencies, which is a list of frequencies where the filter function is to be sampled,
  • drives, shifts, and drifts, which represent the different terms in the Hamiltonian (at least one term needs to be provided),
  • sample_count, which can be used to adjust the precision of the calculation of the filter functions (setting it to None uses the exact filter function calculation),
  • projection_operator, which projects into the subspace of interest, and
  • result_scope, which specifies whether to return the frequency domain noise operators.

More specifically, the three different types of Hamiltonian terms correspond to:

2. Access the frequency domain noise operators

The frequency domain noise operators are computed as part of the filter function calculation. They are available in the result samples when the result_scope parameter is set appropriately.

3. Visualize filter functions

After their calculation, the numerical values of the filter functions are stored at filter_function_result.samples (where filter_function_result is the object returned by qctrl.functions.calculate_filter_function). Each sample contains a corresponding frequency (our x coordinates), and an inverse_power (our y coordinates).

This information can be plotted using the plot_filter_functions function from the Q-CTRL Python Visualizer package.

The resulting graph represents the "noise admittance" of the control - where it is low, the noise is suppressed and where it is high, the noise is transmitted.

4. Calculate control infidelity due to noise with filter functions

As explained in the reference documentation for qctrl.functions.calculate_filter_function, the overlap integrals of the filter functions with the noise power spectral densities (PSDs) can be used to estimate the operational infidelity in the weak noise regime. If multiple sources of noise are present, the total infidelity $\mathcal{I}$ can be estimated by adding the contributions of each noise as represented by each filter function $F_k(f)$ and power spectral density $S_k(f)$ where $k$ is the index for the list of filter functions:

$$ \mathcal{I} \approx \sum_k \int_{-\infty}^\infty F_k (f) S_k (f) \,\, \mathrm{d}f. $$

With the filter functions we calculated previously, we can find the operational infidelity with respect to power spectral densities that we define. To do this, we just need to numerically calculate the overlap integrals, which can be done using the np.trapz NumPy function for composite trapezoidal rule integration. In general this first-order approximation provides a good quantitative estimate for noise-infidelities $\lesssim 10\%$. At higher noise strengths simple extensions may be used to approximate higher-order contributions (see Experimental noise filtering by quantum control, Supplementary information).

Note that if the pulses were not capable of perfectly creating the target X gates in the absence of noise, we should add the value of the noise-free infidelity, $\mathcal{I}_0$, to the calculation of the total infidelity: $ \mathcal{I} \approx \mathcal{I}_0 + \sum_k \int_{-\infty}^\infty F_k (f) S_k (f) \, \mathrm{d}f$. In the case of pulses that were predefined to generate the desired gate, $\mathcal{I}_0 \approx 0$.

Worked example: Noise sensitivity of composite pulses applied to a single qubit

In this example, we will use filter functions to compare the sensitivity of different composite $\pi$-pulses to time-varying amplitude and dephasing noise. The Hamiltonian of the system we will be considering is:

$$ H(t) = \frac{1 + \beta_\Omega(t)}{2} \left[ \Omega(t) \sigma_- + \Omega^* (t) \sigma_+ \right] + \frac{\Delta(t)}{2} \sigma_z + \frac{\eta(t)}{2} \sigma_z, $$

where $\Omega(t)$ is a time-dependent Rabi rate, $\beta_\Omega(t)$ is a fractional time-dependent amplitude fluctuation process, $\Delta(t)$ is a time-dependent clock shift, $\eta(t)$ is a small, slowly-varying stochastic dephasing noise process, and $\sigma_k$ are the Pauli matrices (with $\sigma_\pm = (\sigma_x \mp i \sigma_y)/2$).

Because there are two different noise channels we will be calculating two different filter functions for each control we study. This will permit a direct comparison of how each control exhibits different responses to different noise channels. We will compare the filter functions of each for a range of noise frequencies from $10^{-8} \Omega_\mathrm{max}$ to $\Omega_\mathrm{max}$, where $\Omega_\mathrm{max}/2\pi = 1~\mathrm{MHz}$ is the maximum Rabi frequency.

Create the pulses

We will consider the following driven control schemes for the controllable $\Omega (t)$ and $\Delta (t)$ terms:

  • Primitive: sensitive to both dephasing and amplitude noise
  • BB1: suppresses amplitude noise
  • CORPSE: suppresses dephasing noise
  • CORPSE in BB1: suppresses both dephasing and amplitude noise

These schemes are available from Q-CTRL Open Controls and described in the reference documentation.

We first set up Python objects representing the pulse, control, and noise. In the Hamiltonian used in this example, the drive corresponds to the term that contains the Rabi rate ($\Omega(t)$ and $\Omega^*(t)$), the shift corresponds to the term that contains the clock shift ($\Delta(t)$), and the drift corresponds to the dephasing term.

The Rabi coupling term and the clock-shift term are defined using a pulse from Q-CTRL Open Controls.

As we are interested in separately studying the robustness against amplitude noise in the complex-valued controls and in the dephasing term, we define two versions of the drive: one with noise, and one without it.

import matplotlib.pyplot as plt
import numpy as np
from attr import asdict
from dataclasses import dataclass
from typing import List

# Predefined pulse imports
from qctrlopencontrols import (
    new_bb1_control,
    new_corpse_control,
    new_corpse_in_bb1_control,
    new_primitive_control,
)
from qctrlopencontrols.driven_controls.driven_control import DrivenControl
from qctrlvisualizer import plot_filter_functions

from qctrl import Qctrl

# Starting a session with the API
qctrl = Qctrl()
# Define standard matrices
identity = np.array([[1, 0], [0, 1]], dtype=complex)
sigma_z = np.array([[1, 0], [0, -1]], dtype=complex)
sigma_x = np.array([[0, 1], [1, 0]], dtype=complex)
sigma_m = np.array([[0, 1], [0, 0]], dtype=complex)

# Define control parameters
omega_max = 2 * np.pi * 1e6  # Hz
total_rotation = np.pi

# Define schemes for driven controls to compare
pulses: List[DrivenControl] = [
    pulse_function(
        rabi_rotation=total_rotation,
        azimuthal_angle=0.0,
        maximum_rabi_rate=omega_max,
        name=name,
    )
    for name, pulse_function in [
        ("primitive", new_primitive_control),
        ("BB1", new_bb1_control),
        ("CORPSE", new_corpse_control),
        ("CORPSE in BB1", new_corpse_in_bb1_control),
    ]
]


@dataclass
class ControlScheme:
    name: str
    amplitude_filter_function_samples: qctrl.types.filter_function.Sample
    dephasing_filter_function_samples: qctrl.types.filter_function.Sample


def create_control_scheme(
    pulse: DrivenControl, frequencies: np.ndarray
) -> ControlScheme:

    # Define noiseless controls
    noiseless_drive = qctrl.types.filter_function.Drive(
        control=[
            qctrl.types.ComplexSegmentInput(duration=duration, value=value)
            for duration, value in zip(
                pulse.durations, pulse.rabi_rates * np.exp(1j * pulse.azimuthal_angles)
            )
        ],
        operator=sigma_m / 2,
    )

    noiseless_shift = qctrl.types.filter_function.Shift(
        control=[
            qctrl.types.RealSegmentInput(duration=duration, value=value)
            for duration, value in zip(pulse.durations, pulse.detunings)
        ],
        operator=sigma_z / 2,
    )

    # Define complex control with multiplicative noise
    noisy_drive = qctrl.types.filter_function.Drive(
        control=noiseless_drive.control, operator=noiseless_drive.operator, noise=True
    )

    # Define additive noise
    dephasing_noise = qctrl.types.filter_function.Drift(
        noise=True, operator=sigma_z / 2
    )

    # The filter function of the amplitude noise excludes the term
    # for the dephasing noise (the drift)
    amplitude_result = qctrl.functions.calculate_filter_function(
        duration=pulse.duration,
        frequencies=frequencies,
        sample_count=3000,
        drives=[noisy_drive],
        shifts=[noiseless_shift],
        result_scope="NO_FREQUENCY_DOMAIN_NOISE_OPERATORS",
    )
    assert not amplitude_result.errors

    # The filter function of the dephasing noise excludes the term
    # for the amplitude noise (the noisy drive)
    dephasing_result = qctrl.functions.calculate_filter_function(
        duration=pulse.duration,
        frequencies=frequencies,
        sample_count=None,
        drives=[noiseless_drive],
        shifts=[noiseless_shift],
        drifts=[dephasing_noise],
        result_scope="ALL",
    )
    assert not dephasing_result.errors

    return ControlScheme(
        name=pulse.name,
        amplitude_filter_function_samples=amplitude_result.samples,
        dephasing_filter_function_samples=dephasing_result.samples,
    )

Calculate and visualize the filter functions

Here we calculate filter functions and set the frequencies in the array to be spaced logarithmically, because we will also be interested in plotting the filter functions on a log-log graph.

# Define filter function parameters
interpolated_frequencies = omega_max * np.logspace(-8, 0, 1000, base=10)

# We expect two calls to calculate_filter_function for each pulse as we are evaluating
# for both drifts and shifts
control_schemes = [
    create_control_scheme(pulse, interpolated_frequencies) for pulse in pulses
]
Your task calculate_filter_function has started.
Your task calculate_filter_function has completed.
Your task calculate_filter_function has completed.
Your task calculate_filter_function has started.
Your task calculate_filter_function has completed.
Your task calculate_filter_function has completed.
Your task calculate_filter_function has started.
Your task calculate_filter_function has completed.
Your task calculate_filter_function has completed.
Your task calculate_filter_function has started.
Your task calculate_filter_function has completed.
Your task calculate_filter_function has completed.

Amplitude noise

plot_filter_functions(
    plt.figure(),
    {
        scheme.name: map(asdict, scheme.amplitude_filter_function_samples)
        for scheme in control_schemes
    },
)

We see that the BB1 and CORPSE-in-BB1 pulses show small values of the filter function at low frequencies, corresponding to their known amplitude-noise-suppressing characteristics. CORPSE and primitive gates do not exhibit any suppression of amplitude noise.

Dephasing noise

plot_filter_functions(
    plt.figure(),
    {
        scheme.name: map(asdict, scheme.dephasing_filter_function_samples)
        for scheme in control_schemes
    },
)

We see that the CORPSE and CORPSE-in-BB1 pulses show small values of the filter function at low frequencies, corresponding to their known dephasing-noise-suppressing characteristics. BB1 and primitive gates do not exhibit any suppression of dephasing noise.

Calculate the infidelity with respect to noise power spectral densities using noise operators

The frequency domain noise operators are computed as part of the filter function calculation. They are available in the result samples when the result_scope parameter is set appropriately. In this case we have configured the filter function calculation for the dephasing noise to include the operators.

With these in hand we can calculate the expected infidelity of a control in the presence of noise channels exhibiting specific quantitative power spectra. We define dimensionful power spectra below and compare the performance of the pulses.

scheme = [s for s in control_schemes if s.name == "primitive"][0]
last_dephasing_sample = scheme.dephasing_filter_function_samples[-1]
print(f"Scheme: {scheme.name}")
print(f"Frequency domain noise operator at f={last_dephasing_sample.frequency}Hz")
print(last_dephasing_sample.frequency_domain_noise_operator)
Scheme: primitive
Frequency domain noise operator at f=6283185.307179586Hz
[[-1.00946617e-08-2.11765689e-08j  1.60661531e-09+3.37035562e-09j]
 [-1.60661531e-09-3.37035562e-09j  1.00946617e-08+2.11765689e-08j]]
# Define amplitude noise power spectral density
psd_amplitude = lambda frequency: 1e-2 / (frequency ** 2 + 1)

# Define dephasing noise power spectral density
psd_dephasing = lambda frequency: 1e12 / (frequency ** 2 + 1)

for scheme in control_schemes:
    # Calculate overlap integral for amplitude noise
    frequencies = np.array(
        [sample.frequency for sample in scheme.amplitude_filter_function_samples]
    )
    inverse_powers = np.array(
        [sample.inverse_power for sample in scheme.amplitude_filter_function_samples]
    )
    amplitude_infidelity = np.trapz(
        y=inverse_powers * psd_amplitude(frequencies),
        x=frequencies,
    )

    # Calculate overlap integral for dephasing noise
    frequencies = np.array(
        [sample.frequency for sample in scheme.dephasing_filter_function_samples]
    )
    inverse_powers = np.array(
        [sample.inverse_power for sample in scheme.dephasing_filter_function_samples]
    )
    dephasing_infidelity = np.trapz(
        y=inverse_powers * psd_dephasing(frequencies),
        x=frequencies,
    )

    # Add contributions from the two noises
    infidelity = amplitude_infidelity + dephasing_infidelity

    # Print results
    print()
    print(f"Scheme: {scheme.name}")
    print(f"Infidelity due to amplitude noise: {amplitude_infidelity:.2e}")
    print(f"Infidelity due to dephasing noise: {dephasing_infidelity:.2e}")
    print(f"Total infidelity: {infidelity:.2e}")
Scheme: primitive
Infidelity due to amplitude noise: 3.72e-02
Infidelity due to dephasing noise: 3.82e-02
Total infidelity: 7.54e-02

Scheme: BB1
Infidelity due to amplitude noise: 7.50e-07
Infidelity due to dephasing noise: 3.82e-02
Total infidelity: 3.82e-02

Scheme: CORPSE
Infidelity due to amplitude noise: 3.73e-02
Infidelity due to dephasing noise: 7.50e-07
Total infidelity: 3.73e-02

Scheme: CORPSE in BB1
Infidelity due to amplitude noise: 2.84e-06
Infidelity due to dephasing noise: 2.19e-06
Total infidelity: 5.02e-06

Summary

The plots and infidelities show that the BB1 controls perform better against (low-frequency) amplitude noise than the primitive $\pi$-pulse, but perform just like the primitive in the case of dephasing noise. This is expected, as BB1 is one of the control-error-compensating driven controls. The inverse is true of the pure CORPSE controls: they perform as poorly as the primitive against amplitude noise, but perform better against (low-frequency) dephasing noise. This behavior is expected as CORPSE is a dephasing-error-compensating driven control. Finally, the dephasing-and-control-error-compensating driven control CORPSE in BB1 performs better than the primitive for low-frequencies of both noises (albeit less so than the controls specialized for one specific kind of noise).

We have thus demonstrated how the Q-CTRL Python package can be used to characterize the sensitivity of different controls to time-dependent noise channels by calculating their corresponding filter functions.