# Designing robust, configurable, parallel gates for large trapped-ion arrays

**Obtaining control solutions for parallel and specifiable multi-qubit gates using Q-CTRL pulses**

Boulder Opal provides a broad toolkit enabling the generation of high-fidelity entangling gates with trapped ions. In particular, it is possible to engineer parallel optimized gates with arbitrary user-defined phases on large ion arrays. This comes from adding independent addressing and individual tuning of control parameters in order to simultaneously manipulate the phase-space trajectory of multiple motional modes.

In this notebook, we demonstrate control solutions for trapped ions robust to frequency-detuning and pulse-timing errors. We will cover:

- Creating optimized, robust, parallel Mølmer–Sørensen gates over multiple ions
- Calculating phase space trajectories for multiple modes
- Simulating entangling phase accumulation in the gate between arbitrary ion pairs
- Validating performance by calculating quasi-static noise susceptibility

Ultimately we will demonstrate how to parallelize gates in order to speed-up quantum circuit execution in trapped-ion quantum computers, while adding robustness against common hardware error sources.

For further discussion see Numeric optimization for configurable, parallel, error-robust entangling gates in large ion registers, published as Advanced Quantum Technology 3, 11 (2020).

## Worked example: Creating Q-CTRL pulses for parallel 2-of-5 and 3-of-5 qubit gates

Using Boulder Opal, you can obtain control solutions for configurable multi-qubit gates. These gates can also be performed simultaneously (in parallel). In this example, you'll learn to perform simultaneous two-qubit and three-qubit entangling operations.

First specify the ion, trap and laser characteristics.

```
import matplotlib.pyplot as plt
import numpy as np
from qctrlvisualizer import QCTRL_STYLE_COLORS, get_qctrl_style, plot_controls
from qctrl import Qctrl
plt.style.use(get_qctrl_style())
# Start a session with the API
qctrl = Qctrl()
# Trap characteristics
number_of_ions = 5
# 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,
)
```

We 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. For more information, see the How to optimize error-robust Mølmer–Sørensen gates for trapped ions user guide.

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

To demonstrate simultaneous two- and three-qubit gates, configure the target operation 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 target relative phase for the ion pair (1, 0) is set to $\pi/4$ for maximal two-qubit entanglement. For the three-qubit gate, the same relative phase is specified for ion pairs (3, 2), (4, 2) and (4, 3).

```
# Operation duration in seconds
duration = 3e-4
# Target phases: element (j,k) is the relative entanglement for ions j and k (k<j)
target = np.array(
[
[0, 0, 0, 0, 0],
[np.pi / 4, 0, 0, 0, 0],
[0, 0, 0, 0, 0],
[0, 0, np.pi / 4, 0, 0],
[0, 0, np.pi / 4, np.pi / 4, 0],
]
)
```

### Optimize drives (optimal control solution)

To control the system, consider separately-tunable complex drives $\gamma_j (t) = \Omega_j e^{i \phi_j}$ that address each ion $j$ in the trap. The drives are piecewise-constant with amplitude and phase modulation. In the optimization, the gate infidelity quantifies the solution performance.

