Designing robust, configurable, parallel gates for large trapped-ion arrays

Obtaining control solutions for parallel and specifiable multi-qubit gates using Q-CTRL pulses

BOULDER OPAL provides a broad toolkit enabling the generation of high-fidelity entangling gates with trapped ions. In particular, it is possible to engineer parallel optimized gates with arbitrary user-defined phases on large ion arrays. This comes from adding independent addressing and individual tuning of control parameters in order to simultaneously manipulate the phase-space trajectory of multiple motional modes.

In this notebook, we demonstrate control solutions for trapped ions robust to frequency-detuning and pulse-timing errors. We will cover:

  • Creating optimized, robust, parallel Mølmer–Sørensen gates over multiple ions
  • Calculating phase space trajectories for multiple modes
  • Simulating entangling phase accumulation in the gate between arbitrary ion pairs
  • Validating performance by calculating quasistatic noise susceptibility

Ultimately we will demonstrate how to parallelize gates in order to speed-up quantum circuit execution in trapped-ion quantum computers, while adding robustness against common hardware error sources.

For further discussion see Numeric optimization for configurable, parallel, error-robust entangling gates in large ion registers, published as Advanced Quantum Technology 3, 11 (2020).

Worked example: Creating Q-CTRL pulses for parallel 2-of-5 and 3-of-5 qubit gates

Using the Q-CTRL Python package, you can obtain control solutions for configurable multi-qubit gates. These gates can also be performed simultaneously (in parallel). In this example, you'll learn to perform simultaneous two-qubit and three-qubit entangling operations.

First specify the ion, trap and laser characteristics.

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

from qctrl import Qctrl

plt.style.use(get_qctrl_style())


# Start a session with the API
qctrl = Qctrl()

# Trap characteristics
number_of_ions = 5

# Laser characteristics
maximum_rabi_rate = 2 * np.pi * 100e3
laser_detuning = 1.6e6 + 4.7e3

# Derived setup parameters
ion_chain_properties = qctrl.functions.calculate_ion_chain_properties(
    atomic_mass=171,  # Yb ions
    ion_count=number_of_ions,
    # Center of mass frequencies
    radial_x_center_of_mass_frequency=1.6e6,
    radial_y_center_of_mass_frequency=1.5e6,
    axial_center_of_mass_frequency=0.3e6,
    # Laser wavevector
    radial_x_wave_number=(2 * np.pi) / 355e-9,
    radial_y_wave_number=(2 * np.pi) / 355e-9,
    axial_wave_number=0,
)
Your task calculate_ion_chain_properties has completed.

We can then collect the Lamb–Dicke parameters and relative detunings from the result, which will be used in the phase and displacement calculations for ions. For more information, see the How to optimize error-robust Mølmer–Sørensen gates for trapped ions user guide.

# Collect Lamb Dicke parameters in the shape
# [<axis>, <collective_mode>, <ion>]
lamb_dicke_parameters = np.asarray(
    [
        [
            mode.lamb_dicke_parameters
            for mode in ion_chain_properties.radial_x_mode_properties
        ],
        [
            mode.lamb_dicke_parameters
            for mode in ion_chain_properties.radial_y_mode_properties
        ],
        [
            mode.lamb_dicke_parameters
            for mode in ion_chain_properties.axial_mode_properties
        ],
    ]
)

# Collect relative detunings in the shape
# [<axis>, <collective_mode>]
relative_detunings = (
    np.asarray(
        [
            [mode.frequency for mode in ion_chain_properties.radial_x_mode_properties],
            [mode.frequency for mode in ion_chain_properties.radial_y_mode_properties],
            [mode.frequency for mode in ion_chain_properties.axial_mode_properties],
        ]
    )
    - laser_detuning
)

To demonstrate simultaneous two- and three-qubit gates, configure the target operation by specifying target relative phases between each ion pair. Element ($j$, $k$) of the target phase array is the relative phase between ions $j$ and $k$, with $k<j$. In the cell below, the target relative phase for the ion pair (1, 0) is set to $\pi/4$ for maximal two-qubit entanglement. For the three-qubit gate, the same relative phase is specified for ion pairs (3, 2), (4, 2) and (4, 3).

# Operation duration in seconds
duration = 3e-4

# Target phases: element (j,k) is the relative entanglement for ions j and k (k<j)
target = np.array(
    [
        [0, 0, 0, 0, 0],
        [np.pi / 4, 0, 0, 0, 0],
        [0, 0, 0, 0, 0],
        [0, 0, np.pi / 4, 0, 0],
        [0, 0, np.pi / 4, np.pi / 4, 0],
    ]
)

Optimize drives (optimal control solution)

To control the system, consider separately-tunable complex drives $\gamma_j (t) = \Omega_j e^{i \phi_j}$ that address each ion $j$ in the trap. The drives are piecewise-constant with amplitude and phase modulation. In the optimization, the gate infidelity quantifies the solution performance.

number_of_segments = 128
sample_times = np.linspace(0, duration, number_of_segments)

