Superconducting qubits: improving the performance of the cross resonance gate

Increasing robustness against crosstalk

BOULDER OPAL enables you to optimize controls to achieve a target operation on your quantum hardware. In this Application note, you'll find the steps needed to obtain robust control solutions for superconducting qubits coupled by a common transmission line. These controls exhibit increased robustness of the cross resonance gate (target) against noise that arises due to crosstalk. First, you'll learn to set up the system and to optimize the control pulses. You can then verify the optimized controls through simulation of the resulting time evolution. Finally, you can examine the robustness of the controls using a quasi-static scan.

This notebook contains preview features that are not currently available in BOULDER OPAL. Please contact us if you want a live demonstration for your hardware. Preview features are included through the package qctrlcore, which is not publicly available: cells with qctrlcore will not be executable.

Imports and initialization

# Essential imports
import numpy as np
import tensorflow as tf
import qctrlcore as qc
from scipy.linalg import expm
from tensorflow.python.util import deprecation
deprecation._PRINT_DEPRECATION_WARNINGS = False

# Plotting imports
import matplotlib.pyplot as plt
from qctrlvisualizer import plot_controls
from qctrlvisualizer import display_two_qubit_visualization

# Helper functions
def average_gate_infidelity(final_unitary: np.array,
                            target: qc.Target) -> float:
    """Calculates the measure of system performance, using the final and target gates."""
    subspace_target = target.projector.dot(
        target.operator).dot(target.projector)
    subspace_dimension = np.sum(np.abs(subspace_target)**2)

    op = np.matmul(
        subspace_target.T.conj(),
        final_unitary)

    infidelity = np.abs(
        1 -
        (np.abs(np.sum(op.conj() * op)) +
         np.abs(np.sum(subspace_target.conj() * final_unitary))**2) /
        ((1 + subspace_dimension) * subspace_dimension))

    return infidelity

def crosstalk_4(constant, result):
    """Calculates the additional drives due to crosstalk."""
    return [
        {'operator': drive4_2x_operator,
         'control': [qc.ComplexSegment(d['duration'], constant*d['value']) for d in
                     result['output']['epsilon_1x']]},
        {'operator': drive4_2y_operator,
         'control': [qc.ComplexSegment(d['duration'], constant*d['value']) for d in
                     result['output']['epsilon_1y']]},
    ]

def get_robustness_scan_4(detunings, scheme_name):
    """Calculates the operational infidelity of a scheme, at specific noise."""
    infidelities = []
    average_gate_infidelities = []

    if scheme_name is 'Robust':
        phase_shift_operators = phase_shift_op_dic['phase_shift_operators_robust_4']
        result = result_robust_4
        shifts = controls_dic['robust_controls_4']
    elif scheme_name is 'Optimal':
        phase_shift_operators = phase_shift_op_dic['phase_shift_operators_optimal_4']
        result = result_optimal_4
        shifts = controls_dic['optimal_controls_4']
    elif scheme_name is 'Primitive':
        phase_shift_operators = phase_shift_op_dic['phase_shift_operators_primitive_4']
        result = result_primitive_4
        shifts = controls_dic['primitive_controls_4']

    for detuning in detunings:
        simulation_result = qc.calculate_unitary_evolution(
            duration=gate_duration_cr,
            shifts=shifts,
            drifts=drifts_4,
            drives=crosstalk_4(detuning, result))
        final_unitary_with_phase_shift = phase_shift_operators[1].dot(
            simulation_result[-1]['evolution_operator'].operator.dot(
                phase_shift_operators[0]))
        infidelities.append(
            1 - qc.calculate_operational_fidelity(
                final_unitary_with_phase_shift,
                zx_half_4*10/4)['fidelity'])
        average_gate_infidelities.append(
            average_gate_infidelity(
                final_unitary_with_phase_shift,
                qc.Target(zx_half_4)))

    return [np.array(infidelities), np.array(average_gate_infidelities)]

colors = {"Optimal": "#BF04DC", "Robust": "#680CE9", "Primitive": "k"}

