# Performing noise spectroscopy in superconducting hardware

Reconstructing noise power spectrum density in transmon qubits using dynamical decoupling sequences

Boulder Opal enables you to reconstruct power spectral densities (PSD) of different noise sources affecting your system. The characterization of these noise PSD provides valuable information that can be used, for example, in the design of controls and hardware diagnostics.

In this application note, we focus on the reconstruction of noise PSD for superconducting qubits, a system widely used in quantum computing platforms. The notebook starts by defining the model of a superconducting transmon system subject to a dephasing noise whose PSD we want to reconstruct. The reconstruction process employs a carefully designed set of dynamical decoupling (DD) pulse sequences which, in this example, will be CPMG sequences composed of customized $\pi$ and $\pi/2$ pulses. Once defined, the sequences are used in experiments to probe the system and obtain the measurements necessary to perform the PSD reconstruction. In the notebook, the experimental measurements will be replaced by outputs of the appropriate simulation of the physical system. Finally, we will input these measurements directly into Boulder Opal’s calculate_noise_reconstruction function to obtain the reconstructed noise spectrum.

## Imports and initialization

# Boulder Opal
from qctrl import Qctrl

qctrl = Qctrl()

from qctrlvisualizer import (
plot_controls,
get_qctrl_style,
display_bloch_sphere_from_density_matrices,
)
from qctrlopencontrols import new_drag_control

import numpy as np

import matplotlib.pyplot as plt

plt.style.use(get_qctrl_style())


## Model of the system: transmon subject to dephasing noise

We consider a single transmon subject to a dephasing noise described by the following Hamiltonian:

$H = \omega_c a^\dagger a + \frac{\alpha}{2} (a^\dagger)^2 a^2 + \left(\Omega (t) b + H.c.\right) + \eta(t) a^\dagger a,$

where $\omega_c$ is the transmon frequency, $\alpha$ the anharmonicity, $\Omega (t)$ is the complex Rabi rate and $\eta(t)$ is a time-dependent dephasing noise process consistent with the experimental power spectral density of the noise affecting the system. We can further simplify the Hamiltonian by going into the frame rotating with $\omega_c$: $$H = \frac{\alpha}{2} (a^\dagger)^2 a^2 + \left(\Omega (t) b + H.c.\right) + \eta(t) a^\dagger a.$$ The cell below sets up a function (hamiltonian_setup) that creates the Hamiltonian, assuming a three-level truncation of transmon system. The hamiltonian_setup also includes the dephasing noise operator corresponding to the PSD to be reconstructed as well as the pulsed controls that will take the form of the CPMG sequence subsequently. See this user guide for more information on setting up the Hamiltonian of a qubit under noisy conditions.

# system parameters and operators
dim = 3  # transmon dimension
alpha = -2 * np.pi * 400e6  # rad.Hz, transmon anharmonicity
a = np.diag(np.sqrt(np.arange(1, dim)), k=1)

def hamiltonian_setup(
controls,  # pulse sequence controls
noise_power_densities=None,  # noise PSD
noise_frequency_step=None,  # sampling rate of noise PSD
noise_ensemble_count=1000,  # number of noise instances
graph=None,
):

"""
This function sets up the transmon Hamiltonian for the given
set of pulse controls and the noise PSD.
"""

ham = 0
hamiltonian_duration = 0

# control terms
if "Omega" in controls:
drive_signal = graph.pwc(
values=controls["Omega"]["values"],
durations=controls["Omega"]["durations"],
time_dimension=-1,
)
ham += graph.hermitian_part(drive_signal * a)

duration_omega = np.sum(controls["Omega"]["durations"])
if duration_omega > hamiltonian_duration:
hamiltonian_duration = duration_omega

# noise term
if (noise_power_densities is not None) and (noise_frequency_step is not None):

segment_count = (int)(np.round(hamiltonian_duration * np.max(frequencies)) * 10)

eta_stf = graph.random_colored_noise_stf_signal(
power_spectral_density=noise_power_densities,
frequency_step=noise_frequency_step,
batch_shape=(noise_ensemble_count,),
)
eta_signal = graph.discretize_stf(
stf=eta_stf,
duration=hamiltonian_duration,
segment_count=segment_count,
name="noise_pwc",
)

ham += eta_signal * adag @ a

#  alpha term
alpha_signal = graph.constant_pwc(constant=-alpha, duration=hamiltonian_duration)
ham += alpha_signal * (adag @ adag @ a @ a) / 2

