# Cold ions: obtaining robust, configurable multi-qubit gates

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

BOULDER OPAL enables you to optimize controls to achieve a target operation on your cold ion hardware. In this guide, we obtain and verify control solutions for trapped ions that are optimized for robustness to frequency detuning and pulse timing errors. The notebook is organized as follows: we introduce Q-CTRL control solutions for two qubits, experimentally validate our methodology, then provide examples to demonstrate our highly-configurable control solutions. Our user-specified target operations can include parallel multi-qubit gates and distinct configurable relative phases. Control constraints or degrees of freedom can also be user-specified, such as a shared drive with a maximum modulus slew rate. For each use-case we demonstrate amplitude- and phase-modulated solutions, and these degrees of freedom are similarly specifiable.

This notebook contains preview features that are not currently available in BOULDER OPAL. Please contact us if you want a live demonstration for your hardware. Preview features are included through the package qctrlcore, which is not publicly available: cells with qctrlcore will not be executable.

## Imports and initialization

import numpy as np
import tensorflow as tf
import time
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns
from qctrlvisualizer import plot_controls

# Predefined pulse imports
from qctrlopencontrols import new_predefined_driven_control

# Preview features
from qctrlcore import (optimize,
create_complex_pwc_signal,
OptimizationSpecification,
Pwc,
create_anchored_difference_bounded_variables,
)
from qctrlcore.ions import (get_default_modes,
get_lamb_dicke_parameters,
create_geometric_phase_metrics,
get_pwc_acquired_displacements,
get_pwc_acquired_phase,
get_ms_infidelity,
)

#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 = tf.reverse(moduli, [0])
phases_diff_refl = tf.concat([[0], tf.reverse(phases_diff, [0])], 0)
phases_extra = phases[-1] + tf.cumsum(phases_diff_refl)

else:
moduli_refl = tf.reverse(moduli[:-1], [0])
phases_diff_refl = tf.reverse(phases_diff, [0])
phases_extra = phases[-1] + tf.cumsum(phases_diff_refl)

moduli_comb = tf.concat([moduli, moduli_refl], 0)
phases_comb = tf.concat([phases, phases_extra], 0)
return moduli_comb, phases_comb


## Creating Q-CTRL pulses for two-qubit gates

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

#Trap characteristics
number_of_ions = 2
atomic_mass = 171 # Yb ions
com_frequencies = [1.6*1e6, 1.5*1e6, 0.3*1e6]

#Laser characteristics
maximum_rabi_rate = 2 * np.pi * 100e3
laser_wavevector = (2 * np.pi) / 355e-9 * np.sqrt(2) * np.array([1 / np.sqrt(2), 1 / np.sqrt(2), 0])
laser_detuning = 1.6e6 + 4.7e3

#Derived setup parameters
mode_properties = get_default_modes(atomic_mass=atomic_mass,
number_of_ions=number_of_ions,
com_frequencies=com_frequencies)
lamb_dicke_parameters = get_lamb_dicke_parameters(atomic_mass, number_of_ions,
com_frequencies,
laser_wavevector)
relative_detunings = mode_properties['frequencies'] - laser_detuning


We consider Mølmer-Sørensen-type operations, which 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 = e^{i \sum_{j=0}^{N-1} \sum_{k=0}^{j-1} \psi_{jk} \sigma_x^j \sigma_x^k }.$$

Note that multiple entangling operations can be performed in parallel, including multi-qubit entangling operations. We begin with a two-qubit demonstration, and continue in later sections with more exotic operation demonstrations. We configure the target two-qubit 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 we entangle the ion pair (1, 0) and set the relative phase element to $\pi/4$ for maximal two-qubit entanglement.

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


### Generating optimized pulses

We specify the control pulses according to the degrees of freedom in the hardware. In this case, we consider separately-tunable complex drives $\gamma_j (t) = \Omega_j e^{i \phi_j}$ that address each ion $j$ in the trap. Alternatively, the drive could be specified as shared across the ions. We consider piecewise-constant drives with fixed detuning in this example, where the amplitude and phase take constant values for each pulse segment of the given laser drive.

We also specify the cost, which the optimizer attempts to minimize. The cost in the below cell quantifies 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

def create_optimization_specification(variable_factory):
"""This function specifies the control drive variables and the cost function for the optimizer.
It also specifies the optimizer outputs: final cost, drives and infidelity."""

drives = []
for ion in range(number_of_ions):
ion_drive = create_complex_pwc_signal(
#The drive moduli are free variables here. They can also be restricted or fixed.
moduli=variable_factory.create_bounded(
count=number_of_segments, lower_bound=0,
upper_bound=maximum_rabi_rate),
#The drive phases are free variables here. They can also be restricted or fixed.
phases=variable_factory.create_unbounded(
count=number_of_segments, initial_lower_bound=0,
initial_upper_bound=2*np.pi),
duration=duration)
drives.append(ion_drive)

ms_cost = create_geometric_phase_metrics(drives, target, lamb_dicke_parameters, relative_detunings)

#The returned object encapsulates the information required for optimization.
return OptimizationSpecification(
cost=ms_cost['cost'],
output_callable=lambda session: {
'cost': session.run(ms_cost['cost']),
'drive_pulses': [[{'duration': duration, 'value': value}
for duration, value
in zip(drive.durations,
session.run(drive.values))] for drive in drives],
'infidelity': session.run(ms_cost['infidelity'])
})