Creating a Q-CTRL cross resonance gate robust to crosstalk

Consider a system of two transmon qubits coupled by a common transmission line. The effective coupling between the qubits can be calculated by a Schrieffer-Wolff transformation (up to second order). Qubit 1 is defined as the control qubit, and qubit 2 is the target qubit. The control qubit is driven by the in-phase and quadrature microwave drive with a frequency that matches the target qubit frequency. In a frame rotating with the drive frequency and after applying the rotating wave approximation (RWA) the Hamiltonian is of the form:

$$ H = \hbar(\Delta + \epsilon)n_1+ \frac{\hbar\delta_1}{2} n_1(n_1 - 1)+ \hbar\epsilon n_2+ \frac{\hbar\delta_2}{2}n_2(n_2 - 1)+ \hbar g \left( b_1^\dagger b_2 + b_1 b_2^\dagger \right) +\frac{\hbar}{2}\left(\Omega(t)b_1^\dagger + \text{h.c.}\right)\,,$$

where $n_i$ and $b_i$ are the number and annihilation operators of the $i$th qubit, $\Delta=\omega_1-\omega_2$ is the detuning between the two qubits where $\omega_i$ is the frequency of the $i$th qubit, $\epsilon$ is the mismatch between the drive frequency and the target qubit frequency, $g$ is the effective coupling strength between the qubits, $\Omega$ is the drive modulus and $\delta_i$ is the anharmonicity of the $i$th qubit. The above Hamiltonian is valid when the coupling strength is small compared to the detuning between the two qubits, $\Delta\gg g$.

The target operation is the $\pi/2$ cross resonance gate which operates in the space of the two lowest levels:

$$U_{CR}\left(\frac{\pi}{2}\right)=\frac{1}{\sqrt{2}}\left(\mathbb{I}-i\sigma_z^{(1)}\sigma_x^{(2)}\right)=[ZX]^{1/2}.$$

Since the above operation is defined up to phase factors of the individual qubits, this may lead to a large infidelity. Consequently, a final (optimizable) phase shift on each qubit is introduced to overcome the infidelity inflation due to these phase factors.

This Application note addresses robustness to crosstalk. Crosstalk between the qubits leads to unwanted effects that may mimic the effect of noise. The crosstalk Hamiltonian has the form

$$ H_\text{ct} = c_\text{ct} \frac{\hbar\Omega(t)}{2}b_2^\dagger + \text{h.c.}\,, $$

where $c_\text{ct}$ is the (complex) crosstalk constant.

In this section, you'll learn to optimize the implementation of the target gate under the effect of crosstalk. This is performed within the framework of the four-level approximation, where only the first four levels of each qubit are considered. Moreover, only up to three excitations overall are permitted, which further reduces the Hilbert space from $16$ states to $10$ states.

To set up the optimization, first define operators and constants.

Mega = 1e6
micro = 1e-6

# Constants in the lab frame
omega_1 = 5114 * 2 * np.pi # MHz
omega_2 = 4914 * 2 * np.pi # MHz
delta_1 = -330 * 2 * np.pi # MHz
delta_2 = -330 * 2 * np.pi # MHz
g = 3.8 * 2 * np.pi # MHz
Delta = omega_1 - omega_2
epsilon = 0 # g**2/Delta

# Limits for drive amplitudes
epsilon_1_min = -2*np.pi * 200 # MHz
epsilon_1_max = 2*np.pi * 200 # MHz

# These redefined constants are required due to the Schrieffer-Wolff transformation
C1 = 1/(1 + g**2/Delta**2)
C202 = 1/(1 + 2*g**2/(Delta - delta_2)**2)
C211 = 1/(1 + 2*g**2/(Delta - delta_2)**2 + 2*g**2/(Delta + delta_1)**2)
A220 = 2*g**2/((Delta - delta_2)**2 + 2*g**2) * \
    (Delta - delta_2)/(Delta + delta_1)
B220 = np.sqrt(2)*g/((Delta - delta_2)**2 + 2*g**2) * \
    (Delta - delta_2)**2/(Delta + delta_1)
