How to optimize error-robust Mølmer–Sørensen gates for trapped ions

Efficient state preparation using Mølmer–Sørensen-type interactions with in-built convenience functions

The Q-CTRL Python package enables you to efficiently prepare quantum states of trapped ion qubits using Mølmer–Sørensen-type interactions. In this notebook we show how to find both optimal and robust control drives to create a (multi-particle) entangled state using convenience functions designed for trapped-ion experiments.

The interaction Hamiltonian for Mølmer–Sørensen-type operations is in the form of

$$ H(t) = i\frac{\hbar}{2} \sum_{j=1}^N \sigma_{x, j} \sum_{p=1}^{3N} \eta_{pj} \left( \gamma_{j}(t) a_{p}^\dagger - \gamma^\ast_{j}(t) a_p \right) $$

where $N$ is the number of ions, $\sigma_x$ is the Pauli $X$ operator, $a_p$ is the annihilation operator for the mode $p$, $\eta_{pj}$ is the Lamb–Dicke parameter of mode $p$ for the ion $j$, and $\gamma_j(t)$ is the drive applied ion the ion $j$.

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. For these 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). $$

The objective of the gate-design process is to find definitions for the different drives $\gamma_j(t)$ which ensure high-fidelity gates which decouple all motional modes from the qubit degree of freedom. Gates may be designed to accumulate entangling phase either using interaction with only a single mode, or multiple modes. Note that multiple entangling operations can be performed in parallel, including multi-qubit entangling operations.

Summary Mølmer–Sørensen gate design workflow

1. Define and calculate ion properties

To find optimized drives, you first need to know the physical properties of the ions, namely the Lamb–Dicke parameters and mode frequencies. You can calculate these values using the function calculate_ion_chain_properties by providing basic parameters of your ion trap:

qctrl.functions.calculate_ion_chain_properties(
    atomic_mass, ion_count, radial_x_center_of_mass_frequency,
    radial_y_center_of_mass_frequency, axial_center_of_mass_frequency,
    radial_x_wave_number, radial_y_wave_number, axial_wave_number,
)

You 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.

2. Optimize drives

Next you need to set up the optimization properties. Define drives accessible in your system (e.g. individual-ion beams with programmable moduli and phases) as described in the How to optimize controls in D-dimensional systems using graphs User guide. The drives are typically piecewise-constant, where the amplitude and phase take constant values for each pulse segment of the given laser drive. A programmable drive can be specified within a graph as follows:

drive_moduli = qctrl.operations.optimization_variable(
    count, lower_bound, upper_bound
)
drive_phases = qctrl.operations.optimization_variable(
    count, lower_bound, upper_bound, is_lower_unbounded, is_upper_unbounded
)
ion_drive = qctrl.operations.complex_pwc_signal(
    moduli=drive_moduli, phases=drive_phases, duration=duration
)

Using the Lamb–Dicke parameters and relative detunings calculated above, particular drives give rise to relative phases and phase space displacements during an operation. You can calculate these gate properties using qctrl.operations.ms_phases and qctrl.operations.ms_displacements, which can be obtained for a specified array of sample_times:

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

The operational infidelity of the gate is determined the relative phases and displacements at the final time of the drives. You can obtain the infidelity using qctrl.operations.ms_infidelity:

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

Note that you need to format the target phases as a strictly lower triangular matrix (see the notes of qctrl.operations.ms_phases for details). With this infidelity as the cost function, you can find the optimal drives using the Q-CTRL optimization engine.

Visualize your optimized drives using plot_controls.

Note: a piecewise-constant control pulse can produce control frequencies that bring undesired transitions into resonance. To suppress these transitions for piecewise-constant controls, a low-pass filter can be applied to a given solution to eliminate these high-frequency components, typically with negligible fidelity loss. Alternatively, the segment durations of piecewise-constant controls can be selected to be multiples of the carrier detuning period. More generally, smooth controls can be optimized such that control frequency components lie below the carrier detuning, for instance using a sinc filter as demonstrated in the Designing robust, configurable, parallel gates for large trapped-ion arrays Application Note, and in the How to add smoothing and band-limits to control optimization User guide.

3. Optimize drives for error-robustness

