Designing robust Mølmer–Sørensen gates with parametric trap drive amplification

Obtaining control solutions for two-qubit gates with modulation of the confining potential

Boulder Opal supports the development of high-performance entangling gates to enhance the quantum volume of trapped-ion quantum computers. Typically, Mølmer–Sørensen-type operations can be used to entangle multiple ions using pairs of lasers with opposite detunings and a fixed phase relationship. This operation can flexibly incorporate multiple entangling operations performed in parallel, or multi-qubit entangling operations using robust controls designed using Boulder Opal. In the Designing robust, configurable, parallel gates for large trapped-ion arrays Application Note we demonstrate how to combine these various capabilities for complex controls on large multi-ion arrays. Recent work by Burd et al. has demonstrated how modulation on the confining potential can be used to realize faster gates.

In this Application note, we describe how to combine robust-control optimization in Mølmer–Sørensen gates with parametric modulation of the trap drive. We will cover:

  • Simulating gate dynamics with a parametrically modulated trap drive
  • Optimizing Mølmer–Sørensen gates using one and multiple modes
  • Calculating phase-space trajectories for arbitrary controls
  • Validating performance of robust gates by calculating quasistatic noise susceptibility

Ultimately we demonstrate for the first time that it is possible to combine standard Mølmer–Sørensen-gate optimization techniques with parametric trap-drive modulation in order to achieve fast, robust gates.

Imports and initialization

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np
import copy
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()

# 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

Creating Q-CTRL pulses for two-qubit gates mediated by a single mode

For standard Mølmer–Sørensen gate-operations, the laser detuning should be close to a motional mode frequency, and the displacements in the phase-space of each mode should be small. At the end of the operation, the phase-space trajectories for each mode should return to zero (for no residual state-motional entanglement) and the phase acquired between each pair of ions should correspond to target entangling phases $\psi_{jk}$ between each pair of ions ($j$, $k$), indexing from 0:

$$ U = \exp\left(i \sum_{j=1}^{N-1} \sum_{k=0}^{j-1} \psi_{jk} \sigma_x^j \sigma_x^k \right). $$

This operation can incorporate multiple entangling operations performed in parallel, or multi-qubit entangling operations; our multi-qubit gate Application Note equips you to obtain and verify such operations.

Here, we examine how modulation on the confining potential amplifies the mode-mediated gate interaction to provide faster gates as demonstrated recently by Burd et al. The trap modulation produces the Hamiltonian:

$$H_P = \hbar g (a_p + a_p^\dagger)^2 \cos (\omega_P t + \theta)$$

where $a_p$ is the annihilation operator for motional mode $p$, $g$ is the parametric drive strength, $\omega_P$ is the modulation frequency and $\theta$ is a relative phase term (we set $\theta = 0$). When $\omega_P$ is chosen near resonance with the sideband of mode $p$, and for sufficiently small $g$, just a single mode is modified by the parametric drive. A Bogoliubov normal mode transformation can then be applied to map mode $p$ from the lab frame to a transformed frame. This transformation can recover the usual Mølmer–Sørensen-type Hamiltonian, but now written in terms of the transformed mode operators in addition to a frequency transformation and a drive transformation (when just a single mode mediates the gate).

We begin by exploring the impact of the trap modulation on the Mølmer–Sørensen gate time and robustness.

First, specify the number and type of ions, as well as other trap and laser characteristics.

# Trap characteristics
number_of_ions = 2

# Laser characteristics
maximum_rabi_rate = 2 * np.pi * 20e3

# Derived setup parameters
ion_chain_properties = qctrl.functions.calculate_ion_chain_properties(
    atomic_mass=25,  # Mg ions
    ion_count=number_of_ions,
    # Center of mass frequencies
    radial_x_center_of_mass_frequency=5.9e6,
    radial_y_center_of_mass_frequency=6.1e6,
    axial_center_of_mass_frequency=1.5e6,
    # Laser wavevector
    radial_x_wave_number=(2 * np.pi) / 355e-9,
    radial_y_wave_number=0,
    axial_wave_number=0,
)

# 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 mode frequencies in the shape
# [<axis>, <collective_mode>]
mode_frequencies = 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],
    ]
)
Your task calculate_ion_chain_properties has completed.

Next, the target two-qubit operation is configured 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 relative phase for ion pair (1, 0) is set to $\pi/4$ for maximal two-qubit entanglement.

target = np.array(
    [
        [0, 0],
        [np.pi / 4, 0],
    ]
)

Generating standard pulses