start_time = time.time()
number_of_optimizations = 10
result0 = optimize(create_optimization_specification=create_optimization_specification,
number_of_optimizations=number_of_optimizations,
seed=None)['output']
print('Execution time for '+ str(number_of_optimizations)+ ' runs (in sec):', time.time()-start_time)
print('Cost = ', result0['cost'])
print('Infidelity = ', result0['infidelity'])

Execution time for 10 runs (in sec): 2.3925867080688477
Cost =  5.480100056898315e-08
Infidelity =  1.6005774483573987e-09

for ion in range(number_of_ions):
control={f'$\gamma_{ion}$': result0['drive_pulses'][ion]}
plot_controls(plt.figure(), control)


We can display the optimized pulse moduli $\Omega_j (t)$ and phases $\phi_j (t)$ for any ion $j$. The above figure displays the dynamics of these quantities for both ions.

It is important to note that the user-specified form of the control pulses can produce control frequencies that bring undesired transitions into resonance. In particular, the lasers are detuned from the carrier transition between the ground and excited states of the ions, however high-frequency components of the control spectrum can induce 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 very small 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 Flexible Optimization User guide.

### Generating robust optimized pulses

Our (pre-release) features include robustness to mode or laser frequency noise, and pulse timing errors. Here we impose symmetry in each ion's drive as described in (Milne et al., Phys. Rev. Applied, 2020) and optimize the drives such that the centre of mass of each mode's trajectory is at zero. This is realized by changing the cost function: here we include cost terms that quantify the differences between the integrated trajectories and zero, as well as the differences between the realized and target phases.

def create_robust_optimization_specification(variable_factory):
"""This function specifies the control drive variables and the cost function for the optimizer.
Here the cost includes robustness terms."""

drives = []
for ion in range(number_of_ions):
#Specification of free variables and combination with reflected signal
number_of_free_segments = int(np.ceil(number_of_segments / 2))
moduli = variable_factory.create_bounded(
count=number_of_free_segments, lower_bound=0,
upper_bound=maximum_rabi_rate)
phases = variable_factory.create_unbounded(
count=number_of_free_segments, initial_lower_bound=0,
initial_upper_bound=2 * np.pi)
full_moduli, full_phases = reflect_signal(number_of_segments, phases, moduli)

ion_drive = create_complex_pwc_signal(moduli=full_moduli,
phases=full_phases,
duration=duration)
drives.append(ion_drive)

ms_cost = create_geometric_phase_metrics(drives, target, lamb_dicke_parameters, relative_detunings,
robust=True)

return OptimizationSpecification(
cost=ms_cost['cost'],
output_callable=lambda session: {
'cost': session.run(ms_cost['cost']),
'drive_pulses': [[{'duration': duration, 'value': value}
for duration, value
in zip(drive.durations,
session.run(drive.values))] for drive in drives],
'infidelity': session.run(ms_cost['infidelity'])
})

#Robust optimization
start_time = time.time()
number_of_optimizations = 10
result = optimize(create_optimization_specification=create_robust_optimization_specification,
number_of_optimizations=number_of_optimizations,
seed=None)['output']
print('Execution time for '+ str(number_of_optimizations)+ ' runs (in sec):', time.time()-start_time)
print('Cost = ', result['cost'])
print('Infidelity = ', result['infidelity'])

Execution time for 10 runs (in sec): 3.8216259479522705
Cost =  2.0694274766979857e-07
Infidelity =  1.9394623684831913e-09

#Drive modulus and phase
for ion in range(number_of_ions):
control={f'$\gamma_{ion}$': result['drive_pulses'][ion]}
plot_controls(plt.figure(), control)


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

### Operation verification with phase space trajectories

We 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. We first examine the (non-robust) standard optimized control, followed by the robust optimized control.

# Standard optimized control

drive_array = [{'values': [segment['value'] for segment in ion_drive], 'durations': [segment['duration'] for segment in ion_drive]}
for ion_drive in result0['drive_pulses']]
drives = [Pwc(values=drive['values'], durations=drive['durations'])
for drive in drive_array]

trajl=[]
sample_times = np.linspace(0, duration, number_of_segments*100)
for mode in range(number_of_ions*3):
trajl.append(get_pwc_acquired_displacements(mode, drives, number_of_ions,
lamb_dicke_parameters, relative_detunings, sample_times))

gs = gridspec.GridSpec(1, 2)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(2 * 5)

for mode in range(number_of_ions):
ax.plot(np.real(trajl[mode]),np.imag(trajl[mode]), label=str(mode % 2))
ax.plot(np.real(trajl[mode][-1]),np.imag(trajl[mode][-1]), marker='x', color='black', markersize=15)
ax.set_xlim(-0.6, 0.6)
ax.set_ylim(-0.6, 0.6)
ax.set_aspect('equal')
ax.set_xlabel('q')
ax.set_ylabel('p')
ax.set_title('$x$-axis', fontsize=14)

