# Optimization

### Robust and optimal optimizations of pulses that control a quantum system

The Q-CTRL Python Package allows you to create custom optimized controls to implement high-quality gates tailored to your specific hardware. In this guide we show how to set up and run optimizations for different systems using the Q-CTRL Python Package.

```
# Essential imports
import numpy as np
from qctrl import Qctrl
# Plotting imports
from qctrlvisualizer import plot_controls
import matplotlib.pyplot as plt
# Starting a session with the API
qctrl = Qctrl()
```

## Worked example: robust control of a single qubit

We first present a detailed worked example of robust optimization in a single-qubit system. Specifically, we consider a single-qubit system represented by the following Hamiltonian:

\begin{align*} H(t) &= \frac{\nu}{2} \sigma_{z} + \frac{1}{2}\left(\Omega(t)\sigma_{+} + \text{h.c.}\right) + \frac{\Delta(t)}{2} \sigma_{z} + \eta(t) \sigma_{z} \,, \end{align*}where $\nu$ is the qubit detuning, $\Omega(t)$ is a time-dependent Rabi rate, $\Delta(t)$ is a time-dependent clock shift, $\eta(t)$ is a small, slowly-varying stochastic dephasing noise process and $\sigma_{\{+,z\}}$ are the Pauli matrices.

The functions of time $\Omega(t)$ and $\Delta(t)$ are not predetermined, and instead are optimized by the Q-CTRL optimization engine in order to achieve some target operation.

### Creating the system, controls, and pulses

As described in the Setup feature guide, we first set up Python objects representing the system, controls, and pulses. Note that we use "optimum" pulses, indicating that the pulses are free to be optimized by the Q-CTRL optimization engine (up to certain bounds).

```
single_qubit_system = qctrl.factories.systems.create(
name="Single-qubit with noise",
hilbert_space_dimension=2,
)
# Define the standard matrices
identity = np.array([[1, 0], [0, 1]]) # or np.eye(2)
sigma_z = np.array([[1, 0], [0, -1]]) # or np.diag([1, -1])
sigma_p = np.array([[0, 0], [1, 0]]) # or np.diag([1], -1)
Y = np.array([[0, -1j], [1j, 0]])
# Define the physical constants
nu = 2*np.pi * 0.5 * 1e6 #Hz
# Define the control objects
drift = qctrl.factories.drift_controls.create(
name='Detuning',
system=single_qubit_system,
operator=nu*sigma_z/2,
)
drive = qctrl.factories.drive_controls.create(
name='Rabi rate',
system=single_qubit_system,
operator=sigma_p/2,
)
clock = qctrl.factories.shift_controls.create(
name='Clock',
system=single_qubit_system,
operator=sigma_z/2,
)
# Define the noise objects
dephasing = qctrl.factories.additive_noises.create(
name='Dephasing',
system=single_qubit_system,
operator=sigma_z,
)
# Define physical constants
omega_max = 2*np.pi * 0.5e6 #Hz
delta_max = 2*np.pi * 0.5e6 #Hz
segment_count = 5
duration = 10e-6 #s
# Define pulse objects
drive_pulse = qctrl.factories.optimum_pulses.create(
control=drive,
upper_bound=omega_max,
fixed_modulus=False,
segment_count=segment_count,
duration=duration,
)
shift_pulse = qctrl.factories.optimum_pulses.create(
control=clock,
upper_bound=delta_max,
fixed_modulus=False,
segment_count=segment_count,
duration=duration,
)
```

### Creating the target

We will produce a robust control solution that implements a $Y$ operation on the qubit. The form of the solution will be a set of time-varying pulses that define the modulation to be applied to the system in order to achieve the desired operation $Y$,

\begin{align*} Y = P \left( \mathcal{T}\exp\left(-i\int_0^T H(t) dt \right) \right) P, \end{align*}where the projector $P = |0 \rangle \langle 0| + |1 \rangle \langle 1|$ is the identity. We express this target operation as a `target`

object.