```
number_of_segments = 128
sample_times = np.linspace(0, duration, number_of_segments)
drive_names = ["ion_drive" + str(number) for number in range(number_of_ions)]
# Helper function for robustness
def reflect_signal(graph, 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 = graph.reverse(moduli, [0])
phases_diff_refl = graph.concatenate(
[np.array([0]), graph.reverse(phases_diff, [0])], 0
)
phases_extra = phases[-1] + graph.cumulative_sum(phases_diff_refl)
else:
moduli_refl = graph.reverse(moduli[:-1], [0])
phases_diff_refl = graph.reverse(phases_diff, [0])
phases_extra = phases[-1] + graph.cumulative_sum(phases_diff_refl)
moduli_comb = graph.concatenate([moduli, moduli_refl], 0)
phases_comb = graph.concatenate([phases, phases_extra], 0)
return moduli_comb, phases_comb
# 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,
sample_times,
robust,
drive_names,
):
graph = qctrl.create_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 = graph.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 = graph.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(
graph, number_of_segments, drive_phases, drive_moduli
)
drives.append(
graph.complex_pwc_signal(
moduli=full_moduli,
phases=full_phases,
duration=duration,
name=drive_name,
)
)
ms_phases = graph.ms_phases(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings,
sample_times=sample_times,
name="phases",
)
ms_displacements = graph.ms_displacements(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings,
sample_times=sample_times,
name="displacements",
)
infidelity = graph.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 = graph.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,
cost_node_name="cost",
output_node_names=["displacements", "infidelity", "phases"] + drive_names,
)
# Helper function for plotting phase space trajectories
def plot_phase_space_trajectories(total_displacement):
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
plot_range = 1.1 * np.max(np.abs(total_displacement))
fig.suptitle("Phase space trajectories")
for k in range(2):
for mode in range(number_of_ions):
axs[k].plot(
np.real(total_displacement[:, k, mode]),
np.imag(total_displacement[:, k, mode]),
label=f"mode {mode % number_of_ions}",
)
axs[k].plot(
np.real(total_displacement[-1, k, mode]),
np.imag(total_displacement[-1, k, mode]),
"kx",
markersize=15,
)
axs[k].set_xlim(-plot_range, plot_range)
axs[k].set_ylim(-plot_range, plot_range)
axs[k].set_aspect("equal")
axs[k].set_xlabel("q")
axs[0].set_title("$x$-axis")
axs[0].set_ylabel("p")
axs[1].set_title("$y$-axis")
axs[1].yaxis.set_ticklabels([])
hs, ls = axs[0].get_legend_handles_labels()
fig.legend(handles=hs, labels=ls, loc="center left", bbox_to_anchor=(1, 0.5))
plt.tight_layout()
plt.show()
```

```
result_optimal = optimization_with_different_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,
sample_times=sample_times,
robust=False,
drive_names=drive_names,
)
```

```
print("Cost = Infidelity = ", result_optimal.output["infidelity"]["value"])
```

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

The figure displays the optimized pulse modulus and phase dynamics for ion 0.

### Optimize drives (robust control solution)

To obtain robust controls, we impose symmetry on each drive and optimize such that the centre of mass of each mode's trajectory is at zero. This is the same procedure as for the two-qubit case described here. In this case, the solution performance is given by the sum of the gate infidelity and the robustness cost. The robust option is enabled by setting the flag `robust`

to `True`

.

```
# Robust optimization
result_robust = optimization_with_different_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,
sample_times=sample_times,
robust=True,
drive_names=drive_names,
)
```

```
print("Cost = ", result_robust.cost)
print("Infidelity = ", result_robust.output["infidelity"]["value"])
```

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

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

### Performance validation

#### Calculate phase space trajectories

Next, we visualize the trajectory of the center of a coherent state in (rotating) optical phase space for each mode. Note that there are $3N$ modes in general, where $N$ is the ion number. We address only the $2N$ transverse modes in the trap. The closure of these trajectories is a necessary condition for an effective operation. We first examine the (non-robust) standard optimized control, followed by the robust optimized control.

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

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.

