How to evaluate control susceptibility to quasistatic noise

Characterize the robustness of a control pulse to quasi-static noise

The Q-CTRL Python package enables you to simply evaluate control performance by calculating susceptibility to quasi-static noise on multiple simultaneous noise channels. This process can provide a useful characterization of the robustness of candidate controls and would be familiar to anyone from a magnetic resonance background. In this notebook we show how to compute what we refer to as "quasi-static scans" using the Q-CTRL Python package.

It is common to calculate the susceptibility of a control to two different noise processes applied simultaneously. Some of the most common noise terms include

  • the time-dependent multiplicative noise on a drive (commonly denoted $\beta$),
  • the additive dephasing noise which acts independent of any controls (commonly denoted $\eta$).

We may scan across one of these noise processes in a 1D quasistatic scan or both in a two-dimensional grid of values $(\beta,\eta)$, and calculate the infidelity of the resulting gate in each setting.

Summary workflow

1. Establish graph-based noise-susceptibility calculation

Fundamentally the calculation being performed in a quasistatic scan is the calculation of infidelity for an operation in the presence of a static offset term in a Hamiltonian. This simulation is efficiently captured using graphs.

As described in the How to represent quantum systems using graphs user guide, you will build a graph with the qctrl.functions.calculate_graph function.

With the graph built, you can calculate the quasi-static scans by calling the qctrl.functions.calculate_graph function (with the graph you defined and the output_node_names we want to obtain).

2. Batch data for infidelity calculation

1D parameter scan

In a 1D noise-susceptibility scan you create a batch of noise signals over which the simulation is performed using the constant_pwc graph operation and an array of noise coefficient values taken as an input.

noise_signal = graph.constant_pwc(
    constant=noise_coefficients,
    duration=pulse.duration,
    batch_dimension_count=1,
)
noise_term = noise_signal * operator

This signal is added to the Hamiltonian for inclusion in an infidelity calculation via the infidelity_pwc graph operation. You can extract the infidelities from result.output["infidelities"]["value"] (where result is the object returned by qctrl.functions.calculate_graph) for any noise coefficients defined to capture the parameter over which a quasistatic scan is performed.

2D parameter scan

In a 2D parameter scan, calculations are typically parallelized by creating a two-dimensional batch of Hamiltonians, with each element of the batch corresponding to a pair of noise coefficient values. To achieve that, we use graph operations to define PWC signals for each term. We then multiply each of these by their corresponding operator to create the corresponding Hamiltonian term, and add them together to obtain the total Hamiltonian.

When the terms are added, as the batches for both quasi-static noises correspond to different dimensions, we obtain a batch of 2×2 PWC Hamiltonians (when adding PWCs, their value and batch dimensions are broadcasted separately). As a result, the infidelities returned by the infidelity_pwc node are a tensor, where each element gives the infidelity associated with a pair of noise values.

Worked example: 1D quasistatic susceptibility to a single dephasing-noise channel

In this example we consider a 1D quasistatic scan for a single dephasing noise channel affecting a qubit. The system is described by the Hamiltonian:

\begin{align*} H(t) = & \frac{\Omega(t)}{2} \sigma_- + \frac{\Omega^*(t)}{2} \sigma_+ + \frac{\eta(t)}{2} \sigma_z \end{align*}

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

In this case we consider two driven controls, primitive and CORPSE.

As described above, we will build a graph that we will execute with to perform the computation. For convenience, we wrap this calculation in a function that takes a pulse given by a Q-CTRL Open Controls function, defining the piecewise-constant (PWC) pulses for $\Omega(t)$. We extract the infidelities from result.output["infidelities"]["value"] (where result is the object returned by qctrl.functions.calculate_graph) and plot in a 1D graph.

import matplotlib.pyplot as plt
import numpy as np
from qctrlvisualizer import (
    get_qctrl_style,
    plot_controls,
)

plt.style.use(get_qctrl_style())

# Predefined pulse imports
import qctrlopencontrols

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_m = np.array([[0, 1], [0, 0]], dtype=complex)
sqrt_sigma_x = 0.5 * np.array([[1 + 1j, 1 - 1j], [1 - 1j, 1 + 1j]], dtype=complex)

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

# Define coefficient array for the noise values
dephasing_coefficients = np.linspace(-1.0, 1.0, 101) * omega_max

fig = plt.figure(figsize=(10, 5))
fig.suptitle("Dephasing noise susceptibility of different pulses")