The fastest gate is performed by a single displacement loop of the mediating motional mode in phase space. Here we demonstrate this standard operation along with the faster operation produced by adding trap modulation. Note that the gate time is specified by $\pi/\eta \Omega_0$, where $\eta$ is the Lamb Dicke parameter, $\Omega_0$ is the drive strength, and the laser's relative detuning from the motional frequency is $\delta_p = 2 \eta \Omega_0$. The following parameters have been chosen to satisfy these conditions, taking into account the drive and $\delta_p$ transformations when $g>0$ (effectively enhancing the drive strength), as well as the further condition that $|g| < |\delta_p|$.

# Gate and control specifications
drive_g_list = [0, 2.9e3, 50e3]  # Hz
laser_detuning_list = [5.9e6 - 2.93e3, 5.9e6 - 5.005e3, 5.9e6 - 50.901e3]  # Hz
duration_list = [341.3e-6, 245.1e-6, 104.9e-6]  # s

Next we optimize the drive control for the given parameters. In the following optimization, the complex drive $\gamma (t) = \Omega e^{i \phi}$ jointly addresses both ions in the chain such that the gate is mediated only by the center-of-mass mode. The optimizer attempts to minimize a particular cost, which is specified below as the infidelity. Infidelity quantifies the solution performance, and is determined by the final mode displacements, and the difference between the realized and target relative phase.

Note that setting the parameters for the fastest single-loop gate will result in fixed drives close to the maximum Rabi rate to produce this solution. Valid gates can be found for slower durations by exploiting the freedom of the complex drive and modulating the drive controls over time.

# Optimization parameters
number_of_segments = 20
number_of_free_segments = number_of_segments
standard_results = []

for drive_g, laser_detuning, duration in zip(
    drive_g_list, laser_detuning_list, duration_list
):

    sample_times = np.linspace(0, duration, 50)

    # Set up trap modulation transformations
    relative_detunings = np.array(mode_frequencies - laser_detuning)
    mode_rel_det = relative_detunings[0, 1]
    mode_rel_det_driven = np.sqrt(mode_rel_det ** 2 - drive_g ** 2)
    relative_detunings_driven = copy.deepcopy(relative_detunings)
    relative_detunings_driven[0, 1] = mode_rel_det_driven
    r = 1 / 4 * np.log((mode_rel_det + drive_g) / (mode_rel_det - drive_g))

    # Optimize the controls
    with qctrl.create_graph() as graph:

        drive_moduli = qctrl.operations.optimization_variable(
            count=number_of_free_segments,
            lower_bound=0,
            upper_bound=maximum_rabi_rate,
        )
        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,
        )

        ion_drive = qctrl.operations.complex_pwc_signal(
            moduli=drive_moduli,
            phases=drive_phases,
            duration=duration,
            name="ion_drive",
        )
        # Drive transformation due to trap modulation
        ion_drive_mod = ion_drive * np.cosh(r) + qctrl.operations.conjugate(
            ion_drive
        ) * np.sinh(r)

        drives = [
            ion_drive_mod
        ] * number_of_ions  # Identical drives applied to each ion

        ms_phases = qctrl.operations.ms_phases(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings_driven,  # Relative detuning transformation
            sample_times=sample_times,
            name="phases",
        )

        ms_displacements = qctrl.operations.ms_displacements(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings_driven,  # Relative detuning transformation
            sample_times=sample_times,
            name="displacements",
        )

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

    result_optimal = qctrl.functions.calculate_optimization(
        graph=graph,
        optimization_count=10,
        cost_node_name="infidelity",
        output_node_names=["displacements", "infidelity", "ion_drive", "phases"],
    )
    standard_results.append(result_optimal)
    print("Drive strength (kHz):", drive_g / 1e3)
    print("Gate duration (\N{MICRO SIGN}s):", duration * 1e6)
    print("Infidelity:", result_optimal.cost)
Your task calculate_optimization has completed.
Drive strength (kHz): 0.0
Gate duration (µs): 341.3
Infidelity: 2.9287017255796854e-11
Your task calculate_optimization has completed.
Drive strength (kHz): 2.9
Gate duration (µs): 245.1
Infidelity: 7.34827310111541e-10
Your task calculate_optimization has completed.
Drive strength (kHz): 50.0
Gate duration (µs): 104.9
Infidelity: 2.1926549464978962e-10

We now extract the optimized controls.

control = {f"$\gamma$": standard_results[0].output["ion_drive"]}
plot_controls(plt.figure(), control)

The above figure displays an exemplary optimized pulse modulus $\Omega (t)$ and phase $\phi (t)$ that addresses both ions. This control verifies that the maximum Rabi rate is consistently applied to obtain the fastest gate duration.

We can also observe the transformed mode displacements: the trajectory of the center of a coherent state in (rotating) optical phase space for each mode.

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