for mode in range(number_of_ions, number_of_ions*2):
ax.plot(np.real(trajl[mode]),np.imag(trajl[mode]), label=str(mode % 2))
ax.plot(np.real(trajl[mode][-1]),np.imag(trajl[mode][-1]), marker='x', color='black', markersize=15)
ax.set_xlim(-0.6, 0.6)
ax.set_ylim(-0.6, 0.6)
ax.set_aspect('equal')
ax.set_xlabel('q')
ax.set_ylabel('p')
ax.set_title('$y$-axis', fontsize=14)
ax.legend(title='mode number', bbox_to_anchor=(1.4, 0.8))

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.

We again visualize the phase space trajectories, this time for the robust optimized control.

drive_array_rob = [{'values': [segment['value'] for segment in ion_drive], 'durations': [segment['duration'] for segment in ion_drive]}
for ion_drive in result['drive_pulses']]
drives_rob = [Pwc(values=drive['values'], durations=drive['durations'])
for drive in drive_array_rob]

trajl_r=[]
for mode in range(number_of_ions*3):
trajl_r.append(get_pwc_acquired_displacements(mode, drives_rob, number_of_ions,
lamb_dicke_parameters, relative_detunings, sample_times))

gs = gridspec.GridSpec(1, 2)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(2 * 5)

for mode in range(number_of_ions):
ax.plot(np.real(trajl_r[mode]),np.imag(trajl_r[mode]), label=str(mode % 2))
ax.plot(np.real(trajl_r[mode][-1]),np.imag(trajl_r[mode][-1]), marker='x', color='black', markersize=15)
ax.set_xlim(-0.6, 0.6)
ax.set_ylim(-0.6, 0.6)
ax.set_aspect('equal')
ax.set_xlabel('q')
ax.set_ylabel('p')
ax.set_title('$x$-axis', fontsize=14)

for mode in range(number_of_ions, number_of_ions*2):
ax.plot(np.real(trajl_r[mode]),np.imag(trajl_r[mode]), label=str(mode % 2))
ax.plot(np.real(trajl_r[mode][-1]),np.imag(trajl_r[mode][-1]), marker='x', color='black', markersize=15)
ax.set_xlim(-0.6, 0.6)
ax.set_ylim(-0.6, 0.6)
ax.set_aspect('equal')
ax.set_xlabel('q')
ax.set_ylabel('p')
ax.set_title('$y$-axis', fontsize=14)
ax.legend(title='mode number', bbox_to_anchor=(1.4, 0.8))

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.

### Operation verification with relative phase dynamics

We now show the phase accumulation during the gate operation for the pair of ions. The target phase should be achieved at the end of a successful operation. We begin with the standard optimized control.

phases = []
for ion1 in range(number_of_ions):
for ion2 in range(ion1):
phases.append([])
for end in range(number_of_segments+1):
phases[-1].append(get_pwc_acquired_phase(ion1, ion2, drives, number_of_ions, lamb_dicke_parameters,
relative_detunings, end_segment = end))

colors = sns.color_palette()
fig = plt.figure(figsize = (8, 3))
ion1 = 1
ion2 = 0
plt.plot(np.array(phases[0]), color=colors[0], label=str(ion1)+', '+str(ion2))
plt.yticks([0, np.pi/4],[r"0", r"$\frac{\pi}{4}$"], fontsize=15)
plt.plot([np.pi/4]*67, 'k--', linewidth=0.8)
plt.plot([64,64],[-1,1],'k', alpha=0.5)
plt.ylim((-0.01,0.85))
plt.xlabel('Drive segment')
plt.ylabel('Relative phase')
plt.legend(title='Ion pair', bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
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 gray vertical line). The target phase is marked by the horizontal dashed black line.

We also show the phase accumulation between the ion pair for the robust optimized control.

phases_r = []
for ion1 in range(number_of_ions):
for ion2 in range(ion1):
phases_r.append([])
for end in range(number_of_segments+1):
phases_r[-1].append(get_pwc_acquired_phase(ion1, ion2, drives_rob, number_of_ions, lamb_dicke_parameters,
relative_detunings, end_segment = end))