# For each scheme, compute and plot the results of the quasi-static scan
for scheme_name, function in [
    ("CORPSE", qctrlopencontrols.new_corpse_control),
    ("primitive", qctrlopencontrols.new_primitive_control),
]:

    # Define pulse objects using pulses from Q-CTRL Open Controls
    pulse = function(
        rabi_rotation=total_rotation, azimuthal_angle=0.0, maximum_rabi_rate=omega_max
    )

    # Create graph
    graph = qctrl.create_graph()

    # Define Rabi coupling term
    rabi_signal = graph.pwc(
        durations=pulse.durations,
        values=pulse.rabi_rates * np.exp(1j * pulse.azimuthal_angles),
    )
    rabi_coupling_term = graph.pwc_operator_hermitian_part(rabi_signal * sigma_m)

    # Define dephasing term, a [101] batch of operators, created by multiplying
    # a [101] batch of (constant) PWC signals by the corresponding σz/2 operator
    dephasing_signal = graph.constant_pwc(
        constant=dephasing_coefficients,
        duration=pulse.duration,
        batch_dimension_count=1,
    )
    dephasing_term = dephasing_signal * sigma_z / 2

    # Build total Hamiltonian, a [101] batch of operators
    hamiltonian = rabi_coupling_term + dephasing_term

    # Calculate infidelities, a [101] tensor
    graph.infidelity_pwc(
        hamiltonian=hamiltonian, target=graph.target(sqrt_sigma_x), name="infidelities"
    )

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

    # Extract and plot infidelities
    infidelities = result.output["infidelities"]["value"]
    plt.plot(dephasing_coefficients / omega_max, infidelities, label=scheme_name)


plt.ylabel("Infidelity")
plt.xlabel(r"Relative dephasing coefficient $\eta/\Omega_\mathrm{max}$")
plt.legend()
plt.show()
Your task calculate_graph (action_id="625823") has completed.

Your task calculate_graph (action_id="625824") has completed.

Worked example: 2D pulse susceptibility to simultaneous amplitude and dephasing noise

In this example we will compare a series of composite $\pi$ pulses applied to a single qubit under amplitude and dephasing noise. The Hamiltonian of the quantum system is:

\begin{align*} H(t) = &\frac{1+\beta(t)}{2}\Big( \Omega(t) \sigma_- + \Omega^*(t) \sigma_+ \Big) + \frac{\eta(t)}{2} \sigma_z \end{align*}

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

Comparing the variation in infidelity across a 2D grid of $(\beta,\eta)$ gives information about the robustness of the appropriate control to simultaneous quasi-static noise on the two channels.

We consider the following driven control schemes for the controllable $\Omega(t)$ term, which are available from Q-CTRL Open Controls and described in the reference documentation: primitive, BB1, SK1, and CORPSE. Here we invoke these controls, define relevant parameters in the system simulation, and identify the arrays defining the grid of noise coefficients in the quasi-static scan.

# Define standard matrices
identity = np.array([[1, 0], [0, 1]], dtype=complex)
sigma_x = np.array([[0, 1], [1, 0]], dtype=complex)
sigma_z = np.array([[1, 0], [0, -1]], dtype=complex)
sigma_m = np.array([[0, 1], [0, 0]], dtype=complex)

# Define schemes for driven controls to compare
schemes = {
    name: {"function": function}
    for name, function in [
        ("primitive", qctrlopencontrols.new_primitive_control),
        ("BB1", qctrlopencontrols.new_bb1_control),
        ("SK1", qctrlopencontrols.new_sk1_control),
        ("CORPSE", qctrlopencontrols.new_corpse_control),
    ]
}

# Define control parameters
total_rotation = np.pi  # Target rotation angle here is a \pi pulse.
omega_max = 2 * np.pi * 1e6  # Hz, this is the maximum Rabi rate in the system.

# Define coefficient arrays for the noise values.
# We use a grid with 51 points for $\beta\in [-0.5,0.5]$ and 101 points for $\eta\in[-\Omega_{\mathrm{max}}, \Omega_{\mathrm{max}}]$.
amplitude_coefficients = np.linspace(-0.5, 0.5, 51)
dephasing_coefficients = np.linspace(-1.0, 1.0, 101) * omega_max

The infidelities returned by the infidelity_pwc node are a [51, 101] tensor (based on the sizes of $(\beta,\eta)$), where each element gives the infidelity associated with a pair of amplitude and dephasing coefficient values.