C220 = 1/(1 + A220**2 + B220**2)

Next, define the operators in the Hamiltonian and the target gate. Each operator is defined three times, for each level of approximation: two-level, three-level and four-level.

# Lowering operator
b4 = np.diag([1, np.sqrt(2), np.sqrt(3)], 1)
# Number operator
n4 = np.diag([0, 1, 2, 3])
# Projector onto the two lowest states
qubit_proj_4 = np.diag([1, 1, 0, 0])

# Pauli operators
sigma_x = np.array([[0, 1], [1, 0]])
sigma4_x = np.zeros(n4.shape)
sigma4_x[:sigma_x.shape[0], :sigma_x.shape[1]] = sigma_x

sigma_y = 1j*np.array([[0, -1], [1, 0]])
sigma4_y = np.zeros(n4.shape, dtype=np.complex)
sigma4_y[:sigma_y.shape[0], :sigma_y.shape[1]] = sigma_y

sigma_z = np.array([[1, 0], [0, -1]])
sigma4_z = np.zeros(n4.shape)
sigma4_z[:sigma_z.shape[0], :sigma_z.shape[1]] = sigma_z

# Embedding the above single-qubit operators into the two-qubit system of dimension 16
b41 = np.kron(b4, np.eye(4))
b42 = np.kron(np.eye(4), b4)

n41 = np.kron(n4, np.eye(4))
n42 = np.kron(np.eye(4), n4)

proj_4 = np.kron(qubit_proj_4, qubit_proj_4)

sigma41_x = np.kron(sigma4_x, np.eye(4))
sigma42_x = np.kron(np.eye(4), sigma4_x)

sigma41_y = np.kron(sigma4_y, np.eye(4))
sigma42_y = np.kron(np.eye(4), sigma4_y)

sigma41_z = np.kron(sigma4_z, np.eye(4))
sigma42_z = np.kron(np.eye(4), sigma4_z)

# Block-diagonal permutations
perm4 = np.argsort(np.diag(n41 + n42))

# Approximate eigenstates in block-form
e0 = np.array([[1]])
e1 = np.array([
    [np.sqrt(C1), g/Delta*np.sqrt(C1)],
    [-g/Delta*np.sqrt(C1), np.sqrt(C1)]])
e2 = np.array([
    [np.sqrt(C202), np.sqrt(2)*g/(Delta - delta_2)
     * np.sqrt(C211), A220*np.sqrt(C220)],
    [-np.sqrt(2)*g/(Delta - delta_2)*np.sqrt(C202),
     np.sqrt(C211), B220*np.sqrt(C220)],
    [0, -np.sqrt(2)*g/(Delta + delta_1)*np.sqrt(C211), np.sqrt(C220)]])

# Frame transformation HS = SD
S4 = np.eye(16)
S4[1:1+e1.shape[0], 1:1+e1.shape[1]] = e1
S4[1+e1.shape[0]:1+e1.shape[0]+e2.shape[0], 1 +
    e1.shape[1]:1+e1.shape[1]+e2.shape[1]] = e2

# The phase shift operators in the simulation/optimization frame
phase_shift_operator41 = (
    S4.dot(sigma41_z[:, perm4][perm4, :]/2)).dot(S4.T.conj())[:10, :10]
phase_shift_operator42 = (
    S4.dot(sigma42_z[:, perm4][perm4, :]/2)).dot(S4.T.conj())[:10, :10]

# Target in the simulation/optimization frame
zx_half_4 = (S4.dot(
    ((np.kron(qubit_proj_4, qubit_proj_4)
      - 1j*np.kron(sigma4_z, sigma4_x))/np.sqrt(2))[:, perm4][perm4, :]).dot(S4.T.conj()))[:10, :10]

Creating control pulses

Next, define the optimization configuration and the control schemes: primitive, optimal and robust.

# Optimization constraints
gate_duration_cr = 0.4 # us
number_of_segments_primitive = 1
number_of_segments_optimal = 200
crosstalk_constant = 0.1*(1 + 1j)

