How to calculate phase and motion dynamics for arbitrarily modulated Mølmer–Sørensen gates

Calculate the Mølmer–Sørensen gate evolution characteristics for trapped ions

The Q-CTRL Python package enables you to optimize, explore and validate the dynamics induced by Mølmer–Sørensen-type operations. In this notebook, we begin with an optimized control and show how to obtain the evolution dynamics of the entangling phase and the motional trajectories in phase space, using convenience functions designed for trapped-ion experiments.

You can learn how to optimize error-robust Mølmer–Sørensen gates for trapped ions in our user guide, or perform complex operations by designing robust, configurable, parallel gates for large trapped-ion arrays using our application note.

Mølmer–Sørensen gate dynamics 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 qctrl.functions.calculate_ion_chain_properties by providing basic parameters of your ion trap, specifically the mass of the atom, center of mass (COM) frequencies and laser wave numbers along different axes.

You can then collect the Lamb–Dicke parameters and relative detunings from the result, which will be used in the phase and displacement calculations. Note that the Lamb–Dicke parameters need to be a 3D array and the relative detunings need to be a 2D array, as stated in the notes of graph.ms_phases.

2. Calculate motional displacement trajectories and phase dynamics

Given the controls, you can obtain the gate dynamics. In particular, 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 also obtain the phase accumulation during the gate operation for each pair of ions, and the target phases should be achieved at the end of a successful operation.

You can utilize our flexible qctrl.functions.calculate_graph function to track how displacements and phases change during the gate operational time. Within a graph, you can include the graph.ms_phases and graph.ms_displacements operations that take the control drives and ion properties as inputs. You can simply provide the sample_times parameter with your desired time sampling for the evolution, and use qctrl.functions.calculate_graph to return the dynamics.

Note that the results of the graph.ms_displacements function describe contributions to each mode from each ion. To get the overall displacement for each mode, you need to sum over the ion dimension (see the notes of graph.ms_displacements for details).

Worked 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 examine the dynamics induced by drives that have been optimized to generate a relative phase of $\pi/4$ between each ion pair. 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.

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

from qctrl import Qctrl

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

# Helper function
def load_var(file_name):
    # Return a variable from a json file
    f = open(file_name, "r+")
    encoded =
    decoded = jsonpickle.decode(encoded)
    return decoded

data_folder = "./resources/trapped-ions-calculate-dynamics/"
# 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 ion properties
ion_chain_properties = qctrl.functions.calculate_ion_chain_properties(
Your task calculate_ion_chain_properties (action_id="607684") has completed.

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

# Collect relative detunings in the shape
# [<axis>, <collective_mode>]
relative_detunings = (
            [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

Import optimized drives

Next, we import the optimized drive data. The drives have been optimized to generate $\pi/4$ relative phase between each ion pair in a chain of three ions.

drive_data = load_var(data_folder + "optimized_drives")
plot_controls(plt.figure(), drive_data)

The drives are displayed above, where an individual drive, drive_j, is applied to each ion j with optimized modulation in the modulus and phase (or angle).

Calculate relative phase and motional displacement dynamics

Here we apply our qctrl.functions.calculate_graph function to obtain the motional displacement trajectories and the relative phases for the duration of the given drives, with specified sample_times.

drive_names = list(drive_data.keys())
duration = np.sum([seg["duration"] for seg in drive_data[drive_names[0]]])
sample_times = np.linspace(0, duration, 100)

graph = qctrl.create_graph()

# Collect optimized drives
drives = [
        durations=np.array([v["duration"] for v in drive_data[name]]),
        values=np.asarray([v["value"] for v in drive_data[name]]),
    for name in drive_names

ms_phases = graph.ms_phases(

ms_displacements = graph.ms_displacements(

result_optimal = qctrl.functions.calculate_graph(
    graph=graph, output_node_names=["phases", "displacements"]
Your task calculate_graph (action_id="607685") has completed.

Now we can calculate the contribution of each ion to the mode displacements, and plot the mode displacement trajectories.

trajl = np.sum(result_optimal.output["displacements"]["value"], axis=3)
plot_range = 1.1 * np.max(np.abs(trajl))

fig, axs = plt.subplots(1, 2, figsize=(10, 5))
fig.suptitle("Phase space trajectories", y=1.1)

for k in range(2):

    for mode in range(ion_count):
            np.real(trajl[:, k, mode]),
            np.imag(trajl[:, k, mode]),
            label=f"mode {mode % 3}",
            np.real(trajl[-1, k, mode]),
            np.imag(trajl[-1, k, mode]),

    axs[k].set_xlim(-plot_range, plot_range)
    axs[k].set_ylim(-plot_range, plot_range)


hs, ls = axs[0].get_legend_handles_labels()
fig.legend(handles=hs, labels=ls, loc="center", bbox_to_anchor=(0.5, 1), ncol=ion_count)


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$ and its corresponding annihilation operator $a_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 we obtain and visualize the relative phase dynamics for each pair of ions.

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

phases = result_optimal.output["phases"]["value"]

# Plot phase dynamics
fig, ax = plt.subplots(figsize=(10, 5))
for (ion_1, ion_2) in ion_pairs:
        sample_times * 1e3, phases[:, ion_1, ion_2], label=f"Ion pair {ion_2, ion_1}"
plt.plot([0, sample_times[-1] * 1e3], [np.pi / 4, np.pi / 4], "k--")
ax.set_xlabel("Time (ms)")
ax.set_ylabel("Relative phase")
fig.suptitle("Relative phase dynamics")

The figure shows the dynamics of the relative phase between the three pairs of ions, over the duration of the gate. Under the optimized controls, the relative phases evolve from 0 to the target phase of $\pi/4$ at the end of the operation (marked by the horizontal dashed line). This indicates a successful operation.