colors = ["#FD8F35", "#3273DC"]
scale = 1.1 * np.max(np.abs(trajl))

fig, ax = plt.subplots(figsize=(5, 5))
fig.suptitle("Phase space trajectories", fontsize=16, y=1.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 % 2),
        color=colors[mode],
        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.legend(loc="best", fontsize=14)
plt.show()

The above figure displays an example of a phase-space trajectory, chosen for $g=50$ kHz. Here, $q \equiv q_m = (a_m^\dagger + a_m)/\sqrt{2}$ and $p \equiv p_m = i (a_m^\dagger - a_m)/\sqrt{2}$ are the dimensionless quadratures for each mode $m$. Importantly, here only the center-of-mass mode mediates the operation, and the trajectory is displayed for the transformed mode. The trajectories are circular for the fastest gate durations and any parametric drive strength.

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.

Generating robust pulses

Robustness to pulse timing errors and mode or laser frequency noise can be obtained by imposing additional conditions on the system dynamics. Robustness to these noise sources is given by a combination of symmetry in each ion's drive, as described in (Milne et al., Phys. Rev. Applied, 2020), and optimization such that the center of mass of each mode's trajectory is at zero. You can do this using the result of qctrl.operations.ms_dephasing_robust_cost as an extra term in the cost function. For more information, see the Optimizing Mølmer–Sørensen gates user guide.

It is important to note that satisfying these robustness conditions necessarily requires longer gate durations. Additionally, since the mediating motional mode is transformed according to the trap parametric drive, the center of mass condition is calculated below for the transformed mode and does not precisely represent robustness to fluctuations in the lab-frame mode frequency. Using the robust optimization that follows, you can assess the impact of lab frame mode frequency fluctuations in the following section.

In the following cell, the parametric drive and laser detunings are specified as for the standard gates above, with the gate durations extended to accommodate robust solutions.

# Gate and control specifications
drive_g_list = [0, 2.9e3, 50e3]
laser_detuning_list = [5.9e6 - 2.93e3, 5.9e6 - 5.005e3, 5.9e6 - 50.901e3]
robust_duration_list = [390e-6, 295e-6, 150e-6]

# Optimization parameters
number_of_segments = 20
number_of_free_segments = int(np.ceil(number_of_segments / 2))
robust_results = []

for drive_g, laser_detuning, duration in zip(
    drive_g_list, laser_detuning_list, robust_duration_list
):

    sample_times = np.linspace(0, duration, 50)

    # Set up trap modulation transformations
    relative_detunings = np.array(mode_frequencies - laser_detuning)
    mode_rel_det = relative_detunings[0, 1]
    mode_rel_det_driven = np.sqrt(mode_rel_det ** 2 - drive_g ** 2)
    relative_detunings_driven = copy.deepcopy(relative_detunings)
    relative_detunings_driven[0, 1] = mode_rel_det_driven
    r = 1 / 4 * np.log((mode_rel_det + drive_g) / (mode_rel_det - drive_g))

    # Optimize the controls
    with qctrl.create_graph() as graph:

        drive_moduli = qctrl.operations.optimization_variable(
            count=number_of_free_segments,
            lower_bound=0,
            upper_bound=maximum_rabi_rate,
        )
        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, full_phases = reflect_signal(
            number_of_segments, drive_phases, drive_moduli
        )
        ion_drive = qctrl.operations.complex_pwc_signal(
            moduli=full_moduli,
            phases=full_phases,
            duration=duration,
            name="ion_drive",
        )
        # Drive transformation due to trap modulation
        ion_drive_mod = ion_drive * np.cosh(r) + qctrl.operations.conjugate(
            ion_drive
        ) * np.sinh(r)

        drives = [
            ion_drive_mod
        ] * number_of_ions  # Identical drives applied to each ion

        ms_phases = qctrl.operations.ms_phases(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings_driven,  # Relative detuning transformation
            sample_times=sample_times,
            name="phases",
        )

        ms_displacements = qctrl.operations.ms_displacements(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings_driven,  # Relative detuning transformation
            sample_times=sample_times,
            name="displacements",
        )

        infidelity = qctrl.operations.ms_infidelity(
            phases=ms_phases[-1],
            displacements=ms_displacements[-1],
            target_phases=target,
            name="infidelity",
        )
        robust_cost_term = qctrl.operations.ms_dephasing_robust_cost(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters,
            relative_detunings=relative_detunings_driven,
        )
        cost = infidelity + robust_cost_term
        cost.name = "cost"

    result_robust = qctrl.functions.calculate_optimization(
        graph=graph,
        optimization_count=10,
        cost_node_name="cost",
        output_node_names=[
            "displacements",
            "infidelity",
            "cost",
            "ion_drive",
            "phases",
        ],
    )
    robust_results.append(result_robust)
    print("Drive strength (kHz):", drive_g / 1e3)
    print("Gate duration (\N{MICRO SIGN}s):", duration * 1e6)
    print("Infidelity:", result_robust.output["infidelity"]["value"])
    print("Robust cost:", result_robust.cost)
