How to perform noise spectroscopy on arbitrary noise channels

Boulder Opal enables you to characterize and then reconstruct power spectral densities for arbitrary noise processes affecting your system. These power spectral densities can provide useful information for designing robust controls that improve the performance of the system in the presence of noise.

In this notebook we show how to characterize and reconstruct noise spectra using Boulder Opal.

Noise spectroscopy fundamentally relies on performing probe measurements of your quantum system and then inferring information about the system noise from the returned measurement outcomes. For arbitrary controls this challenge is efficiently enabled by calculating their spectral response using the filter function framework (applicable to D-dimensional Hilbert spaces). The measurement associated with the $i^{th}$ applied control can be expressed in discretized form as follows: \begin{align} P_i = \sum_{j=1}^N F^i_j S_j \Delta \nu, \end{align} where the filter function of the control $F^i$ and the spectrum $\mathbf{S}$ are sampled with spectral resolution $\Delta \nu = (\nu_{\rm max}-\nu_{\rm min})/N$.

Spectral reconstruction requires finding the matrix $\mathbf{M}$, such that $\mathbf{P} = \mathbf{M}\mathbf{S}$ and then inverting to find $\mathbf{S}$. Therefore, appropriate selection of controls enables you to interrogate the spectrum window ($\nu_{\rm min}$,$\nu_{\rm max}$) in great detail by overlapping the filter functions to sample the spectrum as densely as needed.

Summary workflow

By following these steps, you'll build a spectrum reconstruction technique customized to your needs:

1. Determine probe controls that are sensitive to noise at different frequencies,
2. Collect experimental measurements $\mathbf{P}$,
3. Construct the mapping $\mathbf{M}$ between measurements $\mathbf{P}$ and the power spectrum $\mathbf{S}$, where $\mathbf{P=MS}$ using the filter-function formalism,
4. Invert the mapping $M$ in order to reconstruct the interrogated section of the spectrum $\mathbf{S} = \mathbf{M^{-1}}\mathbf{P}$.

1. Determine probe controls that are sensitive to noise at different frequencies

In principle any set of measurements can be used to probe a noise spectrum. We provide a framework in which to efficiently define modulated-Gaussian probe controls which give ready access to frequency-resolved information about your system. The assumption here is that you are able to employ shaped pulses in your system.

The specific set of controls you should employ will be determined by the specific reconstruction at hand.

Considerations for efficient spectral sampling.

Let the probe-pulse duration be $T$, the minimum segment duration $t$; generally speaking, the maximum noise frequency to which a characterization pulse is sensitive is $0.5/t$ (the Nyquist frequency). Therefore, to reconstruct a spectrum defined up to some frequency $f_\text{max}$, the minimum segment duration should be chosen as $t\approx 0.5/f_\text{max}$. The bandwidth of each filter function, and thus the resolution of the resulting reconstruction, is roughly $2/T$. Therefore, for a pulse duration $T$, you should not expect to detect features of the noise spectrum any narrower than $2/T$.

It is important that the set of filter functions provides good coverage of the whole spectrum, with no gaps, via appropriate shifting of the band-center of each control's spectral response. This may be accomplished by chooiseing the number of probe controls $m \approx 0.25T/t$. A significantly smaller value will lead to gaps in the coverage of the spectrum, while a significantly larger value does not yield any improvement in coverage and is thus unnecessary.

Note that if the noise spectrum is to be reconstructed in a logarithmic manner (where the desired precision decreases as the frequency increases), it may be beneficial to perform the reconstruction in separate frequency intervals, each using an appropriately tuned control-spectral-response.

Probing multiplicative amplitude noise on a control

To perform spectroscopy on the multiplicative noise present on a control channel (for example, a single or multi-qubit drive), start by generating characterization pulses using new_modulated_gaussian_control from Q-CTRL Open Controls, which requires:

• maximum Rabi rate,
• minimum segment duration (which describes the resolution for pulse shaping),
• duration (which describes the total duration of the pulse),
• modulation frequency.

Note that the modulated Gaussian control pulses are sensitive to amplitude noise at specific frequencies (determined by the modulation frequency), and therefore are suitable for reconstructing amplitude noise. The maximum Rabi rate $\Omega_\text{max}$ generally controls the sensitivity of the pulses; a larger maximum Rabi rate leads to more prominent spectral response, and therefore a higher-quality reconstruction. Information on which Gaussian pulses to select is provided in the description of step 2 below.

Probing additive ambient noise processes