pulse_durations_primitive = np.ones(
    [number_of_segments_primitive])*gate_duration_cr/number_of_segments_primitive
pulse_durations_optimal = np.ones(
    [number_of_segments_optimal])*gate_duration_cr/number_of_segments_optimal

# Gate schemes
detuning_array = np.linspace(-1, 1, 101)
gate_schemes = ['Optimal', 'Robust', 'Primitive']

The optimization specification in the next cell involves several steps:

  • create the (noiseless) Hamiltonian,
  • create the noise operators,
  • define the cost function to be minimized by the optimizer, which is the infidelity of the system's evolution relative to an iSWAP target,
  • introduce the final (optimizable) phase shift on each qubit,
  • compute the infidelity of the Hamiltonian with optimizable phase shifts.

The outputs are the optimized values of $\Omega(t)$ ($x$ and $y$ components) across the gate duration, the trajectory of the coupling qubit frequency, the optimized phase shifts, and the infidelities with and without the phase shifts.

def create_optimization_specification_cr(variable_factory):
    epsilon_1x = variable_factory.create_bounded(
        count=number_of_segments_cr,
        lower_bound=epsilon_1_min,
        upper_bound=epsilon_1_max)
    if primitive_flag:
        epsilon_1y = variable_factory.create_bounded(
            count=number_of_segments_cr,
            lower_bound=epsilon_1_min*1e-6,
            upper_bound=epsilon_1_max*1e-6)
    else:
        epsilon_1y = variable_factory.create_bounded(
            count=number_of_segments_cr,
            lower_bound=epsilon_1_min,
            upper_bound=epsilon_1_max)

    # Build the system Hamiltonian for each level of approximation
    drive_1x = qc.create_pwc_signal(
        values=epsilon_1x,
        duration=gate_duration_cr)
    drive_1y = qc.create_pwc_signal(
        values=epsilon_1y,
        duration=gate_duration_cr)

    if level_count == 4:
        drive_1x_term = qc.create_pwc_operator(
            signal=drive_1x,
            operator=qc.HermitianOperator((b41.T.conj() + b41)[:, perm4[:10]][perm4[:10], :]/2))
        drive_1y_term = qc.create_pwc_operator(
            signal=drive_1y,
            operator=qc.HermitianOperator(1j*(b41.T.conj() - b41)[:, perm4[:10]][perm4[:10], :]/2))

        fixed_frequency_1_term = qc.create_constant_pwc_operator(
            operator=qc.HermitianOperator(
                (Delta - delta_1/2 + epsilon)*n41[:, perm4[:10]][perm4[:10], :]),
            duration=gate_duration_cr)
        fixed_frequency_2_term = qc.create_constant_pwc_operator(
            operator=qc.HermitianOperator(
                (-delta_2/2 + epsilon)*n42[:, perm4[:10]][perm4[:10], :]),
            duration=gate_duration_cr)
        fixed_anharmonicity_1_term = qc.create_constant_pwc_operator(
            operator=qc.HermitianOperator(
                delta_1*np.matmul(n41, n41)[:, perm4[:10]][perm4[:10], :]/2),
            duration=gate_duration_cr)
        fixed_anharmonicity_2_term = qc.create_constant_pwc_operator(
            operator=qc.HermitianOperator(
                delta_2*np.matmul(n42, n42)[:, perm4[:10]][perm4[:10], :]/2),
            duration=gate_duration_cr)
        fixed_coupling_term = qc.create_constant_pwc_operator(
            operator=qc.HermitianOperator(g*(
                np.kron(b4.T.conj(), b4) + np.kron(b4, b4.T.conj()))[:, perm4[:10]][perm4[:10], :]),
            duration=gate_duration_cr)

    # Create noise operators
    if level_count == 4:
        crosstalk_1x = qc.create_pwc_operator(
            signal=drive_1x,
            operator=qc.NonHermitianOperator(crosstalk_constant*b42.T.conj()[:, perm4[:10]][perm4[:10], :]/2))
        crosstalk_1y = qc.create_pwc_operator(
            signal=drive_1y,
            operator=qc.NonHermitianOperator(1j*crosstalk_constant*b42.T.conj()[:, perm4[:10]][perm4[:10], :]/2))

    # Construct the total Hamiltonian
        hamiltonian = qc.get_pwc_operator_sum(
            [drive_1x_term,
             drive_1y_term,
             fixed_frequency_1_term,
             fixed_frequency_2_term,
             fixed_anharmonicity_1_term,
             fixed_anharmonicity_2_term,
             fixed_coupling_term,
             ])

    if noise_flag:
        print("Noise optimization")
        noise_list = [crosstalk_1x, crosstalk_1y]
    else:
        noise_list = []

    # Define the cost function
    if level_count == 4:
        infidelity = qc.create_infidelity(
            hamiltonian=hamiltonian,
            noise_operators=noise_list,
            target=qc.Target(
                zx_half_4))

    # Introduce optimizable phase shifts on each qubit
    phase_shift_1 = variable_factory.create_bounded(
        count=1,
        lower_bound=0,
        upper_bound=2*np.pi)
    phase_shift_2 = variable_factory.create_bounded(
        count=1,
        lower_bound=0,
        upper_bound=2*np.pi)

    if level_count == 4:
        phase_shift_hamiltonian_1 = qc.create_pwc_operator(
            operator=qc.HermitianOperator(-1e6*phase_shift_operator42),
            signal=qc.create_pwc_signal(phase_shift_2, 1e-6))
        phase_shift_hamiltonian_2 = qc.get_pwc_operator_sum(
            [qc.create_pwc_operator(operator=qc.HermitianOperator(1e6*phase_shift_operator41),
                                    signal=qc.create_pwc_signal(phase_shift_1, 1e-6)),
             qc.create_pwc_operator(operator=qc.HermitianOperator(1e6*phase_shift_operator42),
                                    signal=qc.create_pwc_signal(phase_shift_2, 1e-6)),
             ])

    hamiltonian_with_phase_shift = qc.TensorPwc(
        durations=np.concatenate(
            [phase_shift_hamiltonian_1.durations,
             hamiltonian.durations,
             phase_shift_hamiltonian_2.durations]),
        values=tf.concat(
            [phase_shift_hamiltonian_1.values,
             hamiltonian.values,
             phase_shift_hamiltonian_2.values], axis=0))

    # Compute the infidelity of the Hamiltonian with optimizable phase shifts
    if level_count == 4:
        infidelity_with_phase_shift = qc.create_infidelity(
            hamiltonian=hamiltonian_with_phase_shift,
            noise_operators=noise_list,
            target=qc.Target(
                zx_half_4))

    # Define output quantities
    def output_callable(session): return {
        'epsilon_1x': [{'duration': d, 'value': v}
                       for d, v in zip(drive_1x.durations,
                                       session.run(epsilon_1x))],
        'epsilon_1y': [{'duration': d, 'value': v}
                       for d, v in zip(drive_1y.durations,
                                       session.run(epsilon_1y))],
        'phase_shift_1': session.run(phase_shift_1[0]),
        'phase_shift_2': session.run(phase_shift_2[0]),
        'infidelity': session.run(infidelity),
        'infidelity_with_phase_shift': session.run(infidelity_with_phase_shift),
        'final_unitary': session.run(
            qc.tfutil.chain_matmul_sequential(
                qc.cost._create_segment_unitaries_from_hamiltonian(hamiltonian), reverse=True)),
        'final_unitary_with_phase_shift': session.run(
            qc.tfutil.chain_matmul_sequential(
                qc.cost._create_segment_unitaries_from_hamiltonian(hamiltonian_with_phase_shift), reverse=True)),
    }

    return qc.OptimizationSpecification(
        cost=infidelity_with_phase_shift,
        output_callable=output_callable)