def calculate_quasi_static_scan(pulse):

    # Create graph
    graph = qctrl.create_graph()

    # [51, 1, T] array with values of (1+β)Ω(t) on each segment and for each value of β
    complex_rabi_rates_with_noise = (
        (1 + amplitude_coefficients[:, None, None])
        * pulse.rabi_rates
        * np.exp(1j * pulse.azimuthal_angles)
    )
    # [51, 1] batch of PWC signals for (1+β)Ω(t)
    # We specify time_dimension=2, as the PWC values array has two batching dimensions.
    rabi_signal = graph.pwc(
        values=complex_rabi_rates_with_noise,
        durations=pulse.durations,
        time_dimension=2,
    )
    # [51, 1] batch of 2×2 PWC operators for the coupling term, (1+β)[Ω(t)σ- + Ω*(t)σ+]
    rabi_coupling_term = graph.pwc_operator_hermitian_part(rabi_signal * sigma_m)

    # [1, 101] batch of 2×2 constant (with a single segment) PWC signals for η
    dephasing_signal = graph.pwc(
        values=dephasing_coefficients[None, :, None],
        durations=np.array([pulse.duration]),
        time_dimension=2,
    )
    # [1, 101] batch of 2×2 PWC operators for the dephasing drift term, η σz/2
    dephasing_term = dephasing_signal * sigma_z / 2

    # [51, 101] batch of 2×2 PWC operators for the total Hamiltonian,
    # with each element in the batch represents a unique pair of (β, η) values
    # (the [51, 1] and [1, 101] batches get broadcasted to [51, 101]).
    hamiltonian = rabi_coupling_term + dephasing_term

    # [51, 101] tensor with the infidelities
    graph.infidelity_pwc(
        hamiltonian=hamiltonian, target=graph.target(sigma_x), name="infidelities"
    )

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

    # Extract and return infidelities
    return result.output["infidelities"]["value"]

We can now call the function for each Q-CTRL Open Controls scheme and extract the calculated infidelities. We also store the values for the noise coefficients in the quasi-static scan to plot them later.

for scheme_name, scheme_objects in schemes.items():

    # Define pulse objects using pulses from Q-CTRL Open Controls
    pulse = scheme_objects["function"](
        rabi_rotation=total_rotation, azimuthal_angle=0.0, maximum_rabi_rate=omega_max
    )

    print(f"\nObtaining infidelities for {scheme_name} scheme...")

    # Create and execute the graph that calculates the infidelities
    infidelities = calculate_quasi_static_scan(pulse)

    # Save calculated infidelities
    scheme_objects["infidelities"] = infidelities

    # Save relevant quantities for later use
    scheme_objects["dephasing noise values"] = dephasing_coefficients
    scheme_objects["drive noise values"] = amplitude_coefficients
Obtaining infidelities for primitive scheme...
Your task calculate_graph (action_id="625825") has completed.

Obtaining infidelities for BB1 scheme...
Your task calculate_graph (action_id="625826") has completed.

Obtaining infidelities for SK1 scheme...
Your task calculate_graph (action_id="625827") has completed.

Obtaining infidelities for CORPSE scheme...
Your task calculate_graph (action_id="625828") has completed.

For 2D scans, density plots can provide intuitive visualizations of the noise robustness. With the grid of infidelities extracted above, it is simple to create such plots using the Matplotlib library.

for scheme_name, scheme_objects in schemes.items():

    fig, ax = plt.subplots(figsize=(10, 4))

    contours = plt.contour(
        scheme_objects["dephasing noise values"] / omega_max,
        scheme_objects["drive noise values"],
        scheme_objects["infidelities"],
        levels=[0.001, 0.01, 0.1, 0.5],
        colors="#680CEA",
    )
    plt.clabel(contours, inline=True)

    plt.imshow(
        scheme_objects["infidelities"],
        extent=[
            np.min(dephasing_coefficients) / omega_max,
            np.max(dephasing_coefficients) / omega_max,
            np.min(amplitude_coefficients),
            np.max(amplitude_coefficients),
        ],
        origin="lower",
        cmap=plt.cm.get_cmap("gray").reversed(),
        vmin=0,
        vmax=1,
    )

    cbar = plt.colorbar(pad=0.08)
    cbar.set_label("Pulse infidelity", labelpad=-50)
    plt.title(f"Noise susceptibility of {scheme_name} pulse")
    plt.ylabel(r"Amplitude coefficient $\beta$")
    plt.xlabel(r"Relative dephasing coefficient $\eta/\Omega_\mathrm{max}$")
    plt.show()