drive_names = ["ion_drive" + str(number) for number in range(number_of_ions)]

# Helper function for robustness
def reflect_signal(number_of_segments, phases, moduli):
    """This function reflects a drive signal about its temporal midpoint (Milne et al., Phys. Rev. Applied, 2020).
    The combined signal is returned."""

    phases_diff = phases[1:] - phases[:-1]
    if number_of_segments % 2 == 0:
        moduli_refl = qctrl.operations.reverse(moduli, [0])
        phases_diff_refl = qctrl.operations.concatenate(
            [np.array([0]), qctrl.operations.reverse(phases_diff, [0])],
            0,
        )
        phases_extra = phases[-1] + qctrl.operations.cumulative_sum(phases_diff_refl)
    else:
        moduli_refl = qctrl.operations.reverse(moduli[:-1], [0])
        phases_diff_refl = qctrl.operations.reverse(phases_diff, [0])
        phases_extra = phases[-1] + qctrl.operations.cumulative_sum(phases_diff_refl)

    moduli_comb = qctrl.operations.concatenate([moduli, moduli_refl], 0)
    phases_comb = qctrl.operations.concatenate([phases, phases_extra], 0)
    return moduli_comb, phases_comb


# Helper function for optimization with different drives for each ion
def optimization_with_different_drives(
    ion_count,
    number_of_segments,
    duration,
    maximum_rabi_rate,
    lamb_dicke_parameters,
    relative_detunings,
    sample_times,
    robust,
    drive_names,
):
    with qctrl.create_graph() as graph:
        # Specification of free variables and combination with reflected signal
        number_of_free_segments = number_of_segments
        if robust == True:
            number_of_free_segments = int(np.ceil(number_of_segments / 2))

        drives = []
        for drive_name in drive_names:
            # The drive moduli are free variables here. They can also be restricted or fixed.
            drive_moduli = qctrl.operations.optimization_variable(
                count=number_of_free_segments,
                lower_bound=0,
                upper_bound=maximum_rabi_rate,
            )
            # The drive phases are free variables here. They can also be restricted or fixed.
            drive_phases = qctrl.operations.optimization_variable(
                count=number_of_free_segments,
                lower_bound=0,
                upper_bound=2 * np.pi,
                is_lower_unbounded=True,
                is_upper_unbounded=True,
            )

            full_moduli = drive_moduli
            full_phases = drive_phases

            if robust == True:
                full_moduli, full_phases = reflect_signal(
                    number_of_segments, drive_phases, drive_moduli
                )

            drives.append(
                qctrl.operations.complex_pwc_signal(
                    moduli=full_moduli,
                    phases=full_phases,
                    duration=duration,
                    name=drive_name,
                )
            )

        ms_phases = qctrl.operations.ms_phases(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings,
            sample_times=sample_times,
            name="phases",
        )

        ms_displacements = qctrl.operations.ms_displacements(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings,
            sample_times=sample_times,
            name="displacements",
        )

        infidelity = qctrl.operations.ms_infidelity(
            phases=ms_phases[-1],
            displacements=ms_displacements[-1],
            target_phases=target,
            name="infidelity",
        )

        if not robust:
            cost = infidelity + 0
        else:
            robust_cost_term = qctrl.operations.ms_dephasing_robust_cost(
                drives=drives,
                lamb_dicke_parameters=lamb_dicke_parameters,
                relative_detunings=relative_detunings,
            )
            cost = infidelity + robust_cost_term

        cost.name = "cost"

    return qctrl.functions.calculate_optimization(
        graph=graph,
        optimization_count=10,
        cost_node_name="cost",
        output_node_names=["displacements", "infidelity", "phases"] + drive_names,
    )
result_optimal = optimization_with_different_drives(
    ion_count=number_of_ions,
    number_of_segments=number_of_segments,
    duration=duration,
    maximum_rabi_rate=maximum_rabi_rate,
    lamb_dicke_parameters=lamb_dicke_parameters,
    relative_detunings=relative_detunings,
    sample_times=sample_times,
    robust=False,
    drive_names=drive_names,
)
Your task calculate_optimization has started.
Your task calculate_optimization has completed.
print("Cost = Infidelity = ", result_optimal.output["infidelity"]["value"])
Cost = Infidelity =  3.6536125236352746e-09
control = {f"$\gamma$": result_optimal.output[drive_names[0]]}
plot_controls(plt.figure(), control)

The figure displays the optimized pulse modulus and phase dynamics for ion 0.

Optimize drives (robust control solution)

To obtain robust controls, we impose symmetry on each drive and optimize such that the centre of mass of each mode's trajectory is at zero. This is the same procedure as for the two-qubit case described here. In this case, the solution performance is given by the sum of the gate infidelity and the robustness cost. The robust option is enabled by setting the flag robust to True.