Next, perform the optimization within the three-excitation (four-level) approximation. The following cell optimizes a primitive pulse as a benchmark, along with the optimal and the robust controls. The infidelities for each control scheme are presented with and without the phase shift.

level_count = 4
noise_flag = False

print('Primitive')
primitive_flag = True
number_of_segments_cr = number_of_segments_primitive
result_primitive_4 = qc.optimize(create_optimization_specification_cr, number_of_optimizations=40)
print("Cost: " + str(result_primitive_4['cost']))
print("Infidelity: " + str(result_primitive_4['output']['infidelity']))
print("Infidelity with phase shift: " + str(result_primitive_4['output']['infidelity_with_phase_shift']))

level_count = 4
noise_flag = False

print('\nOptimal')
primitive_flag = False
number_of_segments_cr = number_of_segments_optimal
result_optimal_4 = qc.optimize(create_optimization_specification_cr, number_of_optimizations=15)
print("Cost: " + str(result_optimal_4['cost']))
print("Infidelity: " + str(result_optimal_4['output']['infidelity']))
print("Infidelity with phase shift: " + str(result_optimal_4['output']['infidelity_with_phase_shift']))

level_count = 4
noise_flag = True

print('\nRobust')
number_of_segments_cr = number_of_segments_optimal
result_robust_4 = qc.optimize(create_optimization_specification_cr, number_of_optimizations=15)
print("Cost: " + str(result_robust_4['cost']))
print("Infidelity: " + str(result_robust_4['output']['infidelity']))
print("Infidelity with phase shift: " + str(result_robust_4['output']['infidelity_with_phase_shift']))
Primitive
Cost: 0.1139598968360277
Infidelity: 0.9377796375228804
Infidelity with phase shift: 0.1139598968360277