In order to probe ambient dephasing it is generally useful to employ tunable dynamic decoupling sequence such as Carr-Purcell-Meiboom-Gill (CPMG) or XY families from Q-CTRL Open Controls. Here, the peak frequency is set by the base period of the sequence while, as described below, the overall sequence duration determines the resolution of the probe. Typically, one will fix the sequence duration, $T$, and vary the number of control pulses applied within that window.

Note: The existence of high-frequency harmonics in the spectral response of dynamic decoupling sequences can introduce artefacts in spectral reconstruction. A conventional approach using a method developed by G. A. Alvarez, D. Suter places strict constraints on the available probe sequences in order to mitigate these challenges. To see a performance comparison between this algorithm and the methods described here see the Performing narrow-band magnetic-field spectroscopy with NV centers Application Note.

2. Collect experimental measurements

In a noise-free system, each of the generated control pulses would effect an identity gate. In the presence of noise, however, there will be some non-zero infidelity which carries information about the noise in the system. Format controls for your hardware and return measures of infidelity for target operations.

3. Calculate filter functions for controls

With control protocols selected, the next step is to calculate the filter function for each control sequence using graph.filter_function (see the user guide for filter function calculation).

Once calculated the filter functions provide a simple representation of the spectral sensitivity of a set of measurements. It's important to use this opportunity to evaluate the spectral coverage of your estimation.

The combined filter functions and measurements must be packaged for the reconstruction engine using measurement records.

4. Reconstruct spectrum from measurements of your quantum system

Now we are ready to reconstruct the noise PSD by calling the function qctrl.functions.calculate_noise_reconstruction. Boulder Opal offers two different efficient inversion procedures powered by TensorFlow and purpose built for spectral reconstruction. They are based on single value decomposition (SVD) and convex optimization (CVX). Both methods provide parameter free estimation needed for reconstruction of unknown power spectrum densities.

Use the SVD method when noise power spectra are expected to be smoothly varying and numerical efficiency is a key driver. Unphysical oscillations in the reconstructed spectrum may appear when these conditions are not met (see the discussion in this paper). qctrl.functions.calculate_noise_reconstruction uses the SVD method by default.

Use the CVX procedure in the presence of rapidly varying noise terms (for example, discrete spurs). To use the CVX method, you need to set up the qctrl.types.noise_reconstruction.ConvexOptimization object first, which includes the lower and upper bound of the noise PSD to be reconstructed, and a hyperparameter. The hyperparameter is used to control the balance between fitting the data and achieving a smooth noise PSD—a low value leads to a solution that fits the data very well at the expense of smoothness, while a high value leads to a solution is very smooth but does not fit the data as well. We recommend use of the conventional cross-validation method to tune its value.

Results are returned by invoking reconstruction result object.

For further mathematical details and a comparison between the reconstruction methods, refer to Ball et al..

Example: Noise spectroscopy on drive-amplitude noise for a single qubit

In this example we consider a single-qubit system driven by a controllable Rabi rate that is subject to amplitude noise. The Hamiltonian of the system is:

\begin{align} H(t) = &\frac{1+\beta_\Omega(t)}{2} \Omega(t) \sigma_x, \end{align}

where $\Omega(t)$ is the controllable Rabi rate, $\beta_\Omega(t)$ is a fractional time-dependent amplitude fluctuation process and $\sigma_x$ is the Pauli X matrix. We assume that the noise process $\beta_\Omega(t)$ consists of pink noise with a small Gaussian feature, as shown below.

import matplotlib.pyplot as plt
import numpy as np
from qctrlopencontrols import new_modulated_gaussian_control
from qctrlvisualizer import get_qctrl_style, plot_controls, QCTRL_STYLE_COLORS

plt.style.use(get_qctrl_style())

from qctrl import Qctrl

# Start a Boulder Opal session.
qctrl = Qctrl(verbosity="QUIET")

def gaussian(frequencies, offset, width):
    return np.exp(-0.5 * (frequencies - offset) ** 2 / width**2) / (
        np.sqrt(2 * np.pi) * width
    )


def pink(frequencies, frequency_cutoff, power):
    return frequency_cutoff ** (power - 1) / (
        frequencies**power + frequency_cutoff**power
    )


frequencies = np.linspace(0.0, 0.5e6, 1001)

amplitude_noise = 0.5e-11 * (
    25 * pink(frequencies=frequencies, frequency_cutoff=0.1 * 1e6, power=6)
    + gaussian(frequencies=frequencies, offset=0.2 * 1e6, width=10 * 1e3)
)