```
single_qubit_target = qctrl.factories.targets.create(
unitary_operator=Y,
projection_operator=identity,
system=single_qubit_system,
)
```

### Running an optimization

With the system and target set up, we may run the optimization (i.e. calculate optimized control pulses) using the `robust_optimization`

service. The resulting optimized system is returned, and includes both the optimized control pulses and a cost indicating the quality of the optimization result. For a high-quality optimization, the cost should be very close to zero.

```
single_qubit_system = qctrl.services.robust_optimization.run(system=single_qubit_system)
```

```
# Check the cost of the system; it should be << 1.
print("Optimized cost:\t", single_qubit_system.cost)
```

### Visualizing the pulses

Before we can visualize the optimized pulses, we have to extract them from the optimized system via the corresponding control objects. Each pulse is represented as a list of dictionaries corresponding to each pulse segment. We may then use the Q-CTRL Python Visualizer package to plot the optimized pulses.

```
# Get the optimized controls
single_qubit_controls = {control.name: control.pulse.segments
for control in single_qubit_system.controls if control.pulse}
```

```
# Plot the optimized controls
plot_controls(plt.figure(), single_qubit_controls)
plt.show()
```

## Example: two-qubit system with leakage

Next we consider a system consisting of two qutrits, modeling a two-qubit system in which the qubits are subject to leakage to a third energy level. The system is represented by the following Hamiltonian:

\begin{align*} H(t) & = \Delta_1(t) n_1 + \frac{\delta_1}{2} n_1\left(n_1 - 1\right) + \frac{1}{2}\left(\Omega_1(t)a_1^\dagger + \text{h.c.}\right) \\ & + \Delta_2(t) n_2 + \frac{\delta_2}{2} n_2\left(n_2 - 1\right) + \frac{1}{2}\left(\Omega_2(t)a_2^\dagger + \text{h.c.}\right) \\ & + g(t)\left(a_1^\dagger a_2 + \text{h.c.}\right)\,, \end{align*}with $n_1=n\otimes I$, $n_2=I\otimes n$, $a_1=a\otimes I$ and $a_2=I\otimes a$, where $I$ is the identity operator of a single-qubit, $n=|1 \rangle \langle 1 | + 2 |2 \rangle \langle 2 |$ is the number operator and $a = |0 \rangle \langle 1 | + \sqrt{2} |1 \rangle \langle 2 |$. The qutrits are treated as oscillators (truncated to three levels) with the same anharmonicity of $\delta_1=\delta_2=\delta$. Each oscillator $k$ is driven with a controllable Rabi rate $\Omega_k(t)$ and tunable frequency shift of $\Delta_k(t)$, while the controllable coupling between them is given by $g(t)$.

We are going to produce an optimal-control solution that implements an iSWAP operation on the two qubits. The third levels provide pathways for leakage errors, so we target a control that ensures a minimization of the leakage errors. More specifically, our aim is to find control solutions that achieve an iSWAP gate in the projected subspace of the computational states, $|00\rangle$, $|01\rangle$, $|10\rangle$ and $|11 \rangle$. The form of the solution will be a set of time-varying pulses that define the modulation to be applied to the system in order to achieve the desired operation $U_{\rm iSWAP}$,

\begin{align*} U_{\rm iSWAP} = P \left( \mathcal{T}\exp\left(-i\int_0^T H(t) dt \right) \right) P, \end{align*}where $P = |00 \rangle \langle 00| + |01 \rangle \langle 01|+|10 \rangle \langle 10|+|11 \rangle \langle 11|$ represents projection onto the computational subspace.

