# 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 `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 `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 `calculate_graph`

function to track how displacements and phases change during the gate operational time. Within a graph, you can include the `ms_phases`

and `ms_displacements`

functions 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 `calculate_graph`

to return the dynamics.

Note that the results of the `ms_displacement`

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 `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 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()
# Helper function
def load_var(file_name):
# Return a variable from a json file
f = open(file_name, "r+")
encoded = f.read()
decoded = jsonpickle.decode(encoded)
f.close()
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(
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,
)
```

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

```
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 `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)
with qctrl.create_graph() as graph:
# Collect optimized drives
drives = [
qctrl.operations.pwc(
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 = 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",
)
result_optimal = qctrl.functions.calculate_graph(
graph=graph, output_node_names=["phases", "displacements"]
)
```

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)
gs = gridspec.GridSpec(1, 2, wspace=0.1)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(2.0 * 5)
fig.suptitle("Phase space trajectories", fontsize=16, y=1.0)
ax = fig.add_subplot(gs[0])
ax.tick_params(labelsize=14)
for mode in range(ion_count):
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,
)
scale = 1.1 * np.max(np.abs(trajl))
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(ion_count):
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=0, 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$ 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 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_optimal.output["phases"]["value"]
# Plot phase dynamics
fig, ax = plt.subplots(figsize=(10, 5))
for (ion_1, ion_2) in ion_pairs:
ax.plot(
sample_times * 1e3, phases[:, ion_1, ion_2], label=f"Ion pair {ion_1, ion_2}"
)
ax.axhline(np.pi / 4, color="#680CEA", ls="--", label="$\pi/4$") # Target phase
ax.legend()
ax.set_xlabel("Time (ms)", fontsize=14)
ax.set_ylabel("Relative phase", fontsize=14)
ax.tick_params(labelsize=14)
plt.title("Relative phase dynamics", fontsize=16, y=1)
plt.show()
```

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.