Optimal
Cost: 2.0452393290071313e-07
Infidelity: 0.07832160306115121
Infidelity with phase shift: 2.0452393290071313e-07

Robust
Noise optimization
Cost: 0.003352485653389423
Infidelity: 0.8207073789265535
Infidelity with phase shift: 0.003352485653389423

Even in the presence of crosstalk, the infidelity arising from the robust pulse (with phase shifts) is extremely low. Visualize the control schemes using the cell below.

# Primitive
omega_1x_primitive_4 = [
    {'duration': d['duration']*micro, 'value': d['value']*Mega}
    for d in result_primitive_4['output']['epsilon_1x']]
omega_1y_primitive_4 = [
    {'duration': d['duration']*micro, 'value': d['value']*Mega}
    for d in result_primitive_4['output']['epsilon_1y']]

plot_controls(plt.figure(), {"$\\Omega_{rx}$": omega_1x_primitive_4,
                             "$\\Omega_{ry}$": omega_1y_primitive_4,
                            })
# Optimized
omega_1x_optimal_4 = [
    {'duration': d['duration']*micro, 'value': d['value']*Mega}
    for d in result_optimal_4['output']['epsilon_1x']]
omega_1y_optimal_4 = [
    {'duration': d['duration']*micro, 'value': d['value']*Mega}
    for d in result_optimal_4['output']['epsilon_1y']]
plot_controls(plt.figure(), {"$\\Omega_{ox}$": omega_1x_optimal_4,
                             "$\\Omega_{oy}$": omega_1y_optimal_4,
                            })
# Robust
omega_1x_robust_4 = [
    {'duration': d['duration']*micro, 'value': d['value']*Mega}
    for d in result_robust_4['output']['epsilon_1x']]
omega_1y_robust_4 = [
    {'duration': d['duration']*micro, 'value': d['value']*Mega}
    for d in result_robust_4['output']['epsilon_1y']]
plot_controls(plt.figure(), {"$\\Omega_{rx}$": result_robust_4['output']['epsilon_1x'],
                             "$\\Omega_{ry}$": result_robust_4['output']['epsilon_1y'],
                            })

The above plots display the control pulses ($\Omega_x$ and $\Omega_y$) for the primitive, optimal and robust schemes, respectively.

Verifying $U_{CR}\left(\frac{\pi}{2}\right)$ operation through simulation