return ham


### Noise spectrum

For the demonstration in this notebook, we will consider a PSD containing both a $1/f$ noise component and also a thermal broadband term. This PSD is defined in the cell below.

# Define the non-Markovian noise power spectrum
def noise_spectrum(
):
"""
This fuction sets up noise PSD.
"""
return (
scale_factor / (frequencies + cutoff_frequency) ** exponent
)

frequency_step = 1e3
frequencies = np.arange(0, 5e6, frequency_step)
power_densities = noise_spectrum(
)
fig = plt.figure(figsize=(10, 5))
fig.suptitle("Dephasing noise spectral density")
plt.plot(frequencies / 1e6, power_densities / 1e6, color="#EB64B6")
plt.fill_between(frequencies / 1e6, 0, power_densities / 1e6, color="#F790DB")
plt.xlabel("Frequency (MHz)")
plt.ylabel("Power density (1/MHz)")
plt.show()


The noise PSD defined above produces a qubit with a T2 coherence time of about $20\mu$s, which is commensurate with the coherence times often encountered in superconducting systems.

## Reconstructing the noise PSD

### Generate diagnostic DD sequences with custom pulses

To perform the noise reconstruction, we define a series of CPMG sequences with a fixed duration of 10$\mu$s and varying orders (number of $\pi$ pulses in each sequence). This is one of the most straightforward ways of probing the noise PSD. It enables us to control the amplitude and width of the fundamental filter function peaks and keep them relatively constant for CPMG pulse sequences of different orders.

The CPMG pulse sequences are constructed by concatenating a series of $\pi$ pulses appropriately spaced. An ideal sequence would be composed of perfect instantaneous $\pi$ pulses but, in reality, the pulses in the sequences are finite and can be shaped in different ways to satisfy hardware specifications or constraints. In this example, we chose to use DRAG pulses as they help minimize the leakage to the qubit’s excited state and ensure that the quality of the filter functions is maintained. The DRAG $\pi/2$ and $\pi$ pulses defined below are directly imported from the Q-CTRL Open Controls library and converted into the appropriate Boulder Opal format for later use.

def convert_pulse_to_controls_list(pulse):
"""
This function takes in a pulse of Open Controls format and converts it to the Boulder Opal control format.

"""
controls = {}
controls["Omega"] = {}
controls["Omega"]["values"] = pulse.rabi_rates * np.exp(1j * pulse.azimuthal_angles)
controls["Omega"]["durations"] = pulse.durations
return controls

omega_max = 2 * np.pi * 100e6  # rad.Hz, transmon Rabi frequency

# DRAG pulses
segment_count = 100

rabi_rotation = np.pi / 1
pi_pulse = new_drag_control(
rabi_rotation=rabi_rotation,
segment_count=segment_count,
duration=3 * (rabi_rotation / omega_max),
width=(rabi_rotation / omega_max) / 2,
beta=2e-10,
)
pi_controls = convert_pulse_to_controls_list(pi_pulse)

rabi_rotation = np.pi / 2
pi_half_pulse = new_drag_control(
rabi_rotation=rabi_rotation,
segment_count=segment_count,
duration=3 * (rabi_rotation / omega_max),
width=(rabi_rotation / omega_max) / 2,
beta=2e-10,
)
pi_half_controls = convert_pulse_to_controls_list(pi_half_pulse)

plot_controls(plt.figure(), pi_controls, polar=0)
plot_controls(plt.figure(), pi_half_controls, polar=0)
plt.show()


The properties of the pulses can have an impact on the filtering characteristics of the DD sequences. For this reason, when defining these pulses, the filter functions of the sequences should be carefully inspected for systematic anomalies. It’s useful to look out for occurrence of 0Hz peaks across different filter functions. Such anomalous peaks can introduce bias in the PSD estimates, especially in case of exponentially-shaped noise PSDs.

With the custom pulses defined, we can write a function that generates the CPMG sequences using them.

def custom_CPMG(duration, order, pi_controls, pi_half_controls):
"""
This function helps generate CPMG sequences using custom pi and pi/2 pulses.
"""
if order == 0:
seq_values = np.concatenate(
(
pi_half_controls["Omega"]["values"],
[0],
pi_half_controls["Omega"]["values"] * np.exp(1j * np.pi),
)
)
seq_durations = np.concatenate(
(
pi_half_controls["Omega"]["durations"],
[duration],
pi_half_controls["Omega"]["durations"],
)
)