# Robust optimization
number_of_optimizations = 10
result_robust = optimization_with_different_drives(
    ion_count=number_of_ions,
    number_of_segments=number_of_segments,
    duration=duration,
    maximum_rabi_rate=maximum_rabi_rate,
    lamb_dicke_parameters=lamb_dicke_parameters,
    relative_detunings=relative_detunings,
    sample_times=sample_times,
    robust=True,
    drive_names=drive_names,
)
Your task calculate_optimization has started.
Your task calculate_optimization has completed.
print("Cost = ", result_robust.cost)
print("Infidelity = ", result_robust.output["infidelity"]["value"])
Cost =  8.491882491566308e-08
Infidelity =  1.7057765089312227e-08
control = {f"$\gamma$": result_robust.output[drive_names[0]]}
plot_controls(plt.figure(), control)

The above figure displays the dynamics of the modulus and angle of the robust optimized pulse for ion 0. The symmetry in time of the modulus values can be observed.

Performance validation

Calculate phase space trajectories

Next, we visualize the trajectory of the center of a coherent state in (rotating) optical phase space for each mode. Note that there are $3N$ modes in general, where $N$ is the ion number. We address only the $2N$ transverse modes in the trap. The closure of these trajectories is a necessary condition for an effective operation. We first examine the (non-robust) standard optimized control, followed by the robust optimized control.

# Sum over ion index to obtain total displacement of the mode.
trajl = np.sum(result_optimal.output["displacements"]["value"], axis=3)
scale = 1.1 * np.max(np.abs(trajl))

gs = gridspec.GridSpec(1, 2, wspace=0.1)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(2.0 * 5)
fig.suptitle("Phase space trajectories", fontsize=16, y=1.0)

ax = fig.add_subplot(gs[0])
ax.tick_params(labelsize=14)
for mode in range(number_of_ions):
    ax.plot(
        np.real(trajl[:, 0, mode]),
        np.imag(trajl[:, 0, mode]),
        label="mode " + str(mode % 5),
        linewidth=2,
    )
    ax.plot(
        np.real(trajl[-1, 0, mode]),
        np.imag(trajl[-1, 0, mode]),
        marker="x",
        color="black",
        markersize=15,
    )
ax.set_xlim(-scale, scale)
ax.set_ylim(-scale, scale)
ax.set_aspect("equal")
ax.set_xlabel("q", fontsize=14)
ax.set_ylabel("p", fontsize=14)
ax.set_title("$x$-axis", fontsize=14)

ax = fig.add_subplot(gs[1])
ax.tick_params(labelsize=14)
for mode in range(number_of_ions):
    ax.plot(
        np.real(trajl[:, 1, mode]),
        np.imag(trajl[:, 1, mode]),
        label="mode " + str(mode % 5),
        linewidth=2,
    )
    ax.plot(
        np.real(trajl[-1, 1, mode]),
        np.imag(trajl[-1, 1, mode]),
        marker="x",
        color="black",
        markersize=15,
    )
ax.set_xlim(-scale, scale)
ax.set_ylim(-scale, scale)
ax.set_aspect("equal")
ax.set_xlabel("q", fontsize=14)
ax.yaxis.set_ticklabels([])
ax.set_title("$y$-axis", fontsize=14)
ax.legend(loc="best", bbox_to_anchor=(1.45, 0.8), fontsize=14)

plt.show()

The black cross marks the final displacement for each mode. These are overlapping at zero, indicating no residual state-motional entanglement and no motional heating caused by the operations.

Now we visualize the phase space trajectories for the robust optimized control.

# Sum over ion index to obtain total displacement of the mode.
trajl_r = np.sum(result_robust.output["displacements"]["value"], axis=3)
scale = 1.1 * np.max(np.abs(trajl_r))

gs = gridspec.GridSpec(1, 2, wspace=0.1)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(2.0 * 5)
fig.suptitle("Phase space trajectories", fontsize=16, y=1.0)

ax = fig.add_subplot(gs[0])
ax.tick_params(labelsize=14)
for mode in range(number_of_ions):
    ax.plot(
        np.real(trajl_r[:, 0, mode]),
        np.imag(trajl_r[:, 0, mode]),
        label="mode " + str(mode % 5),
        linewidth=2,
    )
    ax.plot(
        np.real(trajl_r[-1, 0, mode]),
        np.imag(trajl_r[-1, 0, mode]),
        marker="x",
        color="black",
        markersize=15,
    )
ax.set_xlim(-scale, scale)
ax.set_ylim(-scale, scale)
ax.set_aspect("equal")
ax.set_xlabel("q", fontsize=14)
ax.set_ylabel("p", fontsize=14)
ax.set_title("$x$-axis", fontsize=14)

ax = fig.add_subplot(gs[1])
ax.tick_params(labelsize=14)
for mode in range(number_of_ions):
    ax.plot(
        np.real(trajl_r[:, 1, mode]),
        np.imag(trajl_r[:, 1, mode]),
        label="mode " + str(mode % 5),
        linewidth=2,
    )
    ax.plot(
        np.real(trajl_r[-1, 1, mode]),
        np.imag(trajl_r[-1, 1, mode]),
        marker="x",
        color="black",
        markersize=15,
    )
