How to optimize Mølmer–Sørensen gates for a multitone global beam

Creating efficient gates for trapped ions without individually addressing the ions

Boulder Opal equips you to efficiently optimize trapped-ion gates that use Mølmer–Sørensen interactions. These gates use the vibrational modes of the ion chain to mediate the interaction between the ions, leaving the qubits in an entangled state and restoring the motional state after the operation.

Mølmer–Sørensen-type operations can be performed using either individual addressing of each ion, or a global beam that interacts with all the ions in the chain. The current user guide shows how to perform the optimization in the case of a global beam with multiple tones. In other words, the control problem we consider here is one where the interaction Hamiltonian is of the form:

\[H(t) = i\frac{\hbar}{2} \sum_{j=1}^N \sigma_{x, j} \sum_{p=1}^{3N} \sum_{k=1}^M \eta_{kpj} \left\{ \gamma^{(k)}(t) e^{i \delta_{kp} t} a_{p}^\dagger - \left[ \gamma^{(k)} (t) \right]^* e^{-i \delta_{kp} t} a_p \right\},\]

where $N$ is the number of ions, $M$ is the number of tones in the global beam, and $a_p$ is the annihilation operator for the p-th mode of oscillation of the chain. The parameters $\delta_{kp}$ indicate the relative detuning between the k-th tone of the global beam and the p-th mode of oscillation, and the $\eta_{kpj}$ are the Lamb–Dicke parameters for the k-th tone, p-th mode of oscillation, and j-th ion. These two set of parameters are provided by you as input for the optimization, in order to obtain the controls $\gamma^{(k)}(t)$ for the k-th tone of the global beam.

For more information about the theory behind Mølmer-Sørensen gates, see the user guide How to optimize error-robust Mølmer–Sørensen gates for trapped ions. There you find an introduction to the Mølmer–Sørensen gates and information about how to optimize them both for the case of a global single-tone beam, and for the more flexible case with individual ion addressing.

Summary workflow

1. Define and calculate ion trap properties

Pass the parameters of the ion chain to qctrl.functions.calculate_ion_chain_properties and obtain the frequencies and Lamb–Dicke parameters associated with each of the vibrational modes. Note that the Lamb–Dicke parameters and frequencies might be different for each of the tones of the global beam; these should be calculated using a separate function evaluation for each tone with a distinct wave number.

2. Optimize drives

Inside a graph, define individual optimizable drives for each of the tones of the system. Use the tones to calculate the relative phases obtained between each ion, and the displacements of each motional mode due to the drives. The node graph.ms_phases_multitone, which returns the relative phases, accepts Lamb-Dicke parameters and relative detunings for all the tones of the beam. In contrast, the node graph.ms_displacements calculates the displacements for one tone at a time. After all the displacements have been calculated, you can add them together to obtain the total mode displacements due to all the tones.

The performance or fidelity of the gate is determined by the final relative phases and mode displacements; the result of the phase and displacement calculations is thus passed to graph.ms_infidelity, in order to obtain the cost function for optimal controls. (You can add the contribution from graph.ms_dephasing_robust_cost if you want to build error-robust gates, as shown in the user guide How to optimize error-robust Mølmer–Sørensen gates for trapped ions.)

Run this optimization with qctrl.functions.calculate_optimization to obtain the optimal drives.

Worked example: Optimal controls to perform a Mølmer–Sørensen gate with a two-tone global beam.

The following example creates optimal two-tone global drives to maximally entangle a two-ion chain. At the end, the drives obtained are plotted.

import matplotlib.pyplot as plt
import numpy as np

from qctrl import Qctrl
from qctrlvisualizer import plot_controls

# Start a session with the API.
qctrl = Qctrl()
# Define parameters of the ion trap.
ion_count = 2

# Define laser parameters.
maximum_rabi_rate = 2 * np.pi * 50e3  # Hz
laser_detuning = 1.6e6 + 4.7e3  # Hz

# Define gate parameters.
duration = 3e-4  # s
target = np.array([[0, 0], [np.pi / 4, 0]])
segment_count = 64
# Obtain setup parameters for the first tone.
ion_chain_properties_1 = qctrl.functions.calculate_ion_chain_properties(
    atomic_mass=171,  # Yb ions
    ion_count=ion_count,
    # 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 * 2.0 / 355e-9,
    radial_y_wave_number=0,
    axial_wave_number=0,
)

lamb_dicke_parameters_1 = np.asarray(
    [
        [
            mode.lamb_dicke_parameters
            for mode in ion_chain_properties_1.radial_x_mode_properties
        ],
        [
            mode.lamb_dicke_parameters
            for mode in ion_chain_properties_1.radial_y_mode_properties
        ],
        [
            mode.lamb_dicke_parameters
            for mode in ion_chain_properties_1.axial_mode_properties
        ],
    ]
)