Summary

The plots demonstrate clearly the different noise susceptibilities of the different controls. For example, the BB1 control exhibits excellent robustness to amplitude noise (it retains a low infidelity across a wide range of amplitude noise values), but low dephasing robustness (infidelity increases rapidly as the dephasing coefficient varies from zero). The CORPSE control is the opposite: it is robust to dephasing noise, but highly susceptible to amplitude noise.

We have thus demonstrated how the Q-CTRL Python package can be used to characterize the robustness of different controls to quasi-static noise on multiple simultaneous noise channels.

Worked example: Combining quasi-static noise susceptibility with open quantum systems

To create quasi-static scans for a system that evolves according to a master equation, you need to employ the graph framework to obtain the points that compose the quasi-static scan. Assuming a noise process $\eta$ which changes sufficiently slowly that it is constant in each run of the gate, each point represents the performance of the pulse under different strengths of $\eta$ in combination with any open-systems dynamics.

To quantify the performance of the target control, you can use the state infidelity as a metric. Using the fact that the target state is a pure state $\rho_{\rm target}$, you can find the infidelity of a pulse that evolves the system to a state $\rho(t)$ through the formula:

$$ \mathcal{I} (t) = 1 - \left( \mathrm{Tr} \left\{ \sqrt{ \sqrt{\rho_{\rm target}} \rho(t) \sqrt{\rho_{\rm target}} } \right\} \right)^2 = 1 - \mathrm{Tr} \left\{ \rho_{\rm target} \rho(t) \right\}. $$

For this example, consider a system under the following system Hamiltonian:

$$ H_{\rm s}(t) = \frac{1}{2} \Omega(t) \sigma_- + \frac{1}{2} \Omega^* (t) \sigma_+ + \eta \sigma_z, $$

where $\sigma_x$, $\sigma_y$, and $\sigma_z$ are the Pauli matrices and $\sigma_\pm = (\sigma_x \mp i \sigma_y)/2$. The signal $\Omega(t)$ represents your complex-valued controls, while $\eta$ is the noise coefficient.

Consider as well that the system can longitudinally relax at a certain rate $T_1$, so that the complete GKS–Lindblad master equation that describes its evolution is:

$$ \frac{\mathrm{d}}{\mathrm{d} t} \rho(t) = - i \left[ H_{\rm s}(t), \rho(t) \right] + \frac{1}{T_1} \sigma_- \rho(t) \sigma_+ - \frac{1}{2} \frac{1}{T_1} \left\{ \rho(t), \sigma_+ \sigma_- \right\}. $$

In this master equation, $\sigma_-$ is the Lindblad operator with an associated rate $1/T_1$.

To see how a set of pulses performs against additive quasistatic $\sigma_z$ noise, when also subject to incoherent $T_1$ processes, you can use the open system tools to calculate the evolution of the master equation for different values of $\eta$. This example shows how you can create an optimized pulse under these conditions by setting the sum of the infidelities under different quasi-static strengths as the cost function of the optimization. The code block below compares the performance of the predefined pulse and the optimized pulse, plotting the shape of the controls and the quasi-static scans for both of them.

# Define standard matrices.
identity = np.identity(2, dtype=complex)
sigma_x = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=complex)
sigma_z = np.array([[1.0, 0.0], [0.0, -1.0]], dtype=complex)
sigma_m = np.array([[0.0, 1.0], [0.0, 0.0]], dtype=complex)

# Define noise and control parameters.
T1 = 100e-6  # s
omega_max = 2 * np.pi * 1e6  # Hz
total_rotation = np.pi / 2
target_operator = (identity - 1j * sigma_x) / np.sqrt(2)
initial_density_matrix = np.array([[1, 0], [0, 0]], dtype=complex)
target_density_matrix = (
    target_operator @ initial_density_matrix @ target_operator.T.conj()
)

# Define scan parameters
scan_point_count = 20
scan_range = np.linspace(-omega_max / 2, omega_max / 2, scan_point_count)

# Obtain predefined pulse from Q-CTRL Open Controls.
corpse_pulse = qctrlopencontrols.new_corpse_control(
    rabi_rotation=total_rotation, azimuthal_angle=0.0, maximum_rabi_rate=omega_max
)
duration = sum(corpse_pulse.durations)


def state_infidelity(graph, density_matrix, name):
    """
    Calculates the state infidelity of the `density_matrix` with respect to the
    `target_density_matrix`, which represents a pure state.
    """
    return graph.abs(1 - graph.trace(target_density_matrix @ density_matrix), name=name)