ax.set_xlim(-scale, scale)
ax.set_ylim(-scale, scale)
ax.set_aspect("equal")
ax.set_xlabel("q", fontsize=14)
ax.yaxis.set_ticklabels([])
ax.set_title("$y$-axis", fontsize=14)
ax.legend(loc="best", bbox_to_anchor=(1.45, 0.8), fontsize=14)

plt.show()

Again, the black crosses at the origin indicate no residual state-motional entanglement, which arises from satisfying the center of mass and symmetry conditions.

Calculate relative phase dynamics

We can also obtain the phase accumulation for each pair of ions. The target phases for each pair of ions should be achieved at the end of a successful operation. First consider the standard optimized control.

phases = result_optimal.output["phases"]["value"]
gs = gridspec.GridSpec(2, 1)
fig = plt.figure()
fig.set_figheight(1.8 * 5)
fig.set_figwidth(10)
fig.suptitle("Relative phase dynamics", fontsize=16, y=0.95)

ax0 = fig.add_subplot(gs[0])
ax0.tick_params(labelsize=14)
ax1 = fig.add_subplot(gs[1])
ax1.tick_params(labelsize=14)
x_timing = sample_times * 1e3

for ion1 in range(number_of_ions):
    for ion2 in range(ion1):
        if target[ion1][ion2] == np.pi / 4:
            ax0.plot(
                x_timing,
                phases[:, ion1, ion2],
                label="(" + str(ion1) + ", " + str(ion2) + ")",
                linewidth=2,
            )
        else:
            ax1.plot(
                x_timing,
                phases[:, ion1, ion2],
                label="(" + str(ion1) + ", " + str(ion2) + ")",
                linewidth=2,
            )

unit = np.pi / 4
ymax = np.pi / 4
ymin = 0 * np.pi / 4
y_tick = np.arange(ymin, ymax + unit, unit)
y_label = [r"0", r"$\pi/4$"]
ax0.set_yticks(y_tick)
ax0.set_yticklabels(y_label)
ax1.set_yticks(y_tick)
ax1.set_yticklabels(y_label)
ax0.plot([0, x_timing[-1]], [np.pi / 4, np.pi / 4], "k--", linewidth=0.8)
ax1.plot([0, x_timing[-1]], [0, 0], "k--", linewidth=0.8)
ax0.plot([x_timing[-1], x_timing[-1]], [-4, 4], "k", alpha=0.5)
ax1.plot([x_timing[-1], x_timing[-1]], [-4, 4], "k", alpha=0.5)
ax0.set_ylim((-0.1, 0.9))
ax1.set_ylim((-0.1, 0.9))
ax0.xaxis.set_ticklabels([])
ax1.set_xlabel("Time (ms)", fontsize=14)
ax0.set_ylabel("Relative phase", fontsize=14)
ax1.set_ylabel("Relative phase", fontsize=14)
ax0.set_title("Pairs with target phase $\psi = \pi/4$", fontsize=14)
ax1.set_title("Pairs with target phase $\psi = 0$", fontsize=14)
ax0.legend(
    bbox_to_anchor=(1.02, 0.82),
    loc="upper left",
    borderaxespad=0.0,
    fontsize=14,
    title="Ion pair",
    title_fontsize=14,
)
ax1.legend(
    bbox_to_anchor=(1.02, 0.82),
    loc="upper left",
    borderaxespad=0.0,
    fontsize=14,
    title="Ion pair",
    title_fontsize=14,
)
plt.show()

These figures display the relative phase dynamics for the duration of the gate operation, and the end of the operation is marked by vertical black lines. For clarity, different plots display the ion pairs with different target relative phases. The two plots are for target relative phases of $\pi/4$ and 0, and these values are marked with black horizontal dashed lines.

Next consider the phase accumulation for the robust optimized control.

phases_r = result_robust.output["phases"]["value"]
gs = gridspec.GridSpec(2, 1)
fig = plt.figure()
fig.set_figheight(1.8 * 5)
fig.set_figwidth(10)
fig.suptitle("Relative phase dynamics", fontsize=16, y=0.95)

ax0 = fig.add_subplot(gs[0])
ax0.tick_params(labelsize=14)
ax1 = fig.add_subplot(gs[1])
ax1.tick_params(labelsize=14)
x_timing = sample_times * 1e3

for ion1 in range(number_of_ions):
    for ion2 in range(ion1):
        if target[ion1][ion2] == np.pi / 4:
            ax0.plot(
                x_timing,
                phases_r[:, ion1, ion2],
                label="(" + str(ion1) + ", " + str(ion2) + ")",
                linewidth=2,
            )
        else:
            ax1.plot(
                x_timing,
                phases_r[:, ion1, ion2],
                label="(" + str(ion1) + ", " + str(ion2) + ")",
                linewidth=2,
            )