Use the optimized controls obtained above to calculate the exact time evolution determined by the Hamiltonian in the four-level approximation. In the cell below, the unitary evolution operators are obtained for each scheme. These are first applied to calculate and output the operational infidelities. Additionally, the evolution operators are used to calculate the evolution of an equal superposition of the four lowest 2-qubit states, $|\psi_0\rangle =\frac{1}{2}\left(|00\rangle-i|01\rangle+i|10\rangle+|11\rangle\right)$. The final state after an ideal operation of the CR gate should be the Bell-state $|\psi_f\rangle =\frac{1}{2}\left(|01\rangle+|10\rangle\right)$. The populations are displayed to assess each scheme.

drifts_4 = [
    {'operator': qc.HermitianOperator(
        (Delta - delta_1/2 + epsilon)*n41[:, perm4[:10]][perm4[:10], :])},
    {'operator': qc.HermitianOperator(
        (-delta_2/2 + epsilon)*n42[:, perm4[:10]][perm4[:10], :])},
    {'operator': qc.HermitianOperator(
        delta_1*np.matmul(n41, n41)[:, perm4[:10]][perm4[:10], :]/2)},
    {'operator': qc.HermitianOperator(
        delta_2*np.matmul(n42, n42)[:, perm4[:10]][perm4[:10], :]/2)},
    {'operator': qc.HermitianOperator(g*(
        np.kron(b4.T.conj(), b4) + np.kron(b4, b4.T.conj()))[:, perm4[:10]][perm4[:10], :])},
]
drive4_1x_operator = qc.HermitianOperator(
    (b41.T.conj() + b41)[:, perm4[:10]][perm4[:10], :]/2)
drive4_1y_operator = qc.HermitianOperator(
    1j*(b41.T.conj() - b41)[:, perm4[:10]][perm4[:10], :]/2)

schemes_results = {'primitive': result_primitive_4,
                   'optimal': result_optimal_4, 'robust': result_robust_4}
phase_shift_op_dic = {}
controls_dic = {}
simulation_dic = {}
final_unitary_dic = {}
final_unitary_with_phase_shift_dic = {}

for scheme in ['primitive', 'optimal', 'robust']:
    phase_shift_op_dic['phase_shift_operators_'+scheme+'_4'] = [
        expm(1j*schemes_results[scheme]['output']
             ['phase_shift_2']*phase_shift_operator42),
        expm(-1j*schemes_results[scheme]['output']['phase_shift_1'] * phase_shift_operator41
             - 1j*schemes_results[scheme]['output']['phase_shift_2']*phase_shift_operator42)]

    controls_dic[scheme+'_controls_4'] = [
        {'operator': drive4_1x_operator,
         'control': [qc.RealSegment(d['duration'], d['value']) for d in
                     schemes_results[scheme]['output']['epsilon_1x']]},
        {'operator': drive4_1y_operator,
         'control': [qc.RealSegment(d['duration'], d['value']) for d in
                     schemes_results[scheme]['output']['epsilon_1y']]},
    ]

    simulation_dic['simulation'+scheme+'4'] = qc.calculate_unitary_evolution(
        duration=gate_duration_cr,
        shifts=controls_dic[scheme+'_controls_4'],
        drifts=drifts_4)

    final_unitary_dic['final_unitary_'+scheme+'_4'] = simulation_dic['simulation' +
                                                                     scheme+'4'][-1]['evolution_operator'].operator
    final_unitary_with_phase_shift_dic['final_unitary_'+scheme + '_with_phase_shift_4'] = phase_shift_op_dic['phase_shift_operators_'+scheme+'_4'][1].dot(
        final_unitary_dic['final_unitary_'+scheme+'_4'].dot(
            phase_shift_op_dic['phase_shift_operators_'+scheme+'_4'][0]))