colors = sns.color_palette()
fig = plt.figure(figsize = (8, 3))
ion1 = 1
ion2 = 0
plt.plot(np.array(phases_r[0]), color=colors[0], label=str(ion1)+', '+str(ion2))
plt.yticks([0, np.pi/4],[r"0", r"$\frac{\pi}{4}$"], fontsize=15)
plt.plot([np.pi/4]*67, 'k--', linewidth=0.8)
plt.plot([64,64],[-1,1],'k', alpha=0.5)
plt.ylim((-0.01,0.85))
plt.xlabel('Drive segment')
plt.ylabel('Relative phase')
plt.legend(title='Ion pair', bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.show()


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

### Pulse timing robustness verification with a quasi-static scan

We assess the robustness of the 'Robust' optimized pulses in comparison to our 'Standard' optimized pulses over 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.

infids=[]
infids_timing_rob=[]
maxshift = 0.02
for shift in np.linspace(-maxshift, maxshift, 21):
drives_var = [Pwc(values=drive['values'], durations=np.array(drive['durations'])*(1+shift))
for drive in drive_array]
infids.append(get_ms_infidelity(drives_var, target, lamb_dicke_parameters, relative_detunings))
drives_rob_var = [Pwc(values=drive['values'], durations=np.array(drive['durations'])*(1+shift))
for drive in drive_array_rob]
infids_timing_rob.append(get_ms_infidelity(drives_rob_var, target, lamb_dicke_parameters, relative_detunings))

plt.clf()
colors = ['#3273DC','#680CEA']
timeshifts=np.linspace(1-maxshift, 1+maxshift, 21)
plt.plot(timeshifts, infids, label='Standard', color=colors[0])
plt.plot(timeshifts, infids_timing_rob, label='Robust', color=colors[1])
plt.xlabel("Pulse timing scaling factor")
plt.ylabel("Infidelity")
plt.legend()
plt.show()


The broader high-fidelity region indicates the benefit of the robust optimized pulses when there is quasi-static noise in the control pulse timings.

### Dephasing robustness verification with a quasi-static scan

We now assess the robustness of the 'Robust' optimized pulses in comparison to our 'Standard' optimized pulses for a systematic scan of the relative detunings (which corresponds to shifting the laser detuning).

infids=[]
infids_rob=[]
minreldet = min(np.abs(relative_detunings))
maxshift = 4000*minreldet/10000
for shift in np.linspace(-maxshift, maxshift, 51):
relative_detuning_var = relative_detunings + shift
infids.append(get_ms_infidelity(drives, target, lamb_dicke_parameters, relative_detuning_var))
infids_rob.append(get_ms_infidelity(drives_rob, target, lamb_dicke_parameters, relative_detuning_var))

plt.clf()
colors = ['#3273DC','#680CEA']
detshifts=np.linspace(-maxshift, maxshift, 51)
plt.plot(detshifts, infids, label='Standard', color=colors[0])
plt.plot(detshifts, infids_rob, label='Robust', color=colors[1])
plt.xlabel("Laser detuning shift $\Delta \delta$ (Hz)")
plt.ylabel("Infidelity")
plt.xlim((-2000, 2000))
plt.ylim((0, 1))
plt.legend()
plt.show()


The broader high-fidelity region indicates the benefit of the robust optimized pulses when faced with dephasing noise.

## Experimental validation of two-ion gates

The optimized 'Standard' control method and the 'Robust' method (using the same robustness conditions) have been performed in experiment (Milne et al., Phys. Rev. Applied, 2020). This work examined phase-modulation for two-qubit gates, where the shared drive has fixed Rabi frequency and variable phase. The target gate operation $U_c$ acts on the initial state $|00\rangle$ to produce a maximally-entangled Bell state: $$U_c |00\rangle = \frac{1}{\sqrt{2}} (|00\rangle -i |11\rangle).$$ The Standard and Robust drives were optimized to perform the operation $U_c$, and were assessed for robustness to laser detuning error both theoretically and in experiment:

Image(filename='Two_qubit_control_legend.png')

Image(filename='Two_qubit_control.png')


Solid lines are from theoretical simulation, and the data points are from the experiment. Experimental data is shown with error bars derived from quantum projection noise on the state population estimates and a fit of the parity contrast. Phase space trajectories are presented in the paper (Milne et al., Phys. Rev. Applied, 2020) for the detuning error marked by the dashed black line.

The figure displays $P_1$, which is the probability of (only) a single ion (either of the two trapped ions) being projected into the state $|1\rangle$. $P_1 > 0$ is thus an indication of error in the gate. Both gates successfully perform the operation with no detuning error. The robust numerical solution has a much broader $P_1 \simeq 0$ region, validated in experiment.

The broad robustness region for the two-ion robust solution is in agreement with the theoretical robustness assessment for the two-qubit gates obtained in the above section. Note that the experimental implementation uses phase modulation with a shared drive, while the above solutions have separate ion-drives each with phase and amplitude modulation. Additionally, the $P_1$ measure is not the same as the operational infidelity plotted in section 1.

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

Preview features in BOULDER OPAL include configurable multi-qubit gates. These gates can also be performed simultaneously (in parallel). We now present simultaneous two-qubit and three-qubit entangling operations to highlight this capability.

As above, we specify the ion, trap and laser characteristics.

#Trap characteristics
number_of_ions = 5
atomic_mass = 171 # Yb ions
com_frequencies = [1.6*1e6, 1.5*1e6, 0.3*1e6]

#Laser characteristics
maximum_rabi_rate = 2 * np.pi * 1000e3
laser_wavevector = (2 * np.pi) / 355e-9 * np.sqrt(2) * np.array([1 / np.sqrt(2), 1 / np.sqrt(2), 0])
laser_detuning = 1.6e6 + 4.7e3

#Derived setup parameters
mode_properties = get_default_modes(atomic_mass=atomic_mass,
number_of_ions=number_of_ions,
com_frequencies=com_frequencies)
lamb_dicke_parameters = get_lamb_dicke_parameters(atomic_mass, number_of_ions,
com_frequencies,
laser_wavevector)
relative_detunings = mode_properties['frequencies'] - laser_detuning


To demonstrate simultaneous two- and three-qubit gates, we configure the target operation by specifiying target relative phases between each ion pair. Recall that element ($j$, $k$) of the target phase array is the relative phase between ions $j$ and $k$, with $k<j$. In the cell below we entangle the ion pair (1, 0) for the two-qubit gate, and pairs (3, 2), (4, 2) and (4, 3) for the three-qubit gate. We set the relative phase element to $\pi/4$ for maximal two-qubit entanglement.

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


### Generating optimized pulses

Here we 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. As above, the optimized cost and infidelity quantify the solution performance.

number_of_segments = 128

def create_optimization_specification(variable_factory):
"""This function specifies the control drive variables and the cost function for the optimizer.
It also specifies the optimizer outputs: final cost, drives and infidelity."""

drives = []
for ion in range(number_of_ions):
ion_drive = create_complex_pwc_signal(
#The drive moduli are free variables here. They can also be restricted or fixed.
moduli=variable_factory.create_bounded(
count=number_of_segments, lower_bound=0,
upper_bound=maximum_rabi_rate),
#The drive phases are free variables here. They can also be restricted or fixed.
phases=variable_factory.create_unbounded(
count=number_of_segments, initial_lower_bound=0,
initial_upper_bound=2*np.pi),
duration=duration)
drives.append(ion_drive)

ms_cost = create_geometric_phase_metrics(drives, target, lamb_dicke_parameters, relative_detunings)

#The returned object encapsulates the information required for optimization.
return OptimizationSpecification(
cost=ms_cost['cost'],
output_callable=lambda session: {
'cost': session.run(ms_cost['cost']),
'drive_pulses': [[{'duration': duration, 'value': value}
for duration, value
in zip(drive.durations,
session.run(drive.values))] for drive in drives],
'infidelity': session.run(ms_cost['infidelity'])
})

start_time = time.time()
number_of_optimizations = 20
result0 = optimize(create_optimization_specification=create_optimization_specification,
number_of_optimizations=number_of_optimizations,
seed=None)['output']
print('Execution time for '+ str(number_of_optimizations)+ ' runs (in sec):', time.time()-start_time)
print('Cost = ', result0['cost'])
print('Infidelity = ', result0['infidelity'])

Execution time for 20 runs (in sec): 32.392768144607544
Cost =  4.734865494953567e-08
Infidelity =  5.460138119417479e-09

ion = 0
control={f'$\gamma_{ion}$': result0['drive_pulses'][ion]}
plot_controls(plt.figure(), control)


We display the optimized pulse modulus and phase dynamics for ion 0.

### Generating robust optimized pulses

As for the two-qubit case, we impose symmetry on each drive and optimize such that the centre of mass of each mode's trajectory is at zero.

def create_robust_optimization_specification(variable_factory):
"""This function specifies the control drive variables and the cost function for the optimizer.
Here the cost includes robustness terms."""

drives = []
for ion in range(number_of_ions):
#Specification of free variables and combination with reflected signal
number_of_free_segments = int(np.ceil(number_of_segments / 2))
moduli = variable_factory.create_bounded(
count=number_of_free_segments, lower_bound=0,
upper_bound=maximum_rabi_rate)
phases = variable_factory.create_unbounded(
count=number_of_free_segments, initial_lower_bound=0,
initial_upper_bound=2 * np.pi)
full_moduli, full_phases = reflect_signal(number_of_segments, phases, moduli)

ion_drive = create_complex_pwc_signal(moduli=full_moduli,
phases=full_phases,
duration=duration)
drives.append(ion_drive)

ms_cost = create_geometric_phase_metrics(drives, target, lamb_dicke_parameters, relative_detunings,
robust=True)

return OptimizationSpecification(
cost=ms_cost['cost'],
output_callable=lambda session: {
'cost': session.run(ms_cost['cost']),
'drive_pulses': [[{'duration': duration, 'value': value}
for duration, value
in zip(drive.durations,
session.run(drive.values))] for drive in drives],
'infidelity': session.run(ms_cost['infidelity'])
})

#Robust optimization
start_time = time.time()
number_of_optimizations = 20
result = optimize(create_optimization_specification=create_robust_optimization_specification,
number_of_optimizations=number_of_optimizations,
seed=None)['output']
print('Execution time for '+ str(number_of_optimizations)+ ' runs (in sec):', time.time()-start_time)
print('Cost = ', result['cost'])
print('Infidelity = ', result['infidelity'])

Execution time for 20 runs (in sec): 104.5357608795166
Cost =  1.9049912894936392e-06
Infidelity =  3.838251994014996e-08

#Drive modulus and phase
ion = 0
control={f'$\gamma_{ion}$': result['drive_pulses'][ion]}
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.

### Operation verification with phase space trajectories

Here we 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. We first examine the (non-robust) standard optimized control, followed by the robust optimized control.

drive_array = [{'values': [segment['value'] for segment in ion_drive], 'durations': [segment['duration'] for segment in ion_drive]}
for ion_drive in result0['drive_pulses']]
drives = [Pwc(values=drive['values'], durations=drive['durations'])
for drive in drive_array]

trajl=[]
sample_times = np.linspace(0, duration, number_of_segments*100)
for mode in range(number_of_ions*3):
trajl.append(get_pwc_acquired_displacements(mode, drives, number_of_ions,
lamb_dicke_parameters, relative_detunings, sample_times))

gs = gridspec.GridSpec(1, 2)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(2 * 5)

for mode in range(number_of_ions):
ax.plot(np.real(trajl[mode]),np.imag(trajl[mode]), label=str(mode % 5))
ax.plot(np.real(trajl[mode][-1]),np.imag(trajl[mode][-1]), marker='x', color='black', markersize=15)
ax.set_xlim(-3, 3)
ax.set_ylim(-3, 3)
ax.set_aspect('equal')
ax.set_xlabel('q')
ax.set_ylabel('p')
ax.set_title('$x$-axis', fontsize=14)

for mode in range(number_of_ions, number_of_ions*2):
ax.plot(np.real(trajl[mode]),np.imag(trajl[mode]), label=str(mode % 5))
ax.plot(np.real(trajl[mode][-1]),np.imag(trajl[mode][-1]), marker='x', color='black', markersize=15)
ax.set_xlim(-3, 3)
ax.set_ylim(-3, 3)
ax.set_aspect('equal')
ax.set_xlabel('q')
ax.set_ylabel('p')
ax.set_title('$y$-axis', fontsize=14)
ax.legend(title='mode number', bbox_to_anchor=(1.4, 0.8))

plt.show()


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.

We again visualize the phase space trajectories, this time for the robust optimized control.

drive_array_rob = [{'values': [segment['value'] for segment in ion_drive], 'durations': [segment['duration'] for segment in ion_drive]}
for ion_drive in result['drive_pulses']]
drives_rob = [Pwc(values=drive['values'], durations=drive['durations'])
for drive in drive_array_rob]

trajl_r=[]
for mode in range(number_of_ions*3):
trajl_r.append(get_pwc_acquired_displacements(mode, drives_rob, number_of_ions,
lamb_dicke_parameters, relative_detunings, sample_times))

gs = gridspec.GridSpec(1, 2)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(2 * 5)

for mode in range(number_of_ions):
ax.plot(np.real(trajl_r[mode]),np.imag(trajl_r[mode]), label=str(mode % 5))
ax.plot(np.real(trajl_r[mode][-1]),np.imag(trajl_r[mode][-1]), marker='x', color='black', markersize=15)
ax.set_xlim(-3, 3)
ax.set_ylim(-3, 3)
ax.set_aspect('equal')
ax.set_xlabel('q')
ax.set_ylabel('p')
ax.set_title('$x$-axis', fontsize=14)

for mode in range(number_of_ions, number_of_ions*2):
ax.plot(np.real(trajl_r[mode]),np.imag(trajl_r[mode]), label=str(mode % 5))
ax.plot(np.real(trajl_r[mode][-1]),np.imag(trajl_r[mode][-1]), marker='x', color='black', markersize=15)
ax.set_xlim(-3, 3)
ax.set_ylim(-3, 3)
ax.set_aspect('equal')
ax.set_xlabel('q')
ax.set_ylabel('p')
ax.set_title('$y$-axis', fontsize=14)
ax.legend(title='mode number', bbox_to_anchor=(1.4, 0.8))

plt.show()


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

### Operation verification with relative phase dynamics

We now show the phase accumulation for each pair of ions. The target phases for each pair of ions should be achieved at the end of a successful operation. We begin with the standard optimized control.

phases = []
for ion1 in range(number_of_ions):
for ion2 in range(ion1):
phases.append([])
for end in range(number_of_segments+1):
phases[-1].append(get_pwc_acquired_phase(ion1, ion2, drives, number_of_ions, lamb_dicke_parameters,
relative_detunings, end_segment = end))

colors = sns.color_palette()
gs = gridspec.GridSpec(2, 1)
fig = plt.figure()
fig.set_figheight(2 * 3)
fig.set_figwidth(8)
ind=0
for ion1 in range(number_of_ions):
for ion2 in range(ion1):
if target[ion1][ion2] == np.pi/4:
ax0.plot(np.array(phases[ind]), color=colors[ind], label=str(ion1)+', '+str(ion2))
else:
ax1.plot(np.array(phases[ind]), color=colors[ind], label=str(ion1)+', '+str(ion2))
ind += 1

unit = np.pi/4
ymax = 4*np.pi/4
ymin = -4*np.pi/4
y_tick = np.arange(ymin, ymax+unit, unit)
y_label = [r"$-\pi$", r"$-3\pi/4$", r"$-\pi/2$", r"$-\pi/4$", r"0"]
y_label += [r"$\pi/4$", r"$\pi/2$", r"$3\pi/4$", r"$\pi$"]
ax0.set_yticks(y_tick)
ax0.set_yticklabels(y_label, fontsize=12)
ax1.set_yticks(y_tick)
ax1.set_yticklabels(y_label, fontsize=12)
ax0.plot([np.pi/4]*130, 'k--', linewidth=0.8)
ax1.plot([0]*130, 'k--', linewidth=0.8)
ax0.plot([128,128],[-4,4],'k', alpha=0.5)
ax1.plot([128,128],[-4,4],'k', alpha=0.5)
ax0.set_ylim((-3.3, 3.15))
ax1.set_ylim((-3.3, 3.15))
ax0.set(xticks=[])
ax1.set_xlabel('Drive segment')
ax0.set_ylabel('Relative phase')
ax1.set_ylabel('Relative phase')
ax0.set_title('Pairs with target phase $\psi = \pi/4$')
ax1.set_title('Pairs with target phase $\psi = 0$')
ax0.legend(title='Ion pair', bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
ax1.legend(title='Ion pair', bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.show()


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

We also show the phase accumulation between each pair of ions for the robust optimized control.

phases_r = []
for ion1 in range(number_of_ions):
for ion2 in range(ion1):
phases_r.append([])
for end in range(number_of_segments+1):
phases_r[-1].append(get_pwc_acquired_phase(ion1, ion2, drives_rob, number_of_ions, lamb_dicke_parameters,
relative_detunings, end_segment = end))

colors = sns.color_palette()
gs = gridspec.GridSpec(2, 1)
fig = plt.figure()
fig.set_figheight(2 * 3)
fig.set_figwidth(8)
ind=0
for ion1 in range(number_of_ions):
for ion2 in range(ion1):
if target[ion1][ion2] == np.pi/4:
ax0.plot(np.array(phases_r[ind]), color=colors[ind], label=str(ion1)+', '+str(ion2))
else:
ax1.plot(np.array(phases_r[ind]), color=colors[ind], label=str(ion1)+', '+str(ion2))
ind += 1

unit = np.pi/4
ymax = 3*np.pi/4
ymin = -3*np.pi/4
y_tick = np.arange(ymin, ymax+unit, unit)
y_label = [r"$-3\pi/4$", r"$-\pi/2$", r"$-\pi/4$", r"0"]
y_label += [r"$\pi/4$", r"$\pi/2$", r"$3\pi/4$"]
ax0.set_yticks(y_tick)
ax0.set_yticklabels(y_label, fontsize=12)
ax1.set_yticks(y_tick)
ax1.set_yticklabels(y_label, fontsize=12)
ax0.plot([np.pi/4]*130, 'k--', linewidth=0.8)
ax1.plot([0]*130, 'k--', linewidth=0.8)
ax0.plot([128,128],[-4,4],'k', alpha=0.5)
ax1.plot([128,128],[-4,4],'k', alpha=0.5)
ax0.set_ylim((-2.5,2.5))
ax1.set_ylim((-2.5,2.5))
ax0.set(xticks=[])
ax1.set_xlabel('Drive segment')
ax0.set_ylabel('Relative phase')
ax1.set_ylabel('Relative phase')
ax0.set_title('Pairs with target phase $\psi = \pi/4$')
ax1.set_title('Pairs with target phase $\psi = 0$')
ax0.legend(title='Ion pair', bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
ax1.legend(title='Ion pair', bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.show()


As above, these figures display the relative phase dynamics for the duration of the gate operation, and the end of the operation is marked by vertical gray lines. For clarity, we display different plots for ion pairs with different target relative phases. 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 every ion is addressed by our chosen drives. 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.

### Pulse timing robustness verification with a quasi-static scan

We assess the robustness of the 'Robust' optimized pulses in comparison to our 'Standard' optimized pulses over 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.

infids=[]
infids_timing_rob=[]
maxshift = 0.004
for shift in np.linspace(-maxshift, maxshift, 21):
drives_var = [Pwc(values=drive['values'], durations=np.array(drive['durations'])*(1+shift))
for drive in drive_array]
infids.append(get_ms_infidelity(drives_var, target, lamb_dicke_parameters, relative_detunings))
drives_rob_var = [Pwc(values=drive['values'], durations=np.array(drive['durations'])*(1+shift))
for drive in drive_array_rob]
infids_timing_rob.append(get_ms_infidelity(drives_rob_var, target, lamb_dicke_parameters, relative_detunings))

plt.clf()
colors = ['#3273DC','#680CEA']
timeshifts=np.linspace(1-maxshift, 1+maxshift, 21)
plt.plot(timeshifts, infids, label='Standard', color=colors[0])
plt.plot(timeshifts, infids_timing_rob, label='Robust', color=colors[1])
plt.xlabel("Pulse timing scaling factor")
plt.ylabel("Infidelity")
plt.legend()
plt.show()


The broader high-fidelity region indicates the benefit of the robust optimized pulses when there is quasi-static noise in the control pulse timings.

### Dephasing robustness verification with a quasi-static scan

We next assess the robustness of the 'Robust' optimized pulses in comparison to our 'Standard' optimized pulses for a systematic scan of the relative detunings (which corresponds to shifting the laser detuning).

infids=[]
infids_rob=[]
minreldet = min(np.abs(relative_detunings))
maxshift = 4*minreldet/10000 * 100
for shift in np.linspace(-maxshift, maxshift, 31):
relative_detuning_var = relative_detunings + shift
infids.append(get_ms_infidelity(drives, target, lamb_dicke_parameters, relative_detuning_var))
infids_rob.append(get_ms_infidelity(drives_rob, target, lamb_dicke_parameters, relative_detuning_var))

plt.clf()
colors = ['#3273DC','#680CEA']
detshifts=np.linspace(-maxshift, maxshift, 31)
plt.plot(detshifts, infids, label='Standard', color=colors[0])
plt.plot(detshifts, infids_rob, label='Robust', color=colors[1])
plt.xlabel("Laser detuning shift $\Delta \delta$ (Hz)")
plt.ylabel("Infidelity")
plt.xlim((-200,200))
plt.ylim((1e-9,1e0))
plt.legend()
plt.show()


The broader high-fidelity region indicates the benefit of the robust optimized pulses when faced with dephasing noise.

The additional robustness for this more general gate is in agreement with the experimental result in section 2 for two-ion gates, using related robustness methodology. These configurable gates are yet to be implemented in experiment: please contact us if you want a live demonstration of these features for your hardware.

## Creating Q-CTRL pulses for 4-of-8 qubit gates

Preview features in BOULDER OPAL include configurable multi-qubit gates. We now demonstrate a 4-of-8 qubit gate solution with specified control constraints: the drive is shared between the four target ions, with a maximum modulus slew rate. We set the target relative phase for each pair within the four ions to $\pi/4$.

#Trap characteristics
number_of_ions = 8
atomic_mass = 171 # Yb ions
com_frequencies = [1.6*1e6, 1.5*1e6, 0.3*1e6]

#Laser characteristics
maximum_rabi_rate = 2 * np.pi * 1000e3
maximum_modulus_slew_rate = 2 * np.pi * 100e3
laser_wavevector = (2 * np.pi) / 355e-9 * np.sqrt(2) * np.array([1 / np.sqrt(2), 1 / np.sqrt(2), 0])
laser_detuning = 1.6e6 + 4.7e3

#Derived setup parameters
mode_properties = get_default_modes(atomic_mass=atomic_mass,
number_of_ions=number_of_ions,
com_frequencies=com_frequencies)
lamb_dicke_parameters = get_lamb_dicke_parameters(atomic_mass, number_of_ions,
com_frequencies,
laser_wavevector)
relative_detunings = mode_properties['frequencies'] - laser_detuning

#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.zeros((8,8))
for ion1 in range(4):
for ion2 in range(ion1):
target[ion1][ion2] = np.pi/4


### Generating optimized pulses

Here we consider a complex drive $\gamma (t) = \Omega e^{i \phi}$ that addresses the four target ions in the chain of eight ions. The drives are piecewise-constant with amplitude and phase modulation. Here, the amplitude modulation is constrained by the maximum modulus slew rate defined above in the laser characteristics.

number_of_segments = 128

def create_optimization_specification(variable_factory):
"""This function specifies the control drive variables and the cost function for the optimizer.
It also specifies the optimizer outputs: final cost, drives and infidelity."""

ion_drive = create_complex_pwc_signal(
moduli=create_anchored_difference_bounded_variables(
variable_factory=variable_factory,
count=number_of_segments,
lower_bound=0.,
upper_bound=maximum_rabi_rate,
difference_bound=maximum_modulus_slew_rate),
phases=variable_factory.create_unbounded(
count=number_of_segments, initial_lower_bound=0,
initial_upper_bound=2*np.pi),
duration=duration)
ion_nodrive = create_complex_pwc_signal(
moduli=tf.constant(np.zeros(1)),
phases=tf.constant(np.zeros(1)),
duration=duration)
drives = [ion_drive]*4 + [ion_nodrive]*4

ms_cost = create_geometric_phase_metrics(drives, target, lamb_dicke_parameters, relative_detunings)

#The returned object encapsulates the information required for optimization.
return OptimizationSpecification(
cost=ms_cost['cost'],
output_callable=lambda session: {
'cost': session.run(ms_cost['cost']),
'drive_pulses': [[{'duration': duration, 'value': value}
for duration, value
in zip(drive.durations,
session.run(drive.values))] for drive in drives],
'infidelity': session.run(ms_cost['infidelity'])
})

start_time = time.time()
number_of_optimizations = 20
result0 = optimize(create_optimization_specification=create_optimization_specification,
number_of_optimizations=number_of_optimizations,
seed=None)['output']
print('Execution time for '+ str(number_of_optimizations)+ ' runs (in sec):', time.time()-start_time)
print('Cost = ', result0['cost'])
print('Infidelity = ', result0['infidelity'])

Execution time for 20 runs (in sec): 162.9632818698883
Cost =  1.0212473294706e-07
Infidelity =  5.572289929478558e-08

ion = 0
control={f'$\gamma$': result0['drive_pulses'][ion]}
plot_controls(plt.figure(), control)


The drive is within the user-defined bounds with the given maximum slew rate.

### Operation verification with phase space trajectories

Here we 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.

drive_array = [{'values': [segment['value'] for segment in ion_drive], 'durations': [segment['duration'] for segment in ion_drive]}
for ion_drive in result0['drive_pulses']]
drives = [Pwc(values=drive['values'], durations=drive['durations'])
for drive in drive_array]

trajl=[]
sample_times = np.linspace(0, duration, number_of_segments*100)
for mode in range(number_of_ions*3):
trajl.append(get_pwc_acquired_displacements(mode, drives, number_of_ions,
lamb_dicke_parameters, relative_detunings, sample_times))

gs = gridspec.GridSpec(1, 2)
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(2 * 5)

for mode in range(number_of_ions):
ax.plot(np.real(trajl[mode]),np.imag(trajl[mode]), label=str(mode % 8))
ax.plot(np.real(trajl[mode][-1]),np.imag(trajl[mode][-1]), marker='x', color='black', markersize=15)
ax.set_xlim(-3, 3)
ax.set_ylim(-3, 3)
ax.set_aspect('equal')
ax.set_xlabel('q')
ax.set_ylabel('p')
ax.set_title