```
# Create the system
two_qubit_system = qctrl.factories.systems.create(
name="Two-qubit with leakage",
hilbert_space_dimension=9,
)
# Define standard matrices
identity = np.eye(3)
a = np.diag([1, np.sqrt(2)], 1)
n = np.diag([0, 1, 2])
qubit_projector = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 0]], dtype=np.complex)
iswap = np.eye(9, dtype=np.complex)
iswap[1,1] = 0
iswap[1,3] = 1j
iswap[3,3] = 0
iswap[3,1] = 1j
# Define physical parameters
delta = 2 * np.pi * -100.0 * 1e6 #Hz
omega_max = 2 * np.pi * 200 * 1e6 #Hz
Delta_max = 2 * np.pi * 50 * 1e6 #Hz
g_max = 2 * np.pi * 10 * 1e6 #Hz
duration = 1e-6 #s
segment_count = 10
# Define control objects
drive_a = qctrl.factories.drive_controls.create(
name='Rabi rate A',
system=two_qubit_system,
operator=np.kron(a, identity)/2,
)
drive_b = qctrl.factories.drive_controls.create(
name='Rabi rate B',
system=two_qubit_system,
operator=np.kron(identity, a)/2,
)
coupling = qctrl.factories.drive_controls.create(
name='Coupling',
system=two_qubit_system,
operator=np.kron(a.T, a),
)
clock_a = qctrl.factories.shift_controls.create(
name='Clock A',
system=two_qubit_system,
operator=np.kron(n, identity),
)
clock_b = qctrl.factories.shift_controls.create(
name='Clock B',
system=two_qubit_system,
operator=np.kron(identity, n),
)
drift = qctrl.factories.drift_controls.create(
name='Anharmonicities',
system=two_qubit_system,
operator=0.5*delta*(
np.kron(n.dot(n) - n, identity) +
np.kron(identity, n.dot(n) - n)),
)
# Define pulse objects
drive_a_pulse = qctrl.factories.optimum_pulses.create(
control=drive_a,
upper_bound=omega_max,
fixed_modulus=False,
segment_count=segment_count,
duration=duration,
)
drive_b_pulse = qctrl.factories.optimum_pulses.create(
control=drive_b,
upper_bound=omega_max,
fixed_modulus=False,
segment_count=segment_count,
duration=duration,
)
coupling_pulse = qctrl.factories.optimum_pulses.create(
control=coupling,
upper_bound=g_max,
fixed_modulus=False,
segment_count=segment_count,
duration=duration,
)
clock_a_pulse = qctrl.factories.optimum_pulses.create(
control=clock_a,
upper_bound=Delta_max,
fixed_modulus=False,
segment_count=segment_count,
duration=duration,
)
clock_b_pulse = qctrl.factories.optimum_pulses.create(
control=clock_b,
upper_bound=Delta_max,
fixed_modulus=False,
segment_count=segment_count,
duration=duration,
)
# Define the target
two_qubit_target = qctrl.factories.targets.create(
unitary_operator=iswap,
projection_operator = np.kron(qubit_projector, qubit_projector),
system=two_qubit_system
)
# Run the optimization
two_qubit_system = qctrl.services.optimal_optimization.run(system=two_qubit_system)
print("Optimized cost:\t", two_qubit_system.cost)
# Get the optimized controls
two_qubit_controls = {control.name: control.pulse.segments
for control in two_qubit_system.controls if control.pulse}
# Plot the optimized controls
plot_controls(plt.figure(), two_qubit_controls)
plt.show()
```

## Example: single-qubit system with dephasing and leakage

In this example, we consider a single-qubit subject to dephasing that also has a third level, whose presence is responsible for leakage errors. Like in the previous section, we will treat this system as an oscillator truncated to three levels, and with anharmonicity $\chi$. The two controls available are the Rabi rate $\Omega(t)$ and the tunable frequency shift $\Delta(t)$, so that the entire system is represented by the Hamiltonian:

$$ H(t) = \frac{\chi}{2} (a^\dagger)^2 a^2 + \Omega(t) a + \Omega^*(t) a^\dagger + \frac{\Delta(t)}{2} a^\dagger a + \eta(t) a^\dagger a, $$where all the operators are the same as defined in the previous sections, i.e. $a = \left| 0 \right\rangle \left\langle 1 \right| + \sqrt{2} \left| 1 \right\rangle \left\langle 2 \right|$. Notice that we are also inserting a term $\eta(t)$ in the Hamiltonian, which represents dephasing errors.