infidelities_4 = [
    1-qc.calculate_operational_fidelity(
        final_unitary_with_phase_shift_dic['final_unitary_primitive_with_phase_shift_4'],
        zx_half_4*10/4)['fidelity'],
    1-qc.calculate_operational_fidelity(
        final_unitary_with_phase_shift_dic['final_unitary_optimal_with_phase_shift_4'],
        zx_half_4*10/4)['fidelity'],
    1-qc.calculate_operational_fidelity(
        final_unitary_with_phase_shift_dic['final_unitary_robust_with_phase_shift_4'],
        zx_half_4*10/4)['fidelity'], ]

print('Primitive\t\t Optimal\t\t\t Robust')
print(infidelities_4[0], '\t', infidelities_4[1], '\t', infidelities_4[2])

initial_state = np.array([1, -1j, 1j, 0, 1, 0, 0, 0, 0, 0])/2
Prob_robust = (
    np.abs(np.dot(final_unitary_with_phase_shift_dic['final_unitary_robust_with_phase_shift_4'], initial_state)))**2
Prob_optimal = (
    np.abs(np.dot(final_unitary_with_phase_shift_dic['final_unitary_optimal_with_phase_shift_4'], initial_state)))**2
Prob_primitive = (
    np.abs(np.dot(final_unitary_with_phase_shift_dic['final_unitary_primitive_with_phase_shift_4'], initial_state)))**2


fig = plt.figure(figsize=(14, 3))
axes = fig.subplots(nrows=1, ncols=3, sharex=False,
                    sharey=True, squeeze=False)[0, :]

langs = ['|01>', '|10>', 'Rest']
population_robust = [Prob_robust[1],
                     Prob_robust[2], 1-Prob_robust[1]-Prob_robust[2]]
axes[0].bar(langs, population_robust,color=colors['Robust'])
axes[0].set_title("Robust", fontsize=18)
axes[0].set_xlabel("States", fontsize=14)
axes[0].set_ylabel("Probability", fontsize=14)

population_optimal = [Prob_optimal[1],
                        Prob_optimal[2], 1-Prob_optimal[1]-Prob_optimal[2]]
axes[1].bar(langs, population_optimal,color=colors['Optimal'])
axes[1].set_title("Optimal", fontsize=18)


population_primitive = [Prob_primitive[1],
                        Prob_primitive[2], 1-Prob_primitive[1]-Prob_primitive[2]]
axes[2].bar(langs, population_primitive,color=colors['Primitive'])
axes[2].set_title("Primitive", fontsize=18)

plt.show()
Primitive		 Optimal			 Robust
0.11395989683593 	 2.0452397042625137e-07 	 0.00016424352656441155

Clearly, there is a non-negligible leakage in the primitive scheme, while both the optimal and robust schemes perform almost perfectly (up to a small asymmetry).

Verifying robustness and accuracy of Q-CTRL cross resonance gate

You can assess the robustness of the different control solutions using a quasi-static scan. Here the scan is over the crosstalk strength.

drive4_2x_operator = qc.NonHermitianOperator(crosstalk_constant*b42.T.conj()[:, perm4[:10]][perm4[:10], :]/2)
drive4_2y_operator = qc.NonHermitianOperator(1j*crosstalk_constant*b42.T.conj()[:, perm4[:10]][perm4[:10], :]/2)

gate_infidelity_4 = {scheme_name: get_robustness_scan_4(detuning_array, scheme_name)
                     for scheme_name in gate_schemes}
fig = plt.figure(figsize=(12,3))
axes = fig.subplots(nrows=1, ncols=1, sharex=False, sharey=True, squeeze=False)[0, :]

for name, infidelity in gate_infidelity_4.items():
    axes[0].plot(detuning_array, infidelity[0], label=name,color=colors[name])
    axes[0].set_xlabel("Relative crosstalk", fontsize=14)
    axes[0].set_ylabel("Infidelities", fontsize=14)
    axes[0].ticklabel_format(style='sci', axis='x', scilimits=(0,0))
    axes[0].legend(loc='best', bbox_to_anchor=(1, 1.0))
plt.show()

The figure demonstrates the different noise (crosstalk) susceptibilities (in terms of operational infidelity) of the primitive, optimal and robust controls. The robust scheme has a much broader low-infidelity region.