Your task calculate_optimization has completed.
Drive strength (kHz): 0.0
Gate duration (µs): 390.0
Infidelity: 0.0010925187871255737
Robust cost: 0.002435084472769894
Your task calculate_optimization has completed.
Drive strength (kHz): 2.9
Gate duration (µs): 295.0
Infidelity: 0.0017940433096121922
Robust cost: 0.007383165603800266
Your task calculate_optimization has completed.
Drive strength (kHz): 50.0
Gate duration (µs): 150.0
Infidelity: 5.867138350312828e-05
Robust cost: 9.561868459156533e-05

We extract the controls and transformed phase-space displacement plots, as above.

control_index = 2
control = {f"$\gamma$": robust_results[control_index].output["ion_drive"]}
plot_controls(plt.figure(), control)

# Sum over ion index to obtain total displacement of the mode.
trajl = np.sum(robust_results[control_index].output["displacements"]["value"], axis=3)

colors = ["#FD8F35", "#3273DC"]
scale = 1.1 * np.max(np.abs(trajl))

fig, ax = plt.subplots(figsize=(5, 5))
fig.suptitle("Phase space trajectories", fontsize=16, y=1.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 % 2),
        color=colors[mode],
        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.legend(loc="best", fontsize=14)
plt.show()

The complex control for the larger parametric drive (50 kHz) displayed above (chosen as an example) exhibits a Rabi rate close to the maximum value for most of its duration, while the phase exhibits more modulation over time. This control produces the symmetric transformed phase-space trajectory that satisfies the robustness condition for the transformed mode.

Verifying robustness to dephasing errors with a quasi-static scan

We assess the robustness of the 'Robust' optimized pulses in comparison to the 'Standard' optimized pulses using 1D quasi-static error-susceptibility scans. This involves calculating the infidelity of the controls over a systematic scan of the relative detunings (which corresponds to shifting the laser detuning). Note that this shift is applied to the lab-frame modes, and the effect of this noise can then be calculated by proceeding with the mode transformations.

maxshift = 500  # Hz
dephasing_scan_points = 31

shifts = np.linspace(-maxshift, maxshift, dephasing_scan_points)
results = standard_results + robust_results
drives_g = drive_g_list + drive_g_list
laser_detunings = laser_detuning_list + laser_detuning_list
durations = duration_list + robust_duration_list

name_list = []
for robust_name in ["standard_", "robust_"]:
    for name in ["no_g", "low_g", "high_g"]:
        name_list.append(
            [
                robust_name + name + str(number)
                for number in range(dephasing_scan_points)
            ]
        )

with qctrl.create_graph() as graph:
    for result, infidelity_names, drive_g, laser_detuning, duration in zip(
        results, name_list, drives_g, laser_detunings, durations
    ):
        for shift, infidelity_name in zip(shifts, infidelity_names):
            ion_drive = qctrl.operations.pwc_signal(
                values=np.array(
                    [segment["value"] for segment in result.output["ion_drive"]]
                ),
                duration=np.sum(
                    [segment["duration"] for segment in result.output["ion_drive"]]
                ),
            )

            # Set up trap modulation transformations
            relative_detunings = np.array(mode_frequencies - laser_detuning) + shift
            mode_rel_det = relative_detunings[0, 1]
            mode_rel_det_driven = np.sqrt(mode_rel_det ** 2 - drive_g ** 2)
            relative_detunings_driven = copy.deepcopy(relative_detunings)
            relative_detunings_driven[0, 1] = mode_rel_det_driven
            r = 1 / 4 * np.log((mode_rel_det + drive_g) / (mode_rel_det - drive_g))

            ion_drive_mod = ion_drive * np.cosh(r) + qctrl.operations.conjugate(
                ion_drive
            ) * np.sinh(r)
            drives = [ion_drive_mod] * number_of_ions

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

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

            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=np.array(name_list).flatten(),
)
Your task calculate_graph has started.
Your task calculate_graph has completed.
detshifts = np.linspace(-maxshift, maxshift, dephasing_scan_points)
infidelity_list = [
    [dephasing_scan.output[name]["value"] for name in names] for names in name_list
]
robust_names = ["Standard ", "Robust "]