Our objective is to find control solutions to achieve a Hadamard gate inside the subspace of the computational states $\left| 0 \right\rangle$ and $\left| 1 \right\rangle$. We want these controls to be robust against leakage errors and also against dephasing. The form of the solution will be a set of time-varying pulses that define the modulation to be applied to the system in order to achieve the desired Hadamard gate $U_{\rm H}$,

\begin{align*} U_{\rm H} = P \left( \mathcal{T}\exp\left(-i\int_0^T H(t) dt \right) \right) P, \end{align*}where $P = |0 \rangle \langle 0| + |1 \rangle \langle 1|$ represents a projection onto the computational subspace.

In this section, we are also exemplifying the use of `maximum_slew_rate`

, an argument passed to `qctrl.factories.optimum_pulses.create`

to specify the maximum variation allowed between two neighboring control segments.
In the case of a drive, whose controls have complex values, this maximum rate represents the maximum variation of the modulus of the control.
For shifts, whose values are real, the maximum rate represents variation in its total value, not just the modulus.
In the controls below, we impose on the drive a `maximum_slew_rate`

that is three times smaller than the `maximum_slew_rate`

of the shift, so that a comparison of the plot of these two controls illustrates how a more strict slew rate affects the form of the pulses.

```
# Define the physical constants
new_duration = 10. * 1e-9 # s
chi = 2*np.pi * -100. * 1e6 # Hz
new_omega_max = 2 * np.pi * 100. * 1e6 # Hz
new_delta_max = 2 * np.pi * 100. * 1e6 # Hz
max_omega_slew_rate = new_omega_max/15.
max_delta_slew_rate = new_omega_max/5.
# Define optimization parameters
new_segment_count = 100
# Define standard matrices
a2 = np.matmul(a, a)
ada = np.matmul(a.T.conj(), a)
ad2a2 = np.matmul(a2.T.conj(), a2)
hadamard = np.array([[1./np.sqrt(2.), 1./np.sqrt(2.), 0.],
[1./np.sqrt(2.), -1./np.sqrt(2.), 0.],
[0., 0., 1.]], dtype=np.complex)
# Create the system
new_single_qubit_system = qctrl.factories.systems.create(
name="Single Qubit with Leakage",
hilbert_space_dimension=3,
)
# Define control objects
new_drive = qctrl.factories.drive_controls.create(
name="Microwave",
operator=a,
system=new_single_qubit_system,
)
new_shift = qctrl.factories.shift_controls.create(
name="Clock",
operator=ada/2.,
system=new_single_qubit_system,
)
new_drift = qctrl.factories.drift_controls.create(
name="Anharmonic Oscillator",
operator=chi*ad2a2/2.,
system=new_single_qubit_system,
)
# Define the noise objects
dephasing = qctrl.factories.additive_noises.create(
name='Dephasing',
operator=ada,
system=new_single_qubit_system,
)
# Define pulse objects
new_drive_pulse = qctrl.factories.optimum_pulses.create(
control=new_drive,
upper_bound=new_omega_max,
fixed_modulus=False,
segment_count=new_segment_count,
duration=new_duration,
maximum_slew_rate=max_omega_slew_rate,
)
new_shift_pulse = qctrl.factories.optimum_pulses.create(
control=new_shift,
upper_bound=new_delta_max,
fixed_modulus=False,
segment_count=new_segment_count,
duration=new_duration,
maximum_slew_rate=max_delta_slew_rate,
)
# Define the target
new_target = qctrl.factories.targets.create(
unitary_operator=hadamard,
projection_operator=qubit_projector,
system=new_single_qubit_system,
)
# Run the optimization
new_single_qubit_system = qctrl.services.robust_optimization.run(
system=new_single_qubit_system)
print("Optimized cost:\t", new_single_qubit_system.cost)
# Get the optimized controls
new_single_qubit_controls = {control.name: control.pulse.segments
for control in new_single_qubit_system.controls if control.pulse}
# Plot the optimized controls
plot_controls(plt.figure(), new_single_qubit_controls)
plt.show()
```