else:

seq_values = np.tile(
np.concatenate(
([0.0], pi_controls["Omega"]["values"] * np.exp(1j * np.pi / 2), [0.0])
),
order,
)
tau = duration / order
seq_durations = np.tile(
np.concatenate(
(
[(tau - np.sum(pi_controls["Omega"]["durations"])) / 2],
pi_controls["Omega"]["durations"],
[(tau - np.sum(pi_controls["Omega"]["durations"])) / 2],
)
),
order,
)

seq_durations[0] = (
seq_durations[0] - np.sum(pi_half_controls["Omega"]["durations"]) / 2
)
seq_durations = np.insert(
seq_durations, 0, pi_half_controls["Omega"]["durations"]
)
seq_values = np.insert(seq_values, 0, pi_half_controls["Omega"]["values"])

seq_durations[-1] = (
seq_durations[-1] - np.sum(pi_half_controls["Omega"]["durations"]) / 2
)
seq_durations = np.append(seq_durations, pi_half_controls["Omega"]["durations"])
seq_values = np.append(
seq_values, pi_half_controls["Omega"]["values"] * np.exp(1j * np.pi)
)

return {"Omega": {"durations": seq_durations, "values": seq_values}}


Using the previous function, we define a list of CPMG sequences and plot one of them as an example.

CPMG_sequence_duration = 10e-6
CPMG_sequence_orders = np.arange(0, 40, 2)

CPMG_controls_list = [
custom_CPMG(CPMG_sequence_duration, order, pi_controls, pi_half_controls)
for order in CPMG_sequence_orders
]

# plot a particular sequence
plot_controls(plt.figure(), CPMG_controls_list[2], polar=True)
plt.show()


Note that various classes of DD sequences can be chosen for noise reconstruction depending on the characteristics of the system. We picked CPMG in this example because of its simplicity and high degree of efficiency in filtering dephasing noise. However, if there were Rabi rate fluctuations in the controls acting on the system, for example, we would pick a sequence with higher robustness to this Rabi rate noise, such as the XY8 sequence.

### Compute filter functions

Next, we calculate the filter function associated with each pulse sequence. See the How to calculate and use filter functions user guide for more details.

with qctrl.parallel():
filter_function_data = [
qctrl.functions.calculate_filter_function(
duration=np.sum(control["Omega"]["durations"]),
frequencies=frequencies,
sample_count=2000,
projection_operator=np.diag([1, 1, 0]),
drives=[
qctrl.types.filter_function.Drive(
control=[
qctrl.types.ComplexSegmentInput(duration=d, value=v)
for d, v in zip(
control["Omega"]["durations"], control["Omega"]["values"]
)
],
operator=a / 2,
)
],
drifts=[
qctrl.types.filter_function.Drift(
),
],
)
for control in CPMG_controls_list
]

Your task calculate_filter_function (action_id="1196215") has completed.

# extract filter functions from data
filter_functions = [
np.asarray([sample.inverse_power for sample in ffd.samples])
for ffd in filter_function_data
]

ff_peaks = np.array(
[frequencies[np.where(ff == np.max(ff))[0][0]] for ff in filter_functions]
)

# Plot
fig, ax = plt.subplots(figsize=(15, 5))
for ff in filter_functions:
ax.plot(frequencies / 1e6, ff, linewidth=0.8)

ax.set_xlabel("Frequency (MHz)")
ax.set_ylabel("Inverse power (Hz)")
fig.suptitle("Filter functions and power spectrum")
ax.set_xlim([0, 5])

ax2 = ax.twinx()
ax2.fill_between(frequencies / 1e6, power_densities / 1e6, 0, alpha=0.1, color="k")
ax2.set_ylabel("Spectrum power (1/MHz)")

plt.show()


Observe the fundamental peaks for the different CPMG filter functions (colored lines) interrogating the spectral window of interest. Note that the 0 order peak corresponds to the Ramsey sequence. The FWHM of filter function peaks acts a natural limit to spectral resolution. Longer sequence durations produce narrower peeks requiring use of higher order sequences to sample the same spectral window. Ultimately, the maximal frequency that can be sampled is limited by the duration of the $\pi$ pulses.