unit = np.pi / 4
ymax = np.pi / 4
ymin = 0 * np.pi / 4
y_tick = np.arange(ymin, ymax + unit, unit)
y_label = [r"0", r"$\pi/4$"]
ax0.set_yticks(y_tick)
ax0.set_yticklabels(y_label)
ax1.set_yticks(y_tick)
ax1.set_yticklabels(y_label)
ax0.plot([0, x_timing[-1]], [np.pi / 4, np.pi / 4], "k--", linewidth=0.8)
ax1.plot([0, x_timing[-1]], [0, 0], "k--", linewidth=0.8)
ax0.plot([x_timing[-1], x_timing[-1]], [-4, 4], "k", alpha=0.5)
ax1.plot([x_timing[-1], x_timing[-1]], [-4, 4], "k", alpha=0.5)
ax0.set_ylim((-0.1, 0.9))
ax1.set_ylim((-0.1, 0.9))
ax0.xaxis.set_ticklabels([])
ax1.set_xlabel("Time (ms)", fontsize=14)
ax0.set_ylabel("Relative phase", fontsize=14)
ax1.set_ylabel("Relative phase", fontsize=14)
ax0.set_title("Pairs with target phase $\psi = \pi/4$", fontsize=14)
ax1.set_title("Pairs with target phase $\psi = 0$", fontsize=14)
ax0.legend(
    bbox_to_anchor=(1.02, 0.9),
    loc="upper left",
    borderaxespad=0.0,
    fontsize=14,
    title="Ion pair",
    title_fontsize=14,
)
ax1.legend(
    bbox_to_anchor=(1.02, 0.9),
    loc="upper left",
    borderaxespad=0.0,
    fontsize=14,
    title="Ion pair",
    title_fontsize=14,
)
plt.show()

As above, these figures display the relative phase dynamics for the duration of the gate operation, and the end of the operation is marked by vertical black lines. The two plots are for target relative phases of $\pi/4$ and 0, and these values are marked with black horizontal dashed lines.

The optimized drives produce nontrivial relative phase dynamics for each ion pair since the ions are individually addressed. The entangling phase is mediated by coupling the ion pairs to shared motional modes, and the optimized drives (both standard and robust) provide the necessary degrees of freedom to achieve the different target relative phases for each ion pair.

Calculate robustness to dephasing and pulse timing errors with quasi-static scans

We can assess the robustness of the 'Robust' optimized pulses in comparison to the 'Standard' optimized pulses using 1D quasi-static scans.

First, we calculate a scan of scaling factors for the pulse timings. The scaling factors are applied uniformly: they scale the timing for the entire operation by the same factor.

maxshift = 0.004

timing_scan_points = 21

shifts = np.linspace(-maxshift, maxshift, timing_scan_points)
optimal_infidelity_names = [
    "infidelity" + str(number) for number in range(timing_scan_points)
]

with qctrl.create_graph() as graph:
    for (
        shift,
        infidelity_name,
    ) in zip(shifts, optimal_infidelity_names):
        drives = [
            qctrl.operations.pwc_signal(
                values=np.array(
                    [segment["value"] for segment in result_optimal.output[drive_name]]
                ),
                duration=duration * (1 + shift),
            )
            for drive_name in drive_names
        ]

        ms_phases = qctrl.operations.ms_phases(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings,
        )

        ms_displacements = qctrl.operations.ms_displacements(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings,
        )

        ms_infidelity = qctrl.operations.ms_infidelity(
            phases=ms_phases,
            displacements=ms_displacements,
            target_phases=target,
            name=infidelity_name,
        )

timing_scan = qctrl.functions.calculate_graph(
    graph=graph,
    output_node_names=optimal_infidelity_names,
)

infids_timing = [timing_scan.output[name]["value"] for name in optimal_infidelity_names]
Your task calculate_graph has started.
Your task calculate_graph has completed.
robust_infidelity_names = [
    "robust_infidelity" + str(number) for number in range(timing_scan_points)
]

with qctrl.create_graph() as graph:
    for shift, robust_infidelity_name in zip(shifts, robust_infidelity_names):
        drives = [
            qctrl.operations.pwc_signal(
                values=np.array(
                    [segment["value"] for segment in result_robust.output[drive_name]]
                ),
                duration=duration * (1 + shift),
            )
            for drive_name in drive_names
        ]

        ms_phases = qctrl.operations.ms_phases(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings,
        )

        ms_displacements = qctrl.operations.ms_displacements(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings,
        )

        ms_infidelity = qctrl.operations.ms_infidelity(
            phases=ms_phases,
            displacements=ms_displacements,
            target_phases=target,
            name=robust_infidelity_name,
        )

timing_scan = qctrl.functions.calculate_graph(
    graph=graph,
    output_node_names=robust_infidelity_names,
)

robust_infids_timing = [
    timing_scan.output[name]["value"] for name in robust_infidelity_names
]
Your task calculate_graph has started.
Your task calculate_graph has completed.

Next, calculate the robustness of the optimized pulses using a systematic scan of the relative detunings (which corresponds to shifting the laser detuning).

minreldet = np.min(np.abs(relative_detunings))
maxshift = 4 * minreldet / 10000 * 100

dephasing_scan_points = 31