You may also seek to produce error-robust gates which exhibit resilience against common noise sources (e.g. laser-detuning noise). Robustness to these noise sources is given by imposing a combination of symmetry in each ion's drive, as described in (Milne et al., Phys. Rev. Applied, 2020), and optimization such that the centre of mass of each mode's trajectory is at zero. A detailed discussion is given in Bentley et al. (2020), and this robust control protocol to create the Bell state has also been recently demonstrated experimentally, as shown in Milne et al. (2020).

You can optimize robust operations by both imposing symmetry in the drives, as demonstrated in the following examples, and by including the optimization cost term qctrl.operations.ms_dephasing_robust_cost to calculate the final displacement of each mode's centre of mass and use it in the Q-CTRL optimization engine to find the robust controls:

robust_cost_term = qctrl.operations.ms_dephasing_robust_cost(
    drives, lamb_dicke_parameters, relative_detunings,
)
cost = infidelity + robust_cost_term

4. Validate performance

Calculate phase-space displacement, relative phase, and infidelity dynamics

Given the controls, you can visualize the trajectory of the center of a coherent state in (rotating) optical phase space for each mode. The closure of these trajectories is a necessary condition for an effective operation.

You can return the mode displacements from the qctrl.operations.ms_displacements function within your optimization, for the specified sample_times. The results of the qctrl.operations.ms_displacements function describe contributions to each mode's displacement from each ion. To get the overall displacement for each mode, you need to sum over the ion dimension (see the notes of qctrl.operations.ms_displacements for details):

trajectories = np.sum(result.output["displacements"]["value"], axis=3)

To change the sample_times with a given optimized control, you can calculate the qctrl.operations.ms_displacements function within our flexible qctrl.functions.calculate_graph function, as shown in the How to calculate phase and motion dynamics for arbitrarily modulated Mølmer–Sørensen gates User guide.

The same procedures can be applied to track how relative phases for each ion pair and infidelity change during the gate operational time, using qctrl.operations.ms_phases and qctrl.operations.ms_infidelity, respectively. For high-fidelity operations, the relative phase obtained at the end of the operation should match the target phase for each ion pair.

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

You can assess the robustness of the 'Robust' optimized pulses in comparison to the 'Standard' optimized pulses using 1D quasi-static scans. By introducing a quasi-static offset of your choice, such as detuning, you can adapt and calculate the qctrl.operations.ms_displacements and qctrl.operations.ms_phases terms in the graph, which determine the error-dependent qctrl.operations.ms_infidelity.

Worked example: Optimal and robust Mølmer–Sørensen gates for two-qubit gates

Here we explore a simple example in which we define optimal and robust control solutions implementing a Mølmer–Sørensen gate on a pair of Ytterbium ions.

Define and calculate ion properties

For multi-qubit operations, first specify the number and type of ions, as well as other trap and laser characteristics.

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

plt.style.use(get_qctrl_style())

from qctrl import Qctrl

# Starting a session with the API
qctrl = Qctrl()
# Trap characteristics
number_of_ions = 2

# 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.

You 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.

# 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 = 2e-4

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

Now we establish the graph that represents this system, which is used in the optimization itself. Note that in this example we have established functions to calculate both optimal and robust gates in the same codeblock.

# 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 identical drives in all ions
def optimization_with_identical_drives(
    ion_count,
    number_of_segments,
    duration,
    maximum_rabi_rate,
    lamb_dicke_parameters,
    relative_detunings,
    target,
    sample_times,
    robust,
):
    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))

        # 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
            )

        ion_drive = qctrl.operations.complex_pwc_signal(
            moduli=full_moduli, phases=full_phases, duration=duration, name="ion_drive"
        )

        drives = [ion_drive] * ion_count

        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", "ion_drive", "phases"],
    )


# 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,
    target,
    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,
    )

Optimize drives (optimal control solution)

Next, specify the control pulses according to the degrees of freedom in the hardware. In this optimization, the complex drive $\gamma (t) = \Omega e^{i \phi}$ jointly addresses both ions in the chain. Alternatively, separate drives can be optimized to individually address each ion; graph encodings for both approaches are provided above. Here, the drives are defined by optimizable moduli and phases for individual segments (but note that, in general, you can define the drives using any of the techniques covered in the How to optimize controls in D-dimensional systems using graphs User guide)

The optimizer attempts to minimize a particular cost. The cost is specified in the below cell to quantify the distances of the trajectory end-points from the origin, and the differences between the realized and target phases. The infidelity also quantifies the solution performance.