The key to noise reconstruction is making the fundamental filter function peaks overlap adequately to ensure they sample the noise PSD uniformly, as in the figure above. Note that in general, DD sequences of different durations or classes can be mixed together in our reconstruction method. However, in such cases the sampling coverage should be carefully managed so as to ensure no part of the interrogation window is left undersampled, that is, devoid of overlapping fundamental filter function peaks.

Also note that the CPMG, as well as other DD sequences, present higher order harmonics in addition to their main peak (see the smaller peaks in the plot). Those with significant presence outside the interrogation region will cause spectral leakage and lead to artifacts in the reconstructed spectrum. This can be especially problematic in reconstruction of power spectra containing complex features. In such cases, spectral leakage can be suppressed by using Slepian controls in place of DD sequences, as described in this application note.

### Apply sequences and perform measurements

In a real experiment, the next step would be to apply the set of diagnostic DD sequences to the qubit and measure its population. Here, however, we will simulate such experiments. The cell below set up the function calculate_single_pulse_under_noise that calculates the evolution of the system and outputs the probabilities of measuring the different transmon states at the end of a sequence averaged over an ensemble of different noise samples.

def calculate_single_pulse_under_noise(
controls,
noise_power_densities=power_densities,
noise_frequency_step=frequency_step,
noise_ensemble_count=1,
):
"""
This performs state evolution under the Hamiltonian.
It assumes that the 3 level system starts in the ground. The return is
the average state for all state evolutions under individual noise instances.
"""

graph = qctrl.create_graph()
controls_duration = np.sum(controls["Omega"]["durations"])

hamiltonian = hamiltonian_setup(
controls=controls,
noise_power_densities=noise_power_densities,
noise_frequency_step=noise_frequency_step,
noise_ensemble_count=noise_ensemble_count,
graph=graph,
)

# calculate the evolution unitaries
unitaries = graph.time_evolution_operators_pwc(
hamiltonian=hamiltonian, sample_times=np.array([controls_duration])
)

initial_state = np.array([1, 0, 0])

evolved_ensemble = unitaries @ initial_state[:, None]
probabilities_ensemble = graph.abs(evolved_ensemble[:, 0, :, 0]) ** 2

probabilities = graph.sum(probabilities_ensemble, axis=0) / noise_ensemble_count
probabilities.name = "probabilities"

graph_result = qctrl.functions.calculate_graph(
graph=graph, output_node_names=["probabilities"]
)
return graph_result


Using the functions defined above, we can simulate a series of experiments where the final populations are measured after the application of CPMG sequences of various orders.

noise_ensemble_count = 4000

with qctrl.parallel():
results_list = [
calculate_single_pulse_under_noise(
controls=controls,
noise_power_densities=power_densities,
noise_frequency_step=frequency_step,
noise_ensemble_count=noise_ensemble_count,
)
for controls in CPMG_controls_list
]

CPMG_probabilities = np.array(
[result.output["probabilities"]["value"] for result in results_list]
)

Your task calculate_graph (action_id="1196239") has completed.
Your task calculate_graph (action_id="1196241") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196242") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196243") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196244") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196245") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196246") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196247") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196248") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196249") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196250") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196251") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196252") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196253") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196254") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196255") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196256") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196257") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_graph (action_id="1196258") has started. You can use the qctrl.get_result method to retrieve previous results.

fig = plt.figure(figsize=(10, 5))
fig.suptitle("State Probability as a fuction of CPMG order")

plt.plot(CPMG_sequence_orders, CPMG_probabilities[:, 0], label=r"$|0\rangle$")
plt.plot(CPMG_sequence_orders, CPMG_probabilities[:, 1], label=r"$|1\rangle$")
plt.plot(CPMG_sequence_orders, CPMG_probabilities[:, 2], label=r"$|2\rangle$")
plt.plot(CPMG_sequence_orders, CPMG_probabilities * 0 + 0.5, ":")

plt.xlabel("CPMG order")
plt.ylabel("State Probability $P_{0,1,2}$")
plt.ylim([0, 1])
plt.legend()
plt.show()


### Process the measurements

We now construct a bounded infidelity metric $I_B=1-P_0$, where $P_0$ is the probability of the $|0\rangle$ state. In experiments, the response of this metric to noise is nonlinear and naturally defined in the interval [0.0,0.5]. Note that the filter function formalism assumes unbounded linear response to noise. To use filter functions with greater accuracy, we first linearize the measurement values, as done below (see Soare et. al for more details).