shifts = np.linspace(-maxshift, maxshift, dephasing_scan_points)
optimal_infidelity_names = [
    "infidelity" + str(number) for number in range(dephasing_scan_points)
]

with qctrl.create_graph() as graph:
    for shift, infidelity_name in zip(shifts, optimal_infidelity_names):
        drives = [
            qctrl.operations.pwc_signal(
                values=np.array(
                    [segment["value"] for segment in result_optimal.output[drive_name]]
                ),
                duration=duration,
            )
            for drive_name in drive_names
        ]

        ms_phases = qctrl.operations.ms_phases(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings + shift,
        )

        ms_displacements = qctrl.operations.ms_displacements(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings + shift,
        )

        ms_infidelity = qctrl.operations.ms_infidelity(
            phases=ms_phases,
            displacements=ms_displacements,
            target_phases=target,
            name=infidelity_name,
        )

dephasing_scan = qctrl.functions.calculate_graph(
    graph=graph,
    output_node_names=optimal_infidelity_names,
)

infids_dephasing = [
    dephasing_scan.output[name]["value"] for name in optimal_infidelity_names
]
Your task calculate_graph has started.
Your task calculate_graph has completed.
robust_infidelity_names = [
    "robust_infidelity" + str(number) for number in range(dephasing_scan_points)
]

with qctrl.create_graph() as graph:
    for shift, robust_infidelity_name in zip(shifts, robust_infidelity_names):
        drives = [
            qctrl.operations.pwc_signal(
                values=np.array(
                    [segment["value"] for segment in result_robust.output[drive_name]]
                ),
                duration=duration,
            )
            for drive_name in drive_names
        ]

        ms_phases = qctrl.operations.ms_phases(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings + shift,
        )
        ms_displacements = qctrl.operations.ms_displacements(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings + shift,
        )
        ms_infidelity = qctrl.operations.ms_infidelity(
            phases=ms_phases,
            displacements=ms_displacements,
            target_phases=target,
            name=robust_infidelity_name,
        )

dephasing_scan = qctrl.functions.calculate_graph(
    graph=graph,
    output_node_names=robust_infidelity_names,
)

robust_infids_dephasing = [
    dephasing_scan.output[name]["value"] for name in robust_infidelity_names
]
Your task calculate_graph has started.
Your task calculate_graph has completed.
gs = gridspec.GridSpec(1, 2)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(3 * 5)
fig.suptitle("Quasi-static scans", fontsize=16, y=1.15)

ax = fig.add_subplot(gs[0])
ax.set_title("Timing noise", fontsize=14)
maxshift = 0.004
timeshifts = np.linspace(1 - maxshift, 1 + maxshift, timing_scan_points)
ax.plot(timeshifts, infids_timing, label="Standard", color="k", linewidth=2)
ax.plot(timeshifts, robust_infids_timing, label="Robust", linewidth=2)
ax.tick_params(labelsize=14)
ax.set_xlabel("Pulse timing scaling factor", fontsize=14)
ax.set_ylabel("Infidelity", fontsize=14)

ax = fig.add_subplot(gs[1])
ax.set_title("Dephasing noise", fontsize=14)
maxshift = 4 * minreldet / 10000 * 100
detshifts = np.linspace(-maxshift, maxshift, dephasing_scan_points)
ax.plot(detshifts / 1000, infids_dephasing, label="Standard", color="k", linewidth=2)
ax.plot(
    detshifts / 1000,
    robust_infids_dephasing,
    label="Robust",
    linewidth=2,
)
ax.tick_params(labelsize=14)
ax.set_xlabel("Laser detuning shift $\Delta \delta$ (kHz)", fontsize=14)
ax.set_ylabel("Infidelity", fontsize=14)
ax.legend(loc="best", bbox_to_anchor=(0.06, 1.28), fontsize=14)
plt.show()

The broader high-fidelity regions indicate the benefit of the robust optimized pulses when there is quasi-static dephasing noise or noise in the control pulse timings. The additional robustness for this more general gate is in agreement with published experimental results, using related robustness methodology.

Worked example: Creating smoothed parallel gates with user-defined entangling phases

In this example we highlight the flexibility of our graph-based computational framework to add several substantial constraints on top of our optimization.

  • We generate six parallel pairwise "2-of-12" entangling operations on a 12-ion chain
  • Each pair is assigned a different target entangling phase, increasing in steps of 0.2 for each pair
  • The high-frequency components of the drive are removed using a sinc filter within the optimizer
# Trap characteristics
number_of_ions = 12

# Laser characteristics
maximum_rabi_rate = 2 * np.pi * 100e3
laser_detuning = 1.6e6 + 4.7e3