gs = gridspec.GridSpec(1, 3)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(4 * 5)
fig.suptitle("Quasi-static dephasing scans", fontsize=16, y=1.05)
colors = ["#514B4F", "#680CE9"]

ax = fig.add_subplot(gs[0])
ax.set_title("Parametric drive strength 0 kHz", fontsize=14)
for ind, infids in enumerate([infidelity_list[0], infidelity_list[3]]):
    d_ind = ind * 3
    ax.plot(
        detshifts / 1000,
        infids,
        label=robust_names[ind] + str(durations[d_ind] * 1e6) + " \N{MICRO SIGN}s",
        color=colors[ind],
        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.set_ylim((-0.05, 1.05))
ax.legend(loc="best", fontsize=14)

ax = fig.add_subplot(gs[1])
ax.set_title("Parametric drive strength 3 kHz", fontsize=14)
for ind, infids in enumerate([infidelity_list[1], infidelity_list[4]]):
    d_ind = ind * 3 + 1
    ax.plot(
        detshifts / 1000,
        infids,
        label=robust_names[ind] + str(durations[d_ind] * 1e6) + " \N{MICRO SIGN}s",
        color=colors[ind],
        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.set_ylim((-0.05, 1.05))
ax.legend(loc="best", fontsize=14)

ax = fig.add_subplot(gs[2])
ax.set_title("Parametric drive strength 50 kHz", fontsize=14)
for ind, infids in enumerate([infidelity_list[2], infidelity_list[5]]):
    d_ind = ind * 3 + 2
    ax.plot(
        detshifts / 1000,
        infids,
        label=robust_names[ind] + str(durations[d_ind] * 1e6) + " \N{MICRO SIGN}s",
        color=colors[ind],
        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.set_ylim((-0.05, 1.05))
ax.legend(loc="best", fontsize=14)
plt.show()

The broader high-fidelity region indicates the benefit of the robust optimized pulses when there is quasi-static dephasing noise, for each of the parametric drive strengths. The parametric drive impacts the mediating mode, and is treated by a transformation on the mode. Higher parametric drive strengths reduce the robustness to detuning (or similarly mode frequency) noise. This can be understood by the robustness condition using the transformed mode providing less robustness with respect to the lab-frame mode as the parametric drive grows.

Finally, notice from the legends that the robust gates are each faster than their standard counterparts that have lower (or no) parametric drive. Depending on the trap specifics and noise levels, the parametric drive strength could be chosen to enhance the gate time while giving rise only to acceptable noise susceptibility.

Creating Q-CTRL pulses for two-qubit gates mediated by multiple modes

Mølmer–Sørensen-type operations can be mediated by more than a single mode; indeed this is necessary for more complex multi-ion or parallel gate operations. This may also involve controls optimized individually for each ion to provide the necessary control degrees of freedom to perform the target operation. In this section, you can generate ion-specific modulated control drives that perform the target entangling gate mediated by multiple motional modes in the presence of confining potential modulation.

As before, consider a parametric trap drive modifying a single motional mode. Applying a normal mode transformation for this driven mode can again give rise to Mølmer–Sørensen dynamics, when combined with other transformations. As for the single-mode-mediated gate, the relative detunings are transformed, however the ion drives are not specific to the modified mode (the drives are coupled to each mode) and cannot be freely transformed. Instead, by enforcing that the ion drives are real, the Lamb–Dicke parameter terms associated with the modified motional mode can be transformed to recover the Mølmer–Sørensen Hamiltonian form.

In this section, we use Boulder Opal to explore the impact of the trap modulation on two-ion gates mediated by both modes.

First, specify the number and type of ions, as well as other trap and laser characteristics. Note that the maximum Rabi rate has been increased from the previous section to enhance the coupling to both modes, and this will also result in faster gates.

# Trap characteristics
number_of_ions = 2

# Laser characteristics
maximum_rabi_rate = 2 * np.pi * 20e3 * 10

# Derived setup parameters
ion_chain_properties = qctrl.functions.calculate_ion_chain_properties(
    atomic_mass=25,  # Mg ions
    ion_count=number_of_ions,
    # Center of mass frequencies
    radial_x_center_of_mass_frequency=5.9e6,
    radial_y_center_of_mass_frequency=6.1e6,
    axial_center_of_mass_frequency=1.5e6,
    # Laser wavevector
    radial_x_wave_number=(2 * np.pi) / 355e-9,
    radial_y_wave_number=0,
    axial_wave_number=0,
)

# 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 mode frequencies in the shape
# [<axis>, <collective_mode>]
mode_frequencies = 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],
    ]
)

# Target operation
target = np.array(
    [
        [0, 0],
        [np.pi / 4, 0],
    ]
)
Your task calculate_ion_chain_properties has completed.