bounded_infidelity = 1 - CPMG_probabilities[:, 0]

unbounded_infidelity = -np.log(1 - bounded_infidelity * 2) / 2


Ideally, the duration of the sequence should be significantly shorter than the qubit’s T1 coherence time to minimize the impact of T1 decay over the evolution of the sequence. For durations comparable to the qubit’s T1 time, the effect should be accounted for by rescaling the signal based on the decay envelope.

### Reconstruct the noise PSD

We now reconstruct the noise in several simple steps. First, we select the frequency band that corresponds to the spectral range being interrogated by the filter function peaks - this is the spectral window over which we have the ability to reconstruct the PSD.

lower_sample_frequency = ff_peaks[0]
upper_sample_frequency = ff_peaks[-1]

sample_selection = (frequencies >= lower_sample_frequency) & (
frequencies <= upper_sample_frequency
)
sample_frequencies = frequencies[sample_selection]
sample_spectrum = power_densities[sample_selection]


Next, wrap the measurements along with their corresponding filter functions in the noise reconstructor format, as described in our user guide on How to perform noise spectroscopy on arbitrary noise channels.

# First build the measurement records with the filter function and infidelities generated above
shot_noise = 1e-4
measurement_records = []
sample_filter_functions = np.asarray(filter_functions)[:, sample_selection]

for idx, filter_function in enumerate(sample_filter_functions):
filter_function_samples = [
qctrl.types.noise_reconstruction.FilterFunctionSample(
frequency=frequency, inverse_power=inverse_power
)
for frequency, inverse_power in zip(sample_frequencies, filter_function)
]
measurement_record = qctrl.types.noise_reconstruction.Measurement(
infidelity=unbounded_infidelity[idx],
infidelity_uncertainty=shot_noise,
filter_functions=[
qctrl.types.noise_reconstruction.FilterFunction(
samples=filter_function_samples
)
],
)
measurement_records.append(measurement_record)


Here we call the qctrl.functions.calculate_noise_reconstruction method to generate the reconstructed noise PSD. It uses convex optimization to perform a non-parametric spectrum reconstruction without specific assumptions about the shape of the noise PSD.

# Configure CVX reconstruction method
cvx_configuration = qctrl.types.noise_reconstruction.ConvexOptimization(
power_density_lower_bound=0,
power_density_upper_bound=1e6,
regularization_hyperparameter=1e-11,
)

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

# Extract reconstructed PSD
psd_samples = noise_reconstruction_result_cvx.power_spectral_densities[0].samples
reconstructed_psd = np.array([sample.power_density for sample in psd_samples])
reconstructed_psd_uncertainty = np.array(
[sample.power_density_uncertainty for sample in psd_samples]
)

Your task calculate_noise_reconstruction (action_id="1196571") is currently in a queue waiting to be processed.
Your task calculate_noise_reconstruction (action_id="1196571") has started. You can use the qctrl.get_result method to retrieve previous results.

# Plot
fig, ax = plt.subplots(figsize=(15, 5))
ax.fill_between(
sample_frequencies / 1e6,
sample_spectrum / 1e6,
0,
color="lightgray",
label="True spectrum",
)
ax.plot(
sample_frequencies / 1e6,
reconstructed_psd / 1e6,
"-",
label="Reconstructed spectrum",
)

ax.set_xlabel("Frequency (MHz)")
ax.set_ylabel("Power density (MHz$^{-1}$)")
ax.legend()
fig.suptitle("Power spectrum reconstruction with Boulder Opal")
plt.show()


Here we can compare the reconstructed spectrum (solid line) with the true noise PSD used to simulate the results (shaded).

Note the discrepancies near each of the extremities of the spectrum. These artifacts originate from a reduced sampling level. The sides of the band are probed by a single filter function in contrast to the middle, where every frequency point is sampled by at least two overlapping filter function peaks. Those artifacts can be reduced by narrowing the spectral width of the filter function peaks. Furthermore, these reconstructed results can be used to propose an analytic model of noise PSD with relatively few degrees of freedom and perform a parametric reconstruction that would resolve such artifacts.

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

Package version
Python 3.9.12
matplotlib 3.5.1
numpy 1.21.5
scipy 1.7.3
qctrl 19.2.0
qctrlcommons 17.1.1
qctrlopencontrols 9.1.1
qctrltoolkit 1.6.1
qctrlvisualizer 3.3.0