plt.plot(frequencies / 1e6, amplitude_noise * 1e6)
plt.fill_between(frequencies / 1e6, 0, amplitude_noise * 1e6, alpha=0.25)
plt.xlabel("Frequency (MHz)")
plt.ylabel("Power density (1/MHz)")
plt.title("Amplitude noise spectrum")
plt.show()

Determine (and visualize) probe controls that are sensitive to noise at different frequencies

# Define system parameters.
omega_max = 2 * np.pi * 250 * 1e6  # Hz
duration = 100 * 1e-6  # s
minimum_segment_duration = 1 * 1e-6  # s
pulse_count = 30

# Generate controls using in-built function for Gaussian probe controls.
noise_characterization_result = [
    new_modulated_gaussian_control(omega_max, minimum_segment_duration, duration, f)
    for f in np.linspace(0, 0.5 / minimum_segment_duration, pulse_count, endpoint=False)
]

# Visualize first five pulses.
plot_controls(
    {
        f"Pulse #{idx}": qctrl.utils.pwc_arrays_to_pairs(
            pulse.durations, pulse.rabi_rates * np.cos(pulse.azimuthal_angles)
        )
        for idx, pulse in enumerate(noise_characterization_result[:5])
    }
)

Collect experimental measurements

In this notebook we include code to simulate noisy measurements under the applied probe controls. In a laboratory setting this would be replaced with experimental measurements

# Define noise parameters.
trajectory_count = 150
noise_segment_count = 100
frequency_step = np.diff(frequencies)[0]

graph = qctrl.create_graph()

# Ensure all pulses use the same segmentation.
durations = noise_characterization_result[0].durations
for pulse in noise_characterization_result:
    assert np.array_equal(durations, pulse.durations)

# Create (pulse_count, T) array with the values of the controls.
values = np.array(
    [
        pulse.rabi_rates * np.cos(pulse.azimuthal_angles)
        for pulse in noise_characterization_result
    ]
)

# Create (pulse_count,)-batch of Pwcs representing the controls.
control = graph.pwc(durations=pulse.durations, values=values, time_dimension=-1)

# Create (trajectory_count, 1)-batch of Pwcs with the noises.
noise = graph.discretize_stf(
    stf=graph.random_colored_noise_stf_signal(
        power_spectral_density=amplitude_noise,
        frequency_step=frequency_step,
        batch_shape=(trajectory_count, 1),
    ),
    duration=np.sum(pulse.durations),
    segment_count=noise_segment_count,
)

# Create (trajectory_count, pulse_count)-batch of Pwcs with the Hamiltonians.
hamiltonian = 0.5 * (1 + noise) * control * graph.pauli_matrix("X")

# Calculate the (trajectory_count, pulse_count)-batch of infidelities.
graph.infidelity_pwc(
    hamiltonian=hamiltonian,
    target=graph.target(operator=graph.pauli_matrix("I")),
    name="infidelities",
)

# Execute graph.
infidelities = qctrl.functions.calculate_graph(
    graph=graph, output_node_names=["infidelities"]
).output["infidelities"]["value"]

# Calculate average/std of infidelities over the trajectories.
simulated_infidelities = np.mean(infidelities, axis=0)
simulated_infidelity_uncertainties = np.std(infidelities, axis=0, ddof=1) / np.sqrt(
    trajectory_count
)

Calculate (and visualize) filter functions for controls

# Filter functions.
sample_count = 200

filter_functions = []

for pulse in noise_characterization_result:
    graph = qctrl.create_graph()

    shift = graph.pwc(
        durations=pulse.durations,
        values=pulse.rabi_rates * np.exp(1j * pulse.azimuthal_angles),
    )

    control_hamiltonian = graph.hermitian_part(shift * graph.pauli_matrix("X") / 2)

    graph.filter_function(
        control_hamiltonian=control_hamiltonian,
        noise_operator=control_hamiltonian,
        frequencies=frequencies,
        sample_count=sample_count,
        name="filter_function",
    )

    result = qctrl.functions.calculate_graph(
        graph=graph, output_node_names=["filter_function"], execution_mode="EAGER"
    )

    filter_functions.append(result.output["filter_function"])

# Plot all filter functions.
for filter_function in filter_functions:
    plt.plot(filter_function["frequencies"] * 1e-6, filter_function["inverse_powers"])
plt.xlabel("Frequency (MHz)")
plt.ylabel("Inverse power")
plt.show()