Now we 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)
plot_phase_space_trajectories(trajl_r)
```

Again, the black crosses at the origin indicate no residual state-motional entanglement, which arises from satisfying the center of mass and symmetry conditions.

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

```
fig, axs = plt.subplots(2, 1, figsize=(10, 10))
fig.suptitle("Relative phase dynamics")
x_timing = sample_times * 1e3
target_phases = [np.pi / 4, 0]
target_phase_names = ["\pi/4", "0"]
for k in range(2):
target_phase = target_phases[k]
for ion1 in range(number_of_ions):
for ion2 in range(ion1):
if target[ion1][ion2] == target_phase:
axs[k].plot(x_timing, phases[:, ion1, ion2], label=f"{ion2, ion1}")
axs[k].set_yticks([0, np.pi / 4])
axs[k].set_yticklabels([0, r"$\pi/4$"])
axs[k].set_ylim((-0.1, 0.9))
axs[k].set_ylabel("Relative phase")
axs[k].legend(title="Ion pair")
axs[k].plot([0, x_timing[-1]], [target_phase, target_phase], "k--")
axs[k].set_title(f"Pairs with target phase $\psi = {target_phase_names[k]}$")
axs[0].xaxis.set_ticklabels([])
axs[1].set_xlabel("Time (ms)")
plt.show()
```

These figures display the relative phase dynamics for the duration of the gate operation. For clarity, different plots display the ion pairs with different target relative phases. The two plots are for target relative phases of $\pi/4$ and 0, marked with black horizontal dashed lines.

Next consider the phase accumulation for the robust optimized control.

```
phases_r = result_robust.output["phases"]["value"]
```

```
fig, axs = plt.subplots(2, 1, figsize=(10, 9))
fig.suptitle("Relative phase dynamics")
x_timing = sample_times * 1e3
target_phases = [np.pi / 4, 0]
target_phase_names = ["\pi/4", "0"]
for k in range(2):
target_phase = target_phases[k]
for ion1 in range(number_of_ions):
for ion2 in range(ion1):
if target[ion1][ion2] == target_phase:
axs[k].plot(x_timing, phases_r[:, ion1, ion2], label=f"{ion2, ion1}")
axs[k].set_yticks([0, np.pi / 4])
axs[k].set_yticklabels([0, r"$\pi/4$"])
axs[k].set_ylim((-0.1, 0.9))
axs[k].set_ylabel("Relative phase")
axs[k].legend(title="Ion pair")
axs[k].plot([0, x_timing[-1]], [target_phase, target_phase], "k--")
axs[k].set_title(f"Pairs with target phase $\psi = {target_phase_names[k]}$")
axs[0].xaxis.set_ticklabels([])
axs[1].set_xlabel("Time (ms)")
plt.show()
```

As above, these figures display the relative phase dynamics for the duration of the gate operation. The two plots are for target relative phases of $\pi/4$ and 0, and these values are marked with black horizontal dashed lines.

The optimized drives produce nontrivial relative phase dynamics for each ion pair since the ions are individually addressed. The entangling phase is mediated by coupling the ion pairs to shared motional modes, and the optimized drives (both standard and robust) provide the necessary degrees of freedom to achieve the different target relative phases for each ion pair.

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

We can assess the robustness of the 'Robust' optimized pulses in comparison to the 'Standard' optimized pulses using 1D quasi-static scans.

First, we calculate 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.004
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)
]
graph = qctrl.create_graph()
for shift, infidelity_name in zip(shifts, optimal_infidelity_names):
drives = [
graph.pwc_signal(
values=np.array(
[segment["value"] for segment in result_optimal.output[drive_name]]
),
duration=duration * (1 + shift),
)
for drive_name in drive_names
]
ms_phases = graph.ms_phases(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings,
)
ms_displacements = graph.ms_displacements(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings,
)
ms_infidelity = graph.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,
)
infids_timing = [timing_scan.output[name]["value"] for name in optimal_infidelity_names]
```