# Derived setup parameters
ion_chain_properties = qctrl.functions.calculate_ion_chain_properties(
    atomic_mass=171,  # Yb ions
    ion_count=number_of_ions,
    # Center of mass frequencies
    radial_x_center_of_mass_frequency=1.6e6,
    radial_y_center_of_mass_frequency=1.5e6,
    axial_center_of_mass_frequency=0.25e6,
    # Laser wavevector
    radial_x_wave_number=(2 * np.pi) / 355e-9,
    radial_y_wave_number=(2 * np.pi) / 355e-9,
    axial_wave_number=0,
)
Your task calculate_ion_chain_properties has completed.
# Collect Lamb Dicke parameters in the shape
# [<axis>, <collective_mode>, <ion>]
lamb_dicke_parameters = np.asarray(
    [
        [
            mode.lamb_dicke_parameters
            for mode in ion_chain_properties.radial_x_mode_properties
        ],
        [
            mode.lamb_dicke_parameters
            for mode in ion_chain_properties.radial_y_mode_properties
        ],
        [
            mode.lamb_dicke_parameters
            for mode in ion_chain_properties.axial_mode_properties
        ],
    ]
)

# Collect relative detunings in the shape
# [<axis>, <collective_mode>]
relative_detunings = (
    np.asarray(
        [
            [mode.frequency for mode in ion_chain_properties.radial_x_mode_properties],
            [mode.frequency for mode in ion_chain_properties.radial_y_mode_properties],
            [mode.frequency for mode in ion_chain_properties.axial_mode_properties],
        ]
    )
    - laser_detuning
)
# Operation duration in seconds
duration = 3e-4

# Define target phases: element (j,k) is the relative entanglement for ions j and k (k<j)
target = np.zeros((12, 12))
for ion1 in range(1, 12, 2):
    ion2 = ion1 - 1
    target[ion1][ion2] = ion1 / 10

Generating optimized pulses

Consider separately-tunable complex drives $\gamma_j (t) = \Omega_j e^{i \phi_j}$ that address each ion $j$ in the trap. The drives are piecewise-constant with amplitude and phase modulation. Additionally, a low-pass (sinc) filter is incorporated into the optimizer to smooth the pulse, as required for many practical implementations.

You can set the number of optimized pulse segments and resampling segments for the smoothed pulse, as well as the sinc cutoff frequency (in Hz) in the cell below.

number_of_segments = 32
number_of_sample_segments = 128
sample_times = np.linspace(0, duration, number_of_sample_segments)
sinc_cutoff_frequency_drives = 2 * np.pi * 0.05e6
drive_names = ["ion_drive" + str(number) for number in range(number_of_ions)]

with qctrl.create_graph() as graph:
    sinc_kernel = qctrl.operations.sinc_convolution_kernel(sinc_cutoff_frequency_drives)

    drives = []

    for drive_name in drive_names:
        drive_real_raw_values = qctrl.operations.optimization_variable(
            count=number_of_segments,
            lower_bound=-maximum_rabi_rate,
            upper_bound=maximum_rabi_rate,
        )
        drive_imag_raw_values = qctrl.operations.optimization_variable(
            count=number_of_segments,
            lower_bound=-maximum_rabi_rate,
            upper_bound=maximum_rabi_rate,
        )

        drive_raw = qctrl.operations.pwc_signal(
            values=drive_real_raw_values + 1j * drive_imag_raw_values,
            duration=duration,
        )

        drive_filtered = qctrl.operations.convolve_pwc(
            pwc=drive_raw, kernel=sinc_kernel
        )

        drives.append(
            qctrl.operations.discretize_stf(
                stf=drive_filtered,
                duration=duration,
                segment_count=number_of_sample_segments,
                name=drive_name,
            )
        )

    ms_phases = qctrl.operations.ms_phases(
        drives=drives,
        lamb_dicke_parameters=lamb_dicke_parameters,
        relative_detunings=relative_detunings,
    )

    ms_displacements = qctrl.operations.ms_displacements(
        drives=drives,
        lamb_dicke_parameters=lamb_dicke_parameters,
        relative_detunings=relative_detunings,
    )

    infidelity = qctrl.operations.ms_infidelity(
        phases=ms_phases,
        displacements=ms_displacements,
        target_phases=target,
        name="infidelity",
    )
number_of_optimizations = 3
result = qctrl.functions.calculate_optimization(
    graph=graph,
    optimization_count=number_of_optimizations,
    cost_node_name="infidelity",
    output_node_names=["infidelity"] + drive_names,
)
Your task calculate_optimization has started.
Your task calculate_optimization has completed.
print("Cost = Infidelity = ", result.output["infidelity"]["value"])
Cost = Infidelity =  4.00294495150888e-06
control = {f"$\gamma$": result.output["ion_drive0"]}
plot_controls(plt.figure(), control)

The above figure displays the optimized smooth pulse modulus and phase dynamics for ion 0, as an example of the ion-specific drives.

Performance validation

Calculate phase space trajectories

For a system as large as this, you can separate the simulation of the displacement from the optimization to increase the efficiency. This procedure is demonstrated in the How to calculate phase and motion dynamics for arbitrarily modulated Mølmer–Sørensen gates User guide.