measurement_records = []
for idx, filter_function in enumerate(filter_functions):
    filter_function_samples = [
        qctrl.types.noise_reconstruction.FilterFunctionSample(
            frequency=frequency, inverse_power=inverse_power
        )
        for frequency, inverse_power in zip(
            filter_function["frequencies"], filter_function["inverse_powers"]
        )
    ]
    measurement_record = qctrl.types.noise_reconstruction.Measurement(
        infidelity=simulated_infidelities[idx],
        infidelity_uncertainty=simulated_infidelity_uncertainties[idx],
        filter_functions=[
            qctrl.types.noise_reconstruction.FilterFunction(
                samples=filter_function_samples
            )
        ],
    )
    measurement_records.append(measurement_record)

Reconstruct (and visualize) noise spectrum

Now we are ready to reconstruct the noise PSD by calling the function qctrl.functions.calculate_noise_reconstruction. This function provides two different methods to perform reconstruction, namely the singular value decomposition (SVD) method and convex optimization (CVX) method. You can check the note part of the function documentation for technical details about these two methods. Here we use a hard-coded hyperparameter for the CVX method.

# Perform reconstruction with singular value decomposition.
svd_result = qctrl.functions.calculate_noise_reconstruction(
    measurements=measurement_records
)

# Configure convex optimization method.
cvx_configuration = qctrl.types.noise_reconstruction.ConvexOptimization(
    power_density_lower_bound=0,
    power_density_upper_bound=1,
    regularization_hyperparameter=1e-2,
)

# Perform reconstruction with convex optimization.
cvx_result = qctrl.functions.calculate_noise_reconstruction(
    measurements=measurement_records,
    method=qctrl.types.noise_reconstruction.Method(
        convex_optimization=cvx_configuration
    ),
)

# Extract data from reconstruction results.
svd_power_spectral_densities = svd_result.power_spectral_densities[0].samples
cvx_power_spectral_densities = cvx_result.power_spectral_densities[0].samples

sample_frequencies = np.array(
    [sample.frequency for sample in svd_power_spectral_densities]
)

sample_psd_svd = np.array(
    [sample.power_density for sample in svd_power_spectral_densities]
)
sample_psd_uncertainties_svd = np.array(
    [sample.power_density_uncertainty for sample in svd_power_spectral_densities]
)

sample_psd_cvx = np.array(
    [sample.power_density for sample in cvx_power_spectral_densities]
)
sample_psd_uncertainties_cvx = np.array(
    [sample.power_density_uncertainty for sample in cvx_power_spectral_densities]
)

# Plot the noise spectral densities with their uncertainties.
plt.plot(frequencies * 1e-6, amplitude_noise * 1e6, label="Original")
plt.plot(sample_frequencies * 1e-6, sample_psd_svd * 1e6, label="Reconstructed SVD")
plt.fill_between(
    frequencies * 1e-6,
    (sample_psd_svd - sample_psd_uncertainties_svd) * 1e6,
    (sample_psd_svd + sample_psd_uncertainties_svd) * 1e6,
    alpha=0.4,
    hatch="/" * 5,
    facecolor="none",
    edgecolor=QCTRL_STYLE_COLORS[1],
    linewidth=0.2,
)
plt.plot(sample_frequencies * 1e-6, sample_psd_cvx * 1e6, label="Reconstructed CVX")
plt.fill_between(
    frequencies * 1e-6,
    (sample_psd_cvx - sample_psd_uncertainties_cvx) * 1e6,
    (sample_psd_cvx + sample_psd_uncertainties_cvx) * 1e6,
    alpha=0.4,
    hatch="\\" * 5,
    facecolor="none",
    edgecolor=QCTRL_STYLE_COLORS[2],
    linewidth=0.2,
)
plt.legend()
plt.xlabel("Frequency (MHz)")
plt.ylabel("Power density (1/MHz)")
plt.show()

Summary

We see that the reconstructed noise spectral density matches the original noise spectral density reasonably closely, especially in terms of the large-scale trends. By taking more precise infidelity data (with lower uncertainties) we would see an even closer match between the two spectra. We have thus demonstrated how to use Boulder Opal to reconstruct the noise spectral density of an amplitude noise process affecting a single qubit system. Similar procedures may be employed to reconstruct noise spectral densities for other noise processes (for example stochastic dephasing noises) in other types of systems.

This notebook was run using the following package versions. It should also be compatible with newer versions of the Q-CTRL Python package.