```
robust_infidelity_names = [
"robust_infidelity" + str(number) for number in range(timing_scan_points)
]
graph = qctrl.create_graph()
for shift, robust_infidelity_name in zip(shifts, robust_infidelity_names):
drives = [
graph.pwc_signal(
values=np.array(
[segment["value"] for segment in result_robust.output[drive_name]]
),
duration=duration * (1 + shift),
)
for drive_name in drive_names
]
ms_phases = graph.ms_phases(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings,
)
ms_displacements = graph.ms_displacements(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings,
)
ms_infidelity = graph.ms_infidelity(
phases=ms_phases,
displacements=ms_displacements,
target_phases=target,
name=robust_infidelity_name,
)
timing_scan = qctrl.functions.calculate_graph(
graph=graph,
output_node_names=robust_infidelity_names,
)
robust_infids_timing = [
timing_scan.output[name]["value"] for name in robust_infidelity_names
]
```

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 = 4 * minreldet / 10000 * 100
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)
]
graph = qctrl.create_graph()
for shift, infidelity_name in zip(shifts, optimal_infidelity_names):
drives = [
graph.pwc_signal(
values=np.array(
[segment["value"] for segment in result_optimal.output[drive_name]]
),
duration=duration,
)
for drive_name in drive_names
]
ms_phases = graph.ms_phases(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings + shift,
)
ms_displacements = graph.ms_displacements(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings + shift,
)
ms_infidelity = graph.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,
)
infids_dephasing = [
dephasing_scan.output[name]["value"] for name in optimal_infidelity_names
]
```

```
robust_infidelity_names = [
"robust_infidelity" + str(number) for number in range(dephasing_scan_points)
]
graph = qctrl.create_graph()
for shift, robust_infidelity_name in zip(shifts, robust_infidelity_names):
drives = [
graph.pwc_signal(
values=np.array(
[segment["value"] for segment in result_robust.output[drive_name]]
),
duration=duration,
)
for drive_name in drive_names
]
ms_phases = graph.ms_phases(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings + shift,
)
ms_displacements = graph.ms_displacements(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings + shift,
)
ms_infidelity = graph.ms_infidelity(
phases=ms_phases,
displacements=ms_displacements,
target_phases=target,
name=robust_infidelity_name,
)
dephasing_scan = qctrl.functions.calculate_graph(
graph=graph,
output_node_names=robust_infidelity_names,
)
robust_infids_dephasing = [
dephasing_scan.output[name]["value"] for name in robust_infidelity_names
]
```

```
fig, axs = plt.subplots(1, 2, figsize=(15, 5))
fig.suptitle("Quasi-static scans", y=1.1)
maxshift = 0.004
timeshifts = np.linspace(1 - maxshift, 1 + maxshift, timing_scan_points)
axs[0].set_title("Timing noise")
axs[0].plot(timeshifts, robust_infids_timing, label="Robust")
axs[0].plot(timeshifts, infids_timing, label="Standard")
axs[0].set_xlabel("Pulse timing scaling factor")
axs[0].set_ylabel("Infidelity")
maxshift = 4 * minreldet / 10000 * 100
detshifts = np.linspace(-maxshift, maxshift, dephasing_scan_points)
axs[1].set_title("Dephasing noise")
axs[1].plot(detshifts / 1000, robust_infids_dephasing, label="Robust")
axs[1].plot(detshifts / 1000, infids_dephasing, label="Standard")
axs[1].set_xlabel("Laser detuning shift $\Delta \delta$ (kHz)")
axs[1].set_ylabel("Infidelity")
hs, ls = axs[0].get_legend_handles_labels()
fig.legend(handles=hs, labels=ls, loc="center", bbox_to_anchor=(0.5, 1.0), ncol=2)
plt.show()
```

The broader high-fidelity regions indicate the benefit of the robust optimized pulses when there is quasi-static dephasing noise or noise in the control pulse timings. The additional robustness for this more general gate is in agreement with published experimental results, using related robustness methodology.

## Worked example: Creating smoothed parallel gates with user-defined entangling phases

In this example we highlight the flexibility of our graph-based computational framework to add several substantial constraints on top of our optimization.

- We generate six parallel pairwise "2-of-12" entangling operations on a 12-ion chain
- Each pair is assigned a different target entangling phase, increasing in steps of 0.2 for each pair
- The high-frequency components of the drive are removed using a sinc filter within the optimizer

```
# Trap characteristics
number_of_ions = 12
# 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.25e6,
# Laser wavevector
radial_x_wave_number=(2 * np.pi) / 355e-9,
radial_y_wave_number=(2 * np.pi) / 355e-9,
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 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 = 3e-4
# Define target phases: element (j,k) is the relative entanglement for ions j and k (k<j)
target = np.zeros((12, 12))
for ion1 in range(1, 12, 2):
ion2 = ion1 - 1
target[ion1][ion2] = ion1 / 10
```

### Generating optimized pulses

Consider separately-tunable complex drives $\gamma_j (t) = \Omega_j e^{i \phi_j}$ that address each ion $j$ in the trap. The drives are piecewise-constant with amplitude and phase modulation. Additionally, a low-pass (sinc) filter is incorporated into the optimizer to smooth the pulse, as required for many practical implementations.

You can set the number of optimized pulse segments and resampling segments for the smoothed pulse, as well as the sinc cutoff frequency (in Hz) in the cell below.

```
number_of_segments = 32
number_of_sample_segments = 128
sample_times = np.linspace(0, duration, number_of_sample_segments)
sinc_cutoff_frequency_drives = 2 * np.pi * 0.05e6
```

```
drive_names = ["ion_drive" + str(number) for number in range(number_of_ions)]
graph = qctrl.create_graph()
sinc_kernel = graph.sinc_convolution_kernel(sinc_cutoff_frequency_drives)
drives = []
for drive_name in drive_names:
drive_real_raw_values = graph.optimization_variable(
count=number_of_segments,
lower_bound=-maximum_rabi_rate,
upper_bound=maximum_rabi_rate,
)
drive_imag_raw_values = graph.optimization_variable(
count=number_of_segments,
lower_bound=-maximum_rabi_rate,
upper_bound=maximum_rabi_rate,
)
drive_raw = graph.pwc_signal(
values=drive_real_raw_values + 1j * drive_imag_raw_values, duration=duration
)
drive_filtered = graph.convolve_pwc(pwc=drive_raw, kernel=sinc_kernel)
drives.append(
graph.discretize_stf(
stf=drive_filtered,
duration=duration,
segment_count=number_of_sample_segments,
name=drive_name,
)
)
ms_phases = graph.ms_phases(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings,
)
ms_displacements = graph.ms_displacements(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings,
)
infidelity = graph.ms_infidelity(
phases=ms_phases,
displacements=ms_displacements,
target_phases=target,
name="infidelity",
)
```

```
result = qctrl.functions.calculate_optimization(
graph=graph,
cost_node_name="infidelity",
output_node_names=["infidelity"] + drive_names,
)
```

```
print("Cost = Infidelity = ", result.output["infidelity"]["value"])
```

```
control = {f"$\gamma$": result.output["ion_drive0"]}
plot_controls(plt.figure(), control)
```

The above figure displays the optimized smooth pulse modulus and phase dynamics for ion 0, as an example of the ion-specific drives.

### Performance validation

#### Calculate phase space trajectories

For a system as large as this, you can separate the simulation of the displacement from the optimization to increase the efficiency. This procedure is demonstrated in the How to calculate phase and motion dynamics for arbitrarily modulated Mølmer–Sørensen gates User guide.

```
graph = qctrl.create_graph()
drives = [
graph.pwc_signal(
values=np.array([segment["value"] for segment in result.output[drive_name]]),
duration=duration,
)
for drive_name in drive_names
]
ms_displacements = graph.ms_displacements(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings,
sample_times=sample_times,
name="displacements",
)
displacement_simulation = qctrl.functions.calculate_graph(
graph=graph, output_node_names=["displacements"]
)
```

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

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.

```
graph = qctrl.create_graph()
drives = [
graph.pwc_signal(
values=np.array([segment["value"] for segment in result.output[drive_name]]),
duration=duration,
)
for drive_name in drive_names
]
ms_phases = graph.ms_phases(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings,
sample_times=sample_times,
name="phases",
)
phase_simulation = qctrl.functions.calculate_graph(
graph=graph, output_node_names=["phases"]
)
phases_6 = phase_simulation.output["phases"]["value"]
```

```
fig, ax = plt.subplots(figsize=(10, 5))
fig.suptitle("Relative phase dynamics")
x_timing = sample_times * 1e3
for ind, ion1 in enumerate(range(1, 12, 2)):
ion2 = ion1 - 1
ax.plot(
x_timing,
np.array(phases_6[:, ion1, ion2]),
label=f"{ion2, ion1}",
color=QCTRL_STYLE_COLORS[ind],
)
ax.plot(
[0, x_timing[-1]],
[target[ion1][ion2], target[ion1][ion2]],
"--",
color=QCTRL_STYLE_COLORS[ind],
)
ax.set_xlabel("Time (ms)")
ax.set_ylabel("Relative phase")
ax.legend(title="Ion pair", loc="center left", bbox_to_anchor=(1, 0.5))
plt.show()
```

The figure displays the relative phase dynamics for each entangled ion pair, with a color matched to the pair's target relative phase (dashed). Note that by the end of the operation, each ion pair achieves its specified, distinct relative phase target.