with qctrl.create_graph() as graph:
    drives = [
        qctrl.operations.pwc_signal(
            values=np.array(
                [segment["value"] for segment in result.output[drive_name]]
            ),
            duration=duration,
        )
        for drive_name in drive_names
    ]

    ms_displacements = qctrl.operations.ms_displacements(
        drives=drives,
        lamb_dicke_parameters=lamb_dicke_parameters,
        relative_detunings=relative_detunings,
        sample_times=sample_times,
        name="displacements",
    )

displacement_simulation = qctrl.functions.calculate_graph(
    graph=graph,
    output_node_names=["displacements"],
)

# Sum over ion index to obtain total displacement of the mode.
trajl = np.sum(displacement_simulation.output["displacements"]["value"], axis=3)
Your task calculate_graph has completed.
scale = 1.1 * np.max(np.abs(trajl))

gs = gridspec.GridSpec(1, 2, wspace=0.1)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(2.0 * 5)
fig.suptitle("Phase space trajectories", fontsize=16, y=1.0)

ax = fig.add_subplot(gs[0])
ax.tick_params(labelsize=14)
for mode in range(number_of_ions):
    ax.plot(
        np.real(trajl[:, 0, mode]),
        np.imag(trajl[:, 0, mode]),
        label="mode " + str(mode % 12),
        linewidth=2,
    )
    ax.plot(
        np.real(trajl[-1, 0, mode]),
        np.imag(trajl[-1, 0, mode]),
        marker="x",
        color="black",
        markersize=15,
    )
ax.set_xlim(-scale, scale)
ax.set_ylim(-scale, scale)
ax.set_aspect("equal")
ax.set_xlabel("q", fontsize=14)
ax.set_ylabel("p", fontsize=14)
ax.set_title("$x$-axis", fontsize=14)

ax = fig.add_subplot(gs[1])
ax.tick_params(labelsize=14)
for mode in range(number_of_ions):
    ax.plot(
        np.real(trajl[:, 1, mode]),
        np.imag(trajl[:, 1, mode]),
        label="mode " + str(mode % 12),
        linewidth=2,
    )
    ax.plot(
        np.real(trajl[-1, 1, mode]),
        np.imag(trajl[-1, 1, mode]),
        marker="x",
        color="black",
        markersize=15,
    )
ax.set_xlim(-scale, scale)
ax.set_ylim(-scale, scale)
ax.set_aspect("equal")
ax.set_xlabel("q", fontsize=14)
ax.yaxis.set_ticklabels([])
ax.set_title("$y$-axis", fontsize=14)
ax.legend(loc="best", bbox_to_anchor=(1.5, 1.0), fontsize=14)

plt.show()

The black cross marks the final displacement for each mode. These are overlapping at zero, indicating no residual state-motional entanglement and no motional heating caused by the operations.

Calculate relative phase dynamics

The target phases for each pair of ions should be achieved at the end of a successful operation. Again, for a system as large as this, you can separate the simulation of the phases from the optimization, to increase the efficiency.

with qctrl.create_graph() as graph:
    drives = [
        qctrl.operations.pwc_signal(
            values=np.array(
                [segment["value"] for segment in result.output[drive_name]]
            ),
            duration=duration,
        )
        for drive_name in drive_names
    ]

    ms_phases = qctrl.operations.ms_phases(
        drives=drives,
        lamb_dicke_parameters=lamb_dicke_parameters,
        relative_detunings=relative_detunings,
        sample_times=sample_times,
        name="phases",
    )

phase_simulation = qctrl.functions.calculate_graph(
    graph=graph,
    output_node_names=["phases"],
)

phases_6 = phase_simulation.output["phases"]["value"]
Your task calculate_graph has started.
Your task calculate_graph has completed.
colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]

fig, ax = plt.subplots()
fig.set_figheight(5)
fig.set_figwidth(10)
plt.title("Relative phase dynamics", fontsize=16, y=1)
ax.tick_params(labelsize=14)

x_timing = sample_times * 1e3
ind = 0
for ion1 in range(1, 12, 2):
    ion2 = ion1 - 1
    ax.plot(
        x_timing,
        np.array(phases_6[:, ion1, ion2]),
        label="(" + str(ion1) + ", " + str(ion2) + ")",
        linewidth=2,
        color=colors[ind],
    )
    ax.plot(
        [0, x_timing[-1]],
        [ion1 / 10, ion1 / 10],
        "--",
        linewidth=0.8,
        color=colors[ind],
    )
    ind += 1

ax.plot([x_timing[-1], x_timing[-1]], [-4, 4], "k", alpha=0.5)
ax.set_ylim((-0.1, 1.2))
ax.set_xlabel("Time (ms)", fontsize=14)
ax.set_ylabel("Relative phase", fontsize=14)
ax.legend(
    bbox_to_anchor=(1.02, 0.9),
    loc="upper left",
    borderaxespad=0.0,
    fontsize=14,
    title="Ion pair",
    title_fontsize=14,
)
plt.show()

The figure displays the relative phase dynamics for each entangled ion pair, with a color matched to the pair's target relative phase (dashed). Note that by the end of the operation, marked by the black vertical line, each ion pair achieves its specified, distinct relative phase target.