Generating standard pulses

As above, here we demonstrate the fastest single-displacement-loop gate operations modified by the trap modulation. The following parameters have been chosen to characterize this operation, taking into account the transformations when a parametric drive with strength $g>0$ modifies a particular mode $p$. Here we specify that $p$ is the $x$-axis COM mode.

# Gate and control specifications
drive_g_list = [0, 25e3]
laser_detuning_list = [5.9e6 - 29.3e3, 5.9e6 - 46.75e3]
duration_list = [34.2e-6, 25.4e-6]

Next we optimize the drive control for the given parameters. In the following optimization, the real drives $\gamma_j (t) = \Omega_j (t)$ are separately optimized and individually address each ion $j$.

# Optimization parameters
number_of_segments = 20
number_of_free_segments = number_of_segments
standard_results = []
drive_names = ["ion_drive" + str(number) for number in range(number_of_ions)]

for drive_g, laser_detuning, duration in zip(
    drive_g_list, laser_detuning_list, duration_list
):

    sample_times = np.linspace(0, duration, 50)

    # Set up trap modulation transformations
    relative_detunings = np.array(mode_frequencies - laser_detuning)
    mode_rel_det = relative_detunings[0, 1]
    mode_rel_det_driven = np.sqrt(mode_rel_det ** 2 - drive_g ** 2)
    relative_detunings_driven = copy.deepcopy(relative_detunings)
    relative_detunings_driven[0, 1] = mode_rel_det_driven
    r = 1 / 4 * np.log((mode_rel_det + drive_g) / (mode_rel_det - drive_g))
    lamb_dicke_parameters_driven = copy.deepcopy(lamb_dicke_parameters)
    lamb_dicke_parameters_driven[0][1] *= np.sinh(r) + np.cosh(
        r
    )  # Transform the x-axis COM mode

    # Optimize the controls
    with qctrl.create_graph() as graph:

        drives = []
        for drive_name in drive_names:
            drive_moduli = qctrl.operations.optimization_variable(
                count=number_of_free_segments,
                lower_bound=-maximum_rabi_rate,
                upper_bound=maximum_rabi_rate,
            )
            drive_phases = np.full((number_of_free_segments), 0)
            drives.append(
                qctrl.operations.pwc_signal(
                    values=drive_moduli,
                    duration=duration,
                    name=drive_name,
                )
            )

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

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

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

    result_optimal = qctrl.functions.calculate_optimization(
        graph=graph,
        optimization_count=10,
        cost_node_name="infidelity",
        output_node_names=["displacements", "infidelity", "phases"] + drive_names,
    )
    standard_results.append(result_optimal)
    print("Drive strength (kHz):", drive_g / 1e3)
    print("Gate duration (\N{MICRO SIGN}s):", duration * 1e6)
    print("Infidelity:", result_optimal.cost)
Your task calculate_optimization has completed.
Drive strength (kHz): 0.0
Gate duration (µs): 34.199999999999996
Infidelity: 8.625211656010379e-11
Your task calculate_optimization has completed.
Drive strength (kHz): 25.0
Gate duration (µs): 25.400000000000002
Infidelity: 1.0805049179940518e-05
control = {
    f"$\gamma_0$": standard_results[-1].output[drive_names[0]],
    f"$\gamma_1$": standard_results[-1].output[drive_names[1]],
}
plot_controls(plt.figure(), control)

# Sum over ion index to obtain total displacement of the mode.
trajl = np.sum(standard_results[-1].output["displacements"]["value"], axis=3)

colors = ["#FD8F35", "#3273DC"]
scale = 1.1 * np.max(np.abs(trajl))

fig, ax = plt.subplots(figsize=(5, 5))
fig.suptitle("Phase space trajectories", fontsize=16, y=1.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 % 2),
        color=colors[mode],
        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.legend(loc="best", fontsize=14)
plt.show()

The real control with a parametric drive (25 kHz) is displayed above as an example. The control is close to the maximal Rabi rate throughout the gate, as expected for the minimal gate duration specified above. The transformed closed phase-space trajectory is also displayed for this control.

Although both modes can be displaced to mediate the gate, the minimal gate duration exploits the maximally-coupling (COM) mode to perform the operation as fast as possible. In the next section, we explore longer gate durations for robust operations.

Generating robust pulses

As for the gate mediated by a single mode, we optimize the gate operation to satisfy robustness conditions. The conditions apply to the transformed mode as well as any additional displaced modes. In the following cell, the parametric drive and laser detunings are specified as for the standard gates above, with the gate durations extended to accommodate robust solutions.