number_of_segments = 64

sample_times = np.linspace(0, duration, number_of_segments)
result_optimal = optimization_with_identical_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,
    target=target,
    sample_times=sample_times,
    robust=False,
)
Your task calculate_optimization has completed.

print(f"Cost = Infidelity = {result_optimal.output['infidelity']['value']:.3e}")
Cost = Infidelity = 5.879e-12
control = {f"$\gamma$": result_optimal.output["ion_drive"]}
plot_controls(plt.figure(), control)

The above figure displays the optimized pulse modulus $\Omega (t)$ and phase $\phi (t)$ that addresses both ions.

Optimize drives (robust control solution)

Next we show how to optimize the drives to obtain the maximally-entangled state for two ions affected by dephasing noise, which is common in experiments due to various noise sources, for example frequency fluctuations in control lasers. The optimization procedure is extended by adding the robustness symmetry constraint and calculation of qctrl.operations.ms_dephasing_robust_cost in the graph-based optimization; this is enabled in the helper function defined above by setting robust=True.

# Robust optimization
result_robust = optimization_with_identical_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,
    target=target,
    sample_times=sample_times,
    robust=True,
)
Your task calculate_optimization has completed.

print(f"Cost = {result_robust.cost:.3e}")
print(f"Infidelity = {result_robust.output['infidelity']['value']:.3e}")
Cost = 3.960e-11
Infidelity = 1.488e-12
ion = 0
control = {f"$\gamma$": result_robust.output["ion_drive"]}
plot_controls(plt.figure(), control)

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

Validate performance

Calculate phase space trajectories

# 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.15)

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 % 2),
        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 % 2),
        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=(0.14, 1.3), fontsize=14)

plt.show()

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$. 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. The z-axis modes are not addressed, and thus have no excursions in phase space.

Next 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.15)

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 % 2),
        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 % 2),
        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=(0.14, 1.3), fontsize=14)

plt.show()

Here, the center of mass (the value of the integrated trajectory) for each mode is optimized to be close to the origin. 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

phases = result_optimal.output["phases"]["value"]
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)
ion1 = 1
ion2 = 0
x_timing = sample_times * 1e3
ax.plot(
    x_timing,
    phases[:, ion1, ion2],
    label="(" + str(ion1) + ", " + str(ion2) + ")",
    linewidth=2,
)
ax.plot([0, x_timing[-1]], [np.pi / 4, np.pi / 4], "k--", linewidth=0.8)
ax.plot([x_timing[-1], x_timing[-1]], [-1, 1], "k", alpha=0.5)
plt.yticks([0, np.pi / 4], [r"0", r"$\frac{\pi}{4}$"], fontsize=14)
plt.ylim((-0.03, 0.85))
plt.xlabel("Time (ms)", fontsize=14)
plt.ylabel("Relative phase", fontsize=14)
ax.legend(
    ncol=1,
    loc="best",
    bbox_to_anchor=(0.2, 0.8),
    fontsize=14,
    title="Ion pair",
    title_fontsize=14,
)
plt.show()

The figure shows the dynamics of the relative phase between the two ions, over the duration of the gate. Under the optimized control, this evolves from 0 to the target maximally-entangled phase of $\pi/4$ at the end of the operation (marked by the black vertical line). The target phase is marked by the horizontal dashed black line.

Next obtain the phase accumulation for the robust optimized control.

phases_r = result_robust.output["phases"]["value"]
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)
ion1 = 1
ion2 = 0
x_timing = sample_times * 1e3
ax.plot(
    x_timing,
    phases_r[:, ion1, ion2],
    label="(" + str(ion1) + ", " + str(ion2) + ")",
    linewidth=2,
)
ax.plot([0, x_timing[-1]], [np.pi / 4, np.pi / 4], "k--", linewidth=0.8)
ax.plot([x_timing[-1], x_timing[-1]], [-1, 1], "k", alpha=0.5)
plt.yticks([0, np.pi / 4], [r"0", r"$\frac{\pi}{4}$"], fontsize=14)
plt.ylim((-0.03, 0.85))
plt.xlabel("Time (ms)", fontsize=14)
plt.ylabel("Relative phase", fontsize=14)
ax.legend(
    ncol=1,
    loc="best",
    bbox_to_anchor=(0.2, 0.8),
    fontsize=14,
    title="Ion pair",
    title_fontsize=14,
)
plt.show()