relative_detunings_1 = (
    np.asarray(
        [
            [
                mode.frequency
                for mode in ion_chain_properties_1.radial_x_mode_properties
            ],
            [
                mode.frequency
                for mode in ion_chain_properties_1.radial_y_mode_properties
            ],
            [mode.frequency for mode in ion_chain_properties_1.axial_mode_properties],
        ]
    )
    - laser_detuning
)
Your task calculate_ion_chain_properties (action_id="1118370") has completed.
# Obtain setup parameters for the second tone.
ion_chain_properties_2 = qctrl.functions.calculate_ion_chain_properties(
    atomic_mass=171,  # Yb ions
    ion_count=ion_count,
    # 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 * 2.4 / 355e-9,
    radial_y_wave_number=0,
    axial_wave_number=0,
)

lamb_dicke_parameters_2 = np.asarray(
    [
        [
            mode.lamb_dicke_parameters
            for mode in ion_chain_properties_2.radial_x_mode_properties
        ],
        [
            mode.lamb_dicke_parameters
            for mode in ion_chain_properties_2.radial_y_mode_properties
        ],
        [
            mode.lamb_dicke_parameters
            for mode in ion_chain_properties_2.axial_mode_properties
        ],
    ]
)

relative_detunings_2 = (
    np.asarray(
        [
            [
                mode.frequency
                for mode in ion_chain_properties_2.radial_x_mode_properties
            ],
            [
                mode.frequency
                for mode in ion_chain_properties_2.radial_y_mode_properties
            ],
            [mode.frequency for mode in ion_chain_properties_2.axial_mode_properties],
        ]
    )
    - laser_detuning
)
Your task calculate_ion_chain_properties (action_id="1118371") has completed.
# Define optimization graph.
graph = qctrl.create_graph()

# Define drive for the first tone.
drive_moduli_1 = graph.optimization_variable(
    count=segment_count, lower_bound=0, upper_bound=maximum_rabi_rate
)
drive_phases_1 = graph.optimization_variable(
    count=segment_count,
    lower_bound=0,
    upper_bound=2 * np.pi,
    is_lower_unbounded=True,
    is_upper_unbounded=True,
)
ion_drive_1 = graph.complex_pwc_signal(
    moduli=drive_moduli_1,
    phases=drive_phases_1,
    duration=duration,
    name=r"$\gamma^{(1)}$",
)

# Define drive for the second tone.
drive_moduli_2 = graph.optimization_variable(
    count=segment_count, lower_bound=0, upper_bound=maximum_rabi_rate
)
drive_phases_2 = graph.optimization_variable(
    count=segment_count,
    lower_bound=0,
    upper_bound=2 * np.pi,
    is_lower_unbounded=True,
    is_upper_unbounded=True,
)
ion_drive_2 = graph.complex_pwc_signal(
    moduli=drive_moduli_2,
    phases=drive_phases_2,
    duration=duration,
    name=r"$\gamma^{(2)}$",
)

# Calculate phases due to both drives.
ms_phases = graph.ms_phases_multitone(
    drives=[ion_drive_1, ion_drive_2],
    lamb_dicke_parameters=np.array([lamb_dicke_parameters_1, lamb_dicke_parameters_2]),
    relative_detunings=np.array([relative_detunings_1, relative_detunings_2]),
)

# Calculate displacements due to each drive, which is applied to both ions.
ms_displacements_1 = graph.ms_displacements(
    drives=[ion_drive_1] * ion_count,
    lamb_dicke_parameters=lamb_dicke_parameters_1,
    relative_detunings=relative_detunings_1,
)
ms_displacements_2 = graph.ms_displacements(
    drives=[ion_drive_2] * ion_count,
    lamb_dicke_parameters=lamb_dicke_parameters_2,
    relative_detunings=relative_detunings_2,
)

# Calculate total infidelity due to the phases and both displacement contributions.
infidelity = graph.ms_infidelity(
    phases=ms_phases,
    displacements=ms_displacements_1 + ms_displacements_2,
    target_phases=target,
    name="infidelity",
)
# Run optimization.
result = qctrl.functions.calculate_optimization(
    graph=graph,
    cost_node_name="infidelity",
    output_node_names=[r"$\gamma^{(1)}$", r"$\gamma^{(2)}$"],
)
Your task calculate_optimization (action_id="1118372") has completed.
# Print cost and plot results.
print(f"Optimization cost: {result.cost:.3e}")

plot_controls(plt.figure(), result.output)
Optimization cost: 3.109e-15

png