# Gate and control specifications
drive_g_list = [0, 25e3]
laser_detuning_list = [5.9e6 - 29.3e3, 5.9e6 - 46.75e3]
robust_duration_list = [50e-6, 40e-6]

# Optimization parameters
number_of_segments = 20
number_of_free_segments = int(np.ceil(number_of_segments / 2))
robust_results = []
drive_names = ["ion_drive" + str(number) for number in range(number_of_ions)]

for drive_g, laser_detuning, duration in zip(
    drive_g_list, laser_detuning_list, robust_duration_list
):

    sample_times = np.linspace(0, duration, 50)

    # Set up trap modulation transformations
    relative_detunings = np.array(mode_frequencies - laser_detuning)
    mode_rel_det = relative_detunings[0, 1]
    mode_rel_det_driven = np.sqrt(mode_rel_det ** 2 - drive_g ** 2)
    relative_detunings_driven = copy.deepcopy(relative_detunings)
    relative_detunings_driven[0, 1] = mode_rel_det_driven
    r = 1 / 4 * np.log((mode_rel_det + drive_g) / (mode_rel_det - drive_g))
    lamb_dicke_parameters_driven = copy.deepcopy(lamb_dicke_parameters)
    lamb_dicke_parameters_driven[0][1] *= np.sinh(r) + np.cosh(
        r
    )  # Transform the x-axis COM mode

    # Optimize the controls
    with qctrl.create_graph() as graph:

        drives = []
        for drive_name in drive_names:
            drive_moduli = qctrl.operations.optimization_variable(
                count=number_of_free_segments,
                lower_bound=-maximum_rabi_rate,
                upper_bound=maximum_rabi_rate,
            )
            drive_phases = np.full((number_of_free_segments), 0)
            full_moduli, full_phases = reflect_signal(
                number_of_segments, drive_phases, drive_moduli
            )
            drives.append(
                qctrl.operations.pwc_signal(
                    values=full_moduli,
                    duration=duration,
                    name=drive_name,
                )
            )

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

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

        infidelity = qctrl.operations.ms_infidelity(
            phases=ms_phases[-1],
            displacements=ms_displacements[-1],
            target_phases=target,
            name="infidelity",
        )
        robust_cost_term = qctrl.operations.ms_dephasing_robust_cost(
            drives=drives,
            lamb_dicke_parameters=lamb_dicke_parameters_driven,
            relative_detunings=relative_detunings_driven,
        )
        cost = infidelity + robust_cost_term
        cost.name = "cost"

    result_robust = qctrl.functions.calculate_optimization(
        graph=graph,
        optimization_count=10,
        cost_node_name="cost",
        output_node_names=["displacements", "infidelity", "cost", "phases"]
        + drive_names,
    )
    robust_results.append(result_robust)
    print("Drive strength (kHz):", drive_g / 1e3)
    print("Gate duration (\N{MICRO SIGN}s):", duration * 1e6)
    print("Infidelity:", result_robust.cost)
    print("Robust cost:", result_robust.cost)
Your task calculate_optimization has completed.
Drive strength (kHz): 0.0
Gate duration (µs): 50.0
Infidelity: 4.8766658333243474e-11
Robust cost: 4.8766658333243474e-11
Your task calculate_optimization has completed.
Drive strength (kHz): 25.0
Gate duration (µs): 40.0
Infidelity: 7.287937349315749e-11
Robust cost: 7.287937349315749e-11
control = {
    f"$\gamma_0$": robust_results[1].output[drive_names[0]],
    f"$\gamma_1$": robust_results[1].output[drive_names[1]],
}
plot_controls(plt.figure(), control)

# Sum over ion index to obtain total displacement of the mode.
trajl = np.sum(robust_results[1].output["displacements"]["value"], axis=3)

colors = ["#FD8F35", "#3273DC"]
scale = 1.1 * np.max(np.abs(trajl))

fig, ax = plt.subplots(figsize=(5, 5))
fig.suptitle("Phase space trajectories", fontsize=16, y=1.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 % 2),
        color=colors[mode],
        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.legend(loc="best", fontsize=14)
plt.show()

The real control (the phase takes values of 0 or $\pi$) with a parametric drive (25 kHz) is displayed above as an example. Here both modes exhibit displacement and thus mediate the gate, with symmetry imposed by the robustness conditions on both the transformed mode (mode 1) and the lab-frame mode (mode 0). The gate is still primarily mediated by the transformed and more strongly-coupled COM mode. Both modes are restored to zero displacement at the end of the operation.

Verifying robustness to dephasing errors with a quasi-static scan