By the end of the robust operation (black vertical line), the relative phase of the ions has reached the target value of $\pi/4$, marked by the dashed black line.

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

You can assess the robustness of the 'Robust' optimized pulses in comparison to the 'Standard' optimized pulses using 1D quasi-static scans. First consider 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.02

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)
]
robust_infidelity_names = [
    "robust_infidelity" + str(number) for number in range(timing_scan_points)
]

with qctrl.create_graph() as graph:
    for result, infidelity_names in [
        (result_optimal, optimal_infidelity_names),
        (result_robust, robust_infidelity_names),
    ]:
        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=duration * (1 + shift),
            )

            drives = [ion_drive] * 2

            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 + robust_infidelity_names
)

infids_timing = [timing_scan.output[name]["value"] for name in optimal_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 = 4000 * minreldet / 10000

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)
]
robust_infidelity_names = [
    "robust_infidelity" + str(number) for number in range(dephasing_scan_points)
]

with qctrl.create_graph() as graph:
    for result, infidelity_names in [
        (result_optimal, optimal_infidelity_names),
        (result_robust, robust_infidelity_names),
    ]:
        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=duration,
            )

            drives = [ion_drive] * 2

            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 + robust_infidelity_names
)

infids_dephasing = [
    dephasing_scan.output[name]["value"] for name in optimal_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.02
timeshifts = np.linspace(1 - maxshift, 1 + maxshift, timing_scan_points)
ax.plot(timeshifts, robust_infids_timing, label="Robust", linewidth=2)
ax.plot(timeshifts, infids_timing, label="Standard", 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 = 4000 * minreldet / 10000
detshifts = np.linspace(-maxshift, maxshift, dephasing_scan_points)
ax.plot(detshifts / 1000, robust_infids_dephasing, label="Robust", linewidth=2)
ax.plot(detshifts / 1000, infids_dephasing, label="Standard", 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 region indicates the benefit of the robust optimized pulses when there is quasi-static dephasing noise or noise in the control pulse timings, demonstrating the utility of the robust solution.

Example: GHZ state on a chain of three ions

In this example we consider a chain of three ions (${}^{171}{\rm Yb}^{+}$), and the goal is to find a set of drives such that, at the end of the drives, the relative phase between each ion pair is $\pi/4$. Note that you can perform $X_{\pi/2}$ rotations for all ions after applying the optimal drives to create the GHZ state, as shown in Kim et al. (2009).

Define and calculate ion properties

We begin by calculating ion properties based on inputs relevant to the system under examination.

# Specify ion trap parameters
ion_count = 3
atomic_mass = 171

# COM 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 characteristics
radial_x_wave_number = 1.8e7
radial_y_wave_number = 1.8e7
axial_wave_number = 0

laser_detuning = 1.6047e6
maximum_rabi_rate = 2 * np.pi * 100e3
# Calculate the ions properties
ion_chain_properties = qctrl.functions.calculate_ion_chain_properties(
    atomic_mass=atomic_mass,
    ion_count=ion_count,
    radial_x_center_of_mass_frequency=radial_x_center_of_mass_frequency,
    radial_y_center_of_mass_frequency=radial_y_center_of_mass_frequency,
    axial_center_of_mass_frequency=axial_center_of_mass_frequency,
    radial_x_wave_number=radial_x_wave_number,
    radial_y_wave_number=radial_y_wave_number,
    axial_wave_number=axial_wave_number,
)
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
)

Optimize drives

In this example, the drive for each ion is defined by optimizable moduli and phases for individual segments (but note that, in general, you can define the drives using any of the techniques covered in the How to optimize controls in D-dimensional systems using graphs User guide). In this example we only seek an optimal solution rather than a robust solution.

# Set up the target phases as a strictly lower triangular matrix
target_phases = np.asarray([[0, 0, 0], [np.pi / 4, 0, 0], [np.pi / 4, np.pi / 4, 0]])

# Number of segments of the drives
segment_count = 128

# Total durations of the drives
duration = 2e-4

# Drive node names for retrieving the results
drive_names = ["drive_" + str(i) for i in range(ion_count)]

optimization_result = optimization_with_different_drives(
    ion_count=ion_count,
    number_of_segments=segment_count,
    duration=duration,
    maximum_rabi_rate=maximum_rabi_rate,
    lamb_dicke_parameters=lamb_dicke_parameters,
    relative_detunings=relative_detunings,
    target=target_phases,
    sample_times=np.array([duration]),
    robust=False,
    drive_names=drive_names,
)
Your task calculate_optimization has started.
Your task calculate_optimization has completed.

print(f"Resulting phases: \n {optimization_result.output['phases']['value']}")
print(f"Target phases: \n {target_phases}")
print("\nOptimal control for ion 0:")
plot_controls(plt.figure(), {"drive_0": optimization_result.output["drive_0"]})
Resulting phases: 
 [[0.         0.         0.        ]
 [0.78540307 0.         0.        ]
 [0.78540164 0.78540332 0.        ]]
Target phases: 
 [[0.         0.         0.        ]
 [0.78539816 0.         0.        ]
 [0.78539816 0.78539816 0.        ]]

Optimal control for ion 0:

Calculate gate dynamics

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

with qctrl.create_graph() as graph:

    # Collect optimal controls from optimization
    optimal_drives = [
        qctrl.operations.pwc(
            durations=np.array(
                [v["duration"] for v in optimization_result.output[name]]
            ),
            values=np.asarray([v["value"] for v in optimization_result.output[name]]),
        )
        for name in drive_names
    ]

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

    infidelity = qctrl.operations.ms_infidelity(
        phases=phases,
        displacements=displacements,
        target_phases=target_phases,
        name="infidelity",
    )

result = qctrl.functions.calculate_graph(
    graph=graph, output_node_names=["displacements", "infidelity", "phases"]
)

mode_displacements = np.sum(result.output["displacements"]["value"], axis=-1)
Your task calculate_graph has completed.

fig, axes = plt.subplots(2, 2, figsize=(16, 12))

axes[0][0].plot(
    sample_times * 1e3, result.output["infidelity"]["value"], color="#680CEA"
)
axes[0][0].set_xlabel("Sample times (ms)")
axes[0][0].set_ylabel("Infidelity")

# Recall that the phases are stored as a strictly lower triangular matrix
# See the notes part of qctrl.operations.ms_phases
ion_pairs = []
for ion_1 in range(3):
    for ion_2 in range(ion_1):
        ion_pairs.append((ion_1, ion_2))
phases = result.output["phases"]["value"]

for (ion_1, ion_2) in ion_pairs:
    axes[0][1].plot(
        sample_times * 1e3, phases[:, ion_1, ion_2], label=f"ion pair {ion_1, ion_2}"
    )
# Target phase
axes[0][1].axhline(np.pi / 4, color="#680CEA", ls="--", label="$\pi/4$")
axes[0][1].legend()
axes[0][1].set_xlabel("Sample times (ms)")
axes[0][1].set_ylabel("Relative phase")

# Mode displacement trajectories for radial dimensions
for idx in range(ion_count):
    axes[1][0].plot(
        np.real(mode_displacements[:, 0, idx]),
        np.imag(mode_displacements[:, 0, idx]),
        label=f"mode_x{idx}",
    )
    axes[1][1].plot(
        np.real(mode_displacements[:, 1, idx]),
        np.imag(mode_displacements[:, 1, idx]),
        label=f"mode_y{idx}",
    )
# Displacement trajectories are expected to back to 0
# at the end of the drive as labelled by X here.
axes[1][0].plot(0, 0, marker="x", color="black", markersize=15)
axes[1][1].plot(0, 0, marker="x", color="black", markersize=15)
axes[1][0].legend()
axes[1][0].set_xlabel("q")
axes[1][0].set_ylabel("p")
axes[1][1].legend()
axes[1][1].set_xlabel("q")
axes[1][1].set_ylabel("p")
plt.show()

The above figure shows the infidelity and relative phases over the gate duration (top), as well as the motional phase space trajectories (bottom), which are closed loops indicating no final displacement for each mode.

Summary

In these examples, you have seen how the Q-CTRL Python package can easily calculate the important physical properties for your ion trap, efficiently find optimal or robust drives for state preparation, and track the ion states in phase space at any sample points during the gate operational time.

You can apply the same optimization strategy to large systems involving tens of qubits. This process is demonstrated in the Designing robust, configurable, parallel gates for large trapped-ion arrays Application note, where you will learn about how the Q-CTRL optimization engine enables you to efficiently find optimal and noise-robust control pulses for large trapped ion systems.