# Create optimization graph.
graph = qctrl.create_graph()

# Define the controls for the predefined pulse.
corpse_omega = graph.pwc(
    values=corpse_pulse.rabi_rates * np.exp(1j * corpse_pulse.azimuthal_angles),
    durations=corpse_pulse.durations,
    name="corpse_pulse",
)

# Define controls for the optimized pulse.
optimized_omega = graph.complex_pwc_signal(
    moduli=graph.optimization_variable(count=10, lower_bound=0, upper_bound=omega_max),
    phases=graph.optimization_variable(
        count=10,
        lower_bound=0,
        upper_bound=2 * np.pi,
        is_lower_unbounded=True,
        is_upper_unbounded=True,
    ),
    duration=duration,
    name="optimized_pulse",
)

# Create control Hamiltonians for the predefined pulse and the optimized pulse.
noiseless_corpse_hamiltonian = graph.pwc_operator_hermitian_part(corpse_omega * sigma_m)
noiseless_optimized_hamiltonian = graph.pwc_operator_hermitian_part(
    optimized_omega * sigma_m
)

# Create list of detuning terms for the quasi-static scan.
scan_terms = [amplitude * sigma_z for amplitude in scan_range]

# Add detuning terms to the control Hamiltonians.
corpse_hamiltonians = [
    noiseless_corpse_hamiltonian + scan_term for scan_term in scan_terms
]
optimized_hamiltonians = [
    noiseless_optimized_hamiltonian + scan_term for scan_term in scan_terms
]

# Define decay term of the master equation.
lindblad_operator = sigma_m

# Calculate the final density matrices from the master equation.
corpse_density_matrices = [
    graph.density_matrix_evolution_pwc(
        initial_density_matrix=initial_density_matrix,
        hamiltonian=corpse_hamiltonian,
        lindblad_terms=[(1 / T1, lindblad_operator)],
        sample_times=np.array([duration]),
    )[-1]
    for corpse_hamiltonian in corpse_hamiltonians
]
optimized_density_matrices = [
    graph.density_matrix_evolution_pwc(
        initial_density_matrix=initial_density_matrix,
        hamiltonian=optimized_hamiltonian,
        lindblad_terms=[(1 / T1, lindblad_operator)],
        sample_times=np.array([duration]),
    )[-1]
    for optimized_hamiltonian in optimized_hamiltonians
]

# Calculate infidelities with respect to the target.
corpse_infidelities = [
    state_infidelity(graph, density_matrix, name=f"corpse_infidelity{number}")
    for number, density_matrix in enumerate(corpse_density_matrices)
]
optimized_infidelities = [
    state_infidelity(graph, density_matrix, name=f"optimized_infidelity{number}")
    for number, density_matrix in enumerate(optimized_density_matrices)
]

# Create cost as the sum of all the infidelities of the optimized pulse.
graph.sum(optimized_infidelities, name="cost")

# Run optimization.
results = qctrl.functions.calculate_optimization(
    graph=graph,
    cost_node_name="cost",
    output_node_names=["corpse_pulse", "optimized_pulse"]
    + [f"corpse_infidelity{number}" for number in range(scan_point_count)]
    + [f"optimized_infidelity{number}" for number in range(scan_point_count)],
)

# Plot pulses.
fig1 = plt.figure()
fig1.suptitle("CORPSE pulse")
plot_controls(fig1, {r"$\Omega$": results.output[f"corpse_pulse"]})

fig2 = plt.figure()
fig2.suptitle("Optimized pulse")
plot_controls(fig2, {r"$\Omega$": results.output[f"optimized_pulse"]})

# Plot quasi-static scan.
fig3 = plt.figure(figsize=(10, 5))
plt.xlabel(r"Detuning $\eta$/(2$\pi$) (MHz)")
plt.ylabel("Infidelity")
plt.plot(
    scan_range * 1e-6 / (2 * np.pi),
    [
        results.output[f"corpse_infidelity{number}"]["value"]
        for number in range(scan_point_count)
    ],
    label="CORPSE pulse",
)
plt.plot(
    scan_range * 1e-6 / (2 * np.pi),
    [
        results.output[f"optimized_infidelity{number}"]["value"]
        for number in range(scan_point_count)
    ],
    label="Optimized pulse",
)
plt.legend()
plt.show()
Your task calculate_optimization (action_id="625829") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_optimization (action_id="625829") has completed.