We assess the robustness of the 'Robust' optimized pulses in comparison to the 'Standard' optimized pulses using 1D quasi-static scans. This involves calculating the infidelity of the controls over a systematic scan of the relative detunings (which corresponds to shifting the laser detuning). Note that this shift is applied to the lab-frame modes, and the effect of this noise can then be calculated by proceeding with the mode transformations.

maxshift = 1500  # Hz
dephasing_scan_points = 31

shifts = np.linspace(-maxshift, maxshift, dephasing_scan_points)
results = standard_results + robust_results
drives_g = drive_g_list + drive_g_list
laser_detunings = laser_detuning_list + laser_detuning_list
durations = duration_list + robust_duration_list

name_list = []
for robust_name in ["standard_", "robust_"]:
    for name in ["no_g", "g"]:
        name_list.append(
            [
                robust_name + name + str(number)
                for number in range(dephasing_scan_points)
            ]
        )

with qctrl.create_graph() as graph:
    for result, infidelity_names, drive_g, laser_detuning, duration in zip(
        results, name_list, drives_g, laser_detunings, durations
    ):
        for shift, infidelity_name in zip(shifts, infidelity_names):

            # Set up trap modulation transformations
            relative_detunings = np.array(mode_frequencies - laser_detuning) + shift
            mode_rel_det = relative_detunings[0, 1]
            mode_rel_det_driven = np.sqrt(mode_rel_det ** 2 - drive_g ** 2)
            relative_detunings_driven = copy.deepcopy(relative_detunings)
            relative_detunings_driven[0, 1] = mode_rel_det_driven
            r = 1 / 4 * np.log((mode_rel_det + drive_g) / (mode_rel_det - drive_g))
            lamb_dicke_parameters_driven = copy.deepcopy(lamb_dicke_parameters)
            lamb_dicke_parameters_driven[0][1] *= np.sinh(r) + np.cosh(
                r
            )  # Transform the x-axis COM mode

            # Extract optimized drives
            drives = []
            drives.append(
                qctrl.operations.pwc_signal(
                    values=np.array(
                        [segment["value"] for segment in result.output[drive_names[0]]]
                    ),
                    duration=np.sum(
                        [
                            segment["duration"]
                            for segment in result.output[drive_names[0]]
                        ]
                    ),
                )
            )
            drives.append(
                qctrl.operations.pwc_signal(
                    values=np.array(
                        [segment["value"] for segment in result.output[drive_names[1]]]
                    ),
                    duration=np.sum(
                        [
                            segment["duration"]
                            for segment in result.output[drive_names[1]]
                        ]
                    ),
                )
            )

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

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

            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=np.array(name_list).flatten(),
)
Your task calculate_graph has started.
Your task calculate_graph has completed.
detshifts = np.linspace(-maxshift, maxshift, dephasing_scan_points)
infidelity_list = [
    [dephasing_scan.output[name]["value"] for name in names] for names in name_list
]
robust_names = ["Standard ", "Robust "]

gs = gridspec.GridSpec(1, 2)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(3 * 5)
fig.suptitle("Quasi-static dephasing scans", fontsize=16, y=1.05)
colors = ["#514B4F", "#680CE9"]

ax = fig.add_subplot(gs[0])
ax.set_title("Parametric drive strength 0 kHz", fontsize=14)
for ind, infids in enumerate([infidelity_list[0], infidelity_list[2]]):
    d_ind = ind * 2
    ax.plot(
        detshifts / 1000,
        infids,
        label=robust_names[ind] + str(durations[d_ind] * 1e6)[:5] + " \N{MICRO SIGN}s",
        color=colors[ind],
        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.set_ylim((-0.001, 0.02))
ax.legend(loc="best", fontsize=14)

ax = fig.add_subplot(gs[1])
ax.set_title("Parametric drive strength 25 kHz", fontsize=14)
for ind, infids in enumerate([infidelity_list[1], infidelity_list[3]]):
    d_ind = ind * 2 + 1
    ax.plot(
        detshifts / 1000,
        infids,
        label=robust_names[ind] + str(durations[d_ind] * 1e6)[:5] + " \N{MICRO SIGN}s",
        color=colors[ind],
        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.set_ylim((-0.001, 0.02))
ax.legend(loc="best", fontsize=14)
plt.show()

The broader high-fidelity region indicates the benefit of the robust optimized pulses when there is quasi-static dephasing noise, both with and without parametric trap modulation. This robustness is achieved by amplitude modulation of the ion drives in this case. Here the impact of the parametric drive on the gate robustness is minimal, and the higher-amplitude ion drives and faster gates have greatly enhanced the robustness over the previous quasi-static scans above. The parametric drive permits faster high-fidelity standard and robust solutions, as for the gate mediated by a single motional mode.