# Performing optimal Fock state generation in superconducting resonators

**Engineering fast cavity state generation in superconducting cavity-qubit systems**

BOULDER OPAL provides a flexible toolkit for control optimization in high-dimensional and complex Hilbert spaces. Superconducting cavity-qubit systems exhibit a variety of interaction terms and rates, such that analytical schemes typically isolate target transitions or regimes to simplify the dynamics. This complexity, however, can be a resource for optimized control schemes that can exploit the interplay of interactions and coupling rates to perform faster or higher-fidelity operations.

In this Application note, we demonstrate the flexibility of BOULDER OPAL optimization tools to directly generate and validate Fock state preparation schemes that fully leverage your available controls, without inducing leakage to higher cavity or transmon levels. This direct Fock state generation is complementary to the approach of optimizing SNAP gates demonstrated in our Application note. Here, we will cover Fock state generation with different control freedoms:

- Optimizing both a transmon and a cavity drive in the dispersive regime,
- Optimizing a transmon drive and tunable transmon frequency to exploit resonant qubit-cavity interactions,
- Validating the optimized scheme performance using simulation.

Ultimately, we demonstrate how model-based numeric optimization provides a means to dramatically enhance gate speed without performance degradation in qubit-oscillator systems.

```
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np
import time
from qctrlvisualizer import get_qctrl_style, plot_controls
# Q-CTRL imports
from qctrl import Qctrl
# Starting a session with the API
qctrl = Qctrl()
# Plotting style
plt.style.use(get_qctrl_style())
```

## Generate a Fock state in the dispersive regime

The system consists of a superconducting transmon system coupled to a cavity in the dispersive limit, as given by the following Hamiltonian:

$$H_0 = \omega_C a^\dagger a + \frac{K}{2} (a^\dagger)^2 a^2 + \omega_T b^\dagger b + \frac{\alpha}{2} (b^\dagger)^2 b^2 + \chi a^\dagger a b^\dagger b, $$where $\omega_C$ is the cavity transition frequency, $a$ is the annihilation operator of a cavity excitation, $b$ is the annihilation operator of the transmon system, $K$ is the Kerr coefficient, $\omega_T$ is the transmon frequency, $\alpha$ is the transmon nonlinearity, and $\chi$ is the dispersive shift. Consider drives applied to both the transmon and the cavity with frequencies $\nu_T$ and $\omega_C$, respectively. Then the Hamiltonian in the interaction picture with respect to $\nu_T b^\dagger b + \omega_C a^\dagger a$ becomes:

$$H = \delta b^\dagger b + \frac{K}{2} (a^\dagger)^2 a^2 + \frac{\alpha}{2} (b^\dagger)^2 b^2 + \chi a^\dagger a b^\dagger b + \left(\gamma_T (t) b + H.c.\right) + \left(\gamma_C (t) a + H.c.\right), $$where $\delta = \omega_T - \nu_T$, and $\gamma_{T(C)}(t)= I_{T(C)}(t) + i Q_{T(C)}(t)$ is the complex drive amplitude on the transmon qubit (cavity). The basis states in the Hilbert space will be denoted by $|i,j\rangle =|i\rangle_T \otimes |j\rangle_C$, for the transmon number state $|i\rangle_T$ and cavity number state $|j\rangle_C$.

### Creating optimized controls

Using the Q-CTRL Python package, we apply the optimizer to obtain smooth controls that produce the target state $|0, 2\rangle$, starting from the initial ground state $|0, 0\rangle$.

First, we set the parameters for the system, the optimizer, and the target operation. Optimization variables include the truncation dimensions for the optimizer, the target `gate_duration`

, and maximum drive strengths. To create a smooth pulse, as described in this user guide, we also define a cutoff frequency for the sinc filter that smooths the drives before they are resampled with the specified `number_of_segments`

. These parameters, together with the operators determining the system dynamics, are sufficient to characterize the optimization procedure.

```
# System parameters
chi = -2 * np.pi * 2194 * 1e3
K = -2 * np.pi * 3.7 * 1e3
delta = -2.5e6
alpha = -2 * np.pi * 236 * 1e6
# Optimizer parameters
dimt = 2 # Transmon dimensions
dimc = 12 # Cavity dimensions
gate_duration = 550e-9 # s
max_rabi_rate_T = 2 * np.pi * 8e6 # Maximum transmon drive
max_rabi_rate_C = 2 * np.pi * 8e6 # Maximum cavity drive
number_of_optimizer_vars = 64
number_of_segments = 128
sinc_cutoff_frequency_drives = 2 * np.pi * 4e7
# Annihilation and creation operators for the transmon and cavity
a = np.kron(np.eye(dimt), np.diag(np.sqrt(np.arange(1, dimc)), k=1))
b = np.kron(np.diag(np.sqrt(np.arange(1, dimt)), k=1), np.eye(dimc))
# System Hamiltonians without drives
H_sys = (
delta * b.T @ b
+ K / 2 * a.T @ a.T @ a @ a
+ alpha / 2 * b.T @ b.T @ b @ b
+ chi * b.T @ b @ a.T @ a
)
# Define target operation
tr_g = np.zeros(dimt)
tr_g[0] = 1.0
cav_g = np.zeros(dimc)
cav_g[0] = 1.0
cav_3 = np.zeros(dimc)
cav_3[3] = 1.0
gstate = np.kron([tr_g], [cav_g]).astype(complex)
target = np.kron([tr_g], [cav_3]).astype(complex)
# Top (truncated) level of the cavity
cavitytop = np.zeros(dimc)
cavitytop[-1] = 1.0
```

Next, we set up and perform the optimization procedure. In the following cell, time-dependent transmon and cavity drive controls are characterized by I and Q piecewise-constant drive variables. The smooth drives that would be applied to the experiment are produced by applying a sinc filter to the drive variables. The state overlap fidelity is obtained using the dynamics induced by these smooth drive controls.

The infidelity is the first term of the cost function for optimization, which is minimized to provide high-performance controls. In the cell below, there is a second term in the cost from the occupation of the highest (truncated) oscillator number state during the operation to prevent leakage and truncation interference in the optimized dynamics.

```
with qctrl.create_graph() as graph:
# Set up a sinc kernel to bandlimit the pulses (if desired)
sinc_kernel = qctrl.operations.sinc_convolution_kernel(sinc_cutoff_frequency_drives)
# Set up optimizer variables for the i component of the transmon drive
drive_iT_vars = qctrl.operations.optimization_variable(
count=number_of_optimizer_vars,
lower_bound=-max_rabi_rate_T,
upper_bound=max_rabi_rate_T,
)
# Build the pulse signal with segments based on the optimizer variables
drive_iT_raw = qctrl.operations.pwc_signal(
values=drive_iT_vars, duration=gate_duration
)
# Apply the filter to the raw pulse
drive_iT_filtered = qctrl.operations.convolve_pwc(
pwc=drive_iT_raw, kernel=sinc_kernel
)
# Discretize the filtered drive into the desired number of segments
drive_iT_signal = qctrl.operations.discretize_stf(
stf=drive_iT_filtered,
duration=gate_duration,
segment_count=number_of_segments,
)
# Set up the q component of the drive in a similar way
drive_qT_vars = qctrl.operations.optimization_variable(
count=number_of_optimizer_vars,
lower_bound=-max_rabi_rate_T,
upper_bound=max_rabi_rate_T,
)
drive_qT_raw = qctrl.operations.pwc_signal(
values=drive_qT_vars, duration=gate_duration
)
drive_qT_filtered = qctrl.operations.convolve_pwc(
pwc=drive_qT_raw, kernel=sinc_kernel
)
drive_qT_signal = qctrl.operations.discretize_stf(
stf=drive_qT_filtered,
duration=gate_duration,
segment_count=number_of_segments,
)
# Set up the cavity drive components
drive_iC_vars = qctrl.operations.optimization_variable(
count=number_of_optimizer_vars,
lower_bound=-max_rabi_rate_C,
upper_bound=max_rabi_rate_C,
)
drive_iC_raw = qctrl.operations.pwc_signal(
values=drive_iC_vars, duration=gate_duration
)
drive_iC_filtered = qctrl.operations.convolve_pwc(
pwc=drive_iC_raw, kernel=sinc_kernel
)
drive_iC_signal = qctrl.operations.discretize_stf(
stf=drive_iC_filtered,
duration=gate_duration,
segment_count=number_of_segments,
)
drive_qC_vars = qctrl.operations.optimization_variable(
count=number_of_optimizer_vars,
lower_bound=-max_rabi_rate_C,
upper_bound=max_rabi_rate_C,
)
drive_qC_raw = qctrl.operations.pwc_signal(
values=drive_qC_vars, duration=gate_duration
)
drive_qC_filtered = qctrl.operations.convolve_pwc(
pwc=drive_qC_raw, kernel=sinc_kernel
)
drive_qC_signal = qctrl.operations.discretize_stf(
stf=drive_qC_filtered,
duration=gate_duration,
segment_count=number_of_segments,
)
# Combine the I and Q components
driveT_signal = drive_iT_signal + 1j * drive_qT_signal
driveT_signal.name = "$\gamma_T$"
driveC_signal = drive_iC_signal + 1j * drive_qC_signal
driveC_signal.name = "$\gamma_C$"
# Build the time-dependent system Hamiltonian terms
H_drive_T = qctrl.operations.pwc_operator_hermitian_part(
driveT_signal * b.T + driveC_signal * a.T
)
# Construct the total Hamiltonian
hamiltonian = H_sys + H_drive_T
noise_list = []
# Construct the state preparation cost
unitaries = qctrl.operations.time_evolution_operators_pwc(
hamiltonian=hamiltonian,
sample_times=np.linspace(
gate_duration / number_of_segments, gate_duration, number_of_segments
),
)
states = unitaries @ gstate.T
states = states[:, :, 0]
infidelity = qctrl.operations.abs(
1
- qctrl.operations.abs(
qctrl.operations.sum(qctrl.operations.conjugate(target) * states[-1])
)
** 2,
name="infidelity",
)
cavity_top_states = states @ (np.kron(np.eye(dimt), [cavitytop]).astype(complex)).T
cavity_top_pops = qctrl.operations.sum(
qctrl.operations.abs(cavity_top_states) ** 2, name="cavity_top"
)
cost = infidelity + cavity_top_pops
cost.name = "cost"
start_time = time.time()
result = qctrl.functions.calculate_optimization(
graph=graph,
cost_node_name="cost",
output_node_names=["$\gamma_T$", "$\gamma_C$", "infidelity", "cost", "cavity_top"],
optimization_count=10,
target_cost=5e-3,
)
print("Run time (s):", time.time() - start_time)
print("Infidelity:", result.output["infidelity"]["value"])
print(
"Sum over sampled top-level cavity population:",
result.output["cavity_top"]["value"],
)
```

```
controls = {
"$\gamma_T$": result.output["$\gamma_T$"],
"$\gamma_C$": result.output["$\gamma_C$"],
}
plot_controls(plt.figure(), controls, polar=False)
```

The plots display the optimised drive pulses $\gamma_{T,C}(t)$ applied on the transmon qubit and the cavity, respectively.

### Simulating optimized controls to evaluate higher cavity levels

We now simulate the operation using higher dimensions of the transmon qubit and the cavity to evaluate the operation performance. This ensures that the optimizer is not exploiting boundaries of the truncated space, and allows us to evaluate the infidelity including errors from leakage out of the lower-dimensional Hilbert space used for the optimization.

```
dimt = 3 # Higher transmon dimensions
dimc = 18 # Higher cavity dimensions
# System operators defined with higher dimensions
a = np.kron(np.eye(dimt), np.diag(np.sqrt(np.arange(1, dimc)), k=1))
b = np.kron(np.diag(np.sqrt(np.arange(1, dimt)), k=1), np.eye(dimc))
H_sys = (
delta * b.T @ b
+ K / 2 * a.T @ a.T @ a @ a
+ alpha / 2 * b.T @ b.T @ b @ b
+ chi * b.T @ b @ a.T @ a
)
# Initial and target states
tr_g = np.zeros(dimt)
tr_g[0] = 1.0
cav_g = np.zeros(dimc)
cav_g[0] = 1.0
cav_3 = np.zeros(dimc)
cav_3[3] = 1.0
gstate = np.kron([tr_g], [cav_g]).astype(complex)
target = np.kron([tr_g], [cav_3]).astype(complex)
# Use optimized drives
drive_T_values = np.array([seg["value"] for seg in result.output["$\gamma_T$"]])
drive_C_values = np.array([seg["value"] for seg in result.output["$\gamma_C$"]])
# Perform the simulation
with qctrl.create_graph() as graph:
# Build the time-dependent system Hamiltonian terms
driveT_signal = qctrl.operations.pwc_signal(
values=drive_T_values,
duration=gate_duration,
)
driveC_signal = qctrl.operations.pwc_signal(
values=drive_C_values,
duration=gate_duration,
)
H_drive_T = qctrl.operations.pwc_operator_hermitian_part(
driveT_signal * b.T + driveC_signal * a.T
)
# Construct the total Hamiltonian
hamiltonian = H_sys + H_drive_T
# Calculate states over time
unitaries = qctrl.operations.time_evolution_operators_pwc(
hamiltonian=hamiltonian,
sample_times=np.linspace(0, gate_duration, int(gate_duration * 1e9) + 1),
)
states = unitaries @ gstate.T
states.name = "states"
simulation = qctrl.functions.calculate_graph(
graph=graph,
output_node_names=["states"],
)
print(
"Infidelity:",
1 - np.abs(target @ simulation.output["states"]["value"][-1])[0][0] ** 2,
)
```

```
cavity_pops = []
for step in simulation.output["states"]["value"]:
cavity_pops.append([])
for n in range(dimc):
cavityn = np.zeros(dimc)
cavityn[n] = 1.0
cavityn_p = np.sum(np.abs(np.kron(np.eye(dimt), [cavityn]).dot(step)) ** 2)
cavity_pops[-1].append(cavityn_p)
fig, ax = plt.subplots(figsize=(12, 20))
im = ax.imshow(np.array(cavity_pops).T, aspect=13)
plt.ylabel("Number state")
plt.xlabel("Time (ns)")
cbar = plt.colorbar(im, shrink=0.2)
cbar.ax.set_ylabel("Population")
plt.title("Cavity dynamics")
plt.show()
```

In the above plot of the cavity number-state dynamics, observe the population converging to the target state $|0,2\rangle$ by the end of the operation.

```
transmon_pops = []
for step in simulation.output["states"]["value"]:
transmon_pops.append([])
for n in range(dimt):
transmn = np.zeros(dimt)
transmn[n] = 1.0
transmn_p = np.sum(np.abs(np.kron([transmn], np.eye(dimc)).dot(step) ** 2))
transmon_pops[-1].append(transmn_p)
fig, ax = plt.subplots(figsize=(8, 5))
for n in range(dimt):
ax.plot(np.array(transmon_pops).T[n], label=n)
plt.legend(title="Number state")
plt.ylabel("Population")
plt.xlabel("Time (ns)")
plt.title("Transmon dynamics")
plt.show()
```

The above plot displays the transmon number-state population dynamics, and for the given drives only the lowest two transmon levels are substantially occupied during the operation.

## Generate a Fock state using a tunable transmon

In contrast to the above section in the dispersive regime, here we generate a Fock state by tuning the qubit frequency in and out of resonance with the cavity resonator. The dynamics in this case is described by the following Hamiltonian, rotated into the frame of the transmon drive frequency $\nu_T$:

$$H = \delta_T (t) b^\dagger b + \delta_C a^\dagger a + \frac{K}{2} (a^\dagger)^2 a^2 + \frac{\alpha}{2} (b^\dagger)^2 b^2 + \left(\frac{\Omega}{2} a b^\dagger + \gamma_T (t) b + H.c. \right) $$such that $\delta_T (t) = \omega_T (t) - \nu_T$, $\delta_C = \omega_C - \nu_T$, and $\Omega$ is the strength of the qubit-resonator interaction. Here we exploit the tunable transmon qubit frequency and drive to prepare the Fock state; a cavity drive could similarly be added for additional control freedom.

In this section we prepare a target state $|0, 3\rangle$, starting from the initial ground state $|0, 0\rangle$. First, we optimize and validate a state preparation scheme with a fixed drive frequency that is off resonance with the resonator as in (Hofheinz et al., 2008). Then, we add another control to leverage the complex interactions in the system: we optimize the (static) drive frequency in addition to the drive modulation to generate and validate faster Fock state preparation.

### Creating optimized controls with an off-resonant transmon drive

In this section we apply the optimizer to obtain smooth controls that produce the target state. First, we set the parameters for drive frequency, the optimizer, and the target operation. Note that in the cell below, the transmon and cavity dimensions are set higher and lower, respectively, than in the previous optimization. This is to accomodate the variable transmon frequency and resulting higher-level transitions, along with lower cavity occupation numbers without a cavity drive.

```
omega = 2 * np.pi * 36e6
nu_T = 175 * omega # Fixed transmon drive frequency
delta_C = 2 * np.pi * 6.57e9 - nu_T # Cavity detuning
dimt = 4 # Transmon dimensions
dimc = 6 # Cavity dimensions
gate_duration = 128e-9 # s
max_rabi_rate_T = 2 * np.pi * 12e6
number_of_optimizer_vars = 64
number_of_segments = 128
sinc_cutoff_frequency_drives = 2 * np.pi * 20e7
# Specify system operators
a = np.kron(np.eye(dimt), np.diag(np.sqrt(np.arange(1, dimc)), k=1))
b = np.kron(np.diag(np.sqrt(np.arange(1, dimt)), k=1), np.eye(dimc))
H_sys = (
delta_C * a.T @ a
+ K / 2 * a.T @ a.T @ a @ a
+ alpha / 2 * b.T @ b.T @ b @ b
+ omega / 2 * (b.T @ a + b @ a.T)
)
detuning_T = b.T @ b
# Define target operation
tr_g = np.zeros(dimt)
tr_g[0] = 1.0
cav_g = np.zeros(dimc)
cav_g[0] = 1.0
cav_3 = np.zeros(dimc)
cav_3[3] = 1.0
gstate = np.kron([tr_g], [cav_g]).astype(complex)
target = np.kron([tr_g], [cav_3]).astype(complex)
```

As before, we now perform the optimization procedure. In the following cell, the time-dependent transmon qubit frequency is characterized and filtered, in addition to the transmon drive. Unlike the optimization above, the oscillator truncation term is omitted from the cost as it is usually not necessary in this case without a cavity drive.

```
with qctrl.create_graph() as graph:
# Set up a sinc kernel to bandlimit the pulses (if desired)
sinc_kernel = qctrl.operations.sinc_convolution_kernel(sinc_cutoff_frequency_drives)
# Set up qubit detuning variable
delta_T_vars = qctrl.operations.optimization_variable(
count=number_of_optimizer_vars,
lower_bound=-0.05e9 * 2 * np.pi,
upper_bound=0.35e9 * 2 * np.pi,
)
delta_T_raw = qctrl.operations.pwc_signal(
values=delta_T_vars,
duration=gate_duration,
)
delta_T_filtered = qctrl.operations.convolve_pwc(
pwc=delta_T_raw,
kernel=sinc_kernel,
)
delta_T_signal = qctrl.operations.discretize_stf(
stf=delta_T_filtered,
duration=gate_duration,
segment_count=number_of_segments,
name="$\delta_T$",
)
# Set up transmon drive
drive_iT_vars = qctrl.operations.optimization_variable(
count=number_of_optimizer_vars,
lower_bound=-max_rabi_rate_T,
upper_bound=max_rabi_rate_T,
)
drive_iT_raw = qctrl.operations.pwc_signal(
values=drive_iT_vars, duration=gate_duration
)
drive_iT_filtered = qctrl.operations.convolve_pwc(
pwc=drive_iT_raw, kernel=sinc_kernel
)
drive_iT_signal = qctrl.operations.discretize_stf(
stf=drive_iT_filtered,
duration=gate_duration,
segment_count=number_of_segments,
)
drive_qT_vars = qctrl.operations.optimization_variable(
count=number_of_optimizer_vars,
lower_bound=-max_rabi_rate_T,
upper_bound=max_rabi_rate_T,
)
drive_qT_raw = qctrl.operations.pwc_signal(
values=drive_qT_vars, duration=gate_duration
)
drive_qT_filtered = qctrl.operations.convolve_pwc(
pwc=drive_qT_raw, kernel=sinc_kernel
)
drive_qT_signal = qctrl.operations.discretize_stf(
stf=drive_qT_filtered,
duration=gate_duration,
segment_count=number_of_segments,
)
driveT_signal = drive_iT_signal + 1j * drive_qT_signal
driveT_signal.name = "$\gamma_T$"
# Construct the Hamiltonian
H_drive_T = qctrl.operations.pwc_operator_hermitian_part(driveT_signal * b.T)
H_detuning_T = delta_T_signal * detuning_T
hamiltonian = H_sys + H_drive_T + H_detuning_T
noise_list = []
# Set up the state preparation cost
unitaries = qctrl.operations.time_evolution_operators_pwc(
hamiltonian=hamiltonian,
sample_times=np.array([gate_duration]),
)
states = unitaries @ gstate.T
states = states[:, :, 0]
infidelity = qctrl.operations.abs(
1
- qctrl.operations.abs(
qctrl.operations.sum(qctrl.operations.conjugate(target) * states[-1])
)
** 2,
name="infidelity",
)
start_time = time.time()
result = qctrl.functions.calculate_optimization(
graph=graph,
cost_node_name="infidelity",
output_node_names=["$\gamma_T$", "$\delta_T$", "infidelity"],
optimization_count=10,
)
print("Run time (s):", time.time() - start_time)
print("Infidelity:", result.output["infidelity"]["value"])
```

```
controls = {
"$\delta_T$": result.output["$\delta_T$"],
"$\gamma_T$": result.output["$\gamma_T$"],
}
plot_controls(plt.figure(), controls, polar=False)
```

The plots display the transmon detuning from the drive frequency $\delta_T(t)$, as well as the optimised transmon drive $\gamma_{T}(t)$.

```
dimt = 5 # Higher transmon dimensions
dimc = 9 # Higher cavity dimensions
# System operators defined with higher dimensions
a = np.kron(np.eye(dimt), np.diag(np.sqrt(np.arange(1, dimc)), k=1))
b = np.kron(np.diag(np.sqrt(np.arange(1, dimt)), k=1), np.eye(dimc))
H_sys = (
delta_C * a.T @ a
+ K / 2 * a.T @ a.T @ a @ a
+ alpha / 2 * b.T @ b.T @ b @ b
+ omega / 2 * (b.T @ a + b @ a.T)
)
detuning_T = b.T @ b
detuning_C = a.T @ a
# Initial and target states
tr_g = np.zeros(dimt)
tr_g[0] = 1.0
cav_g = np.zeros(dimc)
cav_g[0] = 1.0
cav_3 = np.zeros(dimc)
cav_3[3] = 1.0
gstate = np.kron([tr_g], [cav_g]).astype(complex)
target = np.kron([tr_g], [cav_3]).astype(complex)
# Use optimized drives
drive_T_values = np.array([seg["value"] for seg in result.output["$\gamma_T$"]])
delta_T_values = np.array([seg["value"] for seg in result.output["$\delta_T$"]])
# Perform the simulation
with qctrl.create_graph() as graph:
# Build the time-dependent system Hamiltonian terms
driveT_signal = qctrl.operations.pwc_signal(
values=drive_T_values,
duration=gate_duration,
)
deltaT_signal = qctrl.operations.pwc_signal(
values=delta_T_values,
duration=gate_duration,
)
H_drive_T = qctrl.operations.pwc_operator_hermitian_part(driveT_signal * b.T)
H_detuning_T = deltaT_signal * detuning_T
# Construct the total Hamiltonian
hamiltonian = H_sys + H_drive_T + H_detuning_T
# Calculate states over time
unitaries = qctrl.operations.time_evolution_operators_pwc(
hamiltonian=hamiltonian,
sample_times=np.linspace(0, gate_duration, int(gate_duration * 1e9) + 1),
)
states = unitaries @ gstate.T
states.name = "states"
simulation = qctrl.functions.calculate_graph(
graph=graph,
output_node_names=["states"],
)
print(
"Infidelity:",
1 - np.abs(target @ simulation.output["states"]["value"][-1])[0][0] ** 2,
)
```

```
cavity_pops = []
for step in simulation.output["states"]["value"]:
cavity_pops.append([])
for n in range(dimc):
cavityn = np.zeros(dimc)
cavityn[n] = 1.0
cavityn_p = np.sum(np.abs(np.kron(np.eye(dimt), [cavityn]).dot(step)) ** 2)
cavity_pops[-1].append(cavityn_p)
fig, ax = plt.subplots(figsize=(12, 20))
im = ax.imshow(np.array(cavity_pops).T, aspect=5)
plt.ylabel("Number state")
plt.xlabel("Time (ns)")
cbar = plt.colorbar(im, shrink=0.25)
cbar.ax.set_ylabel("Population")
plt.title("Cavity dynamics")
plt.show()
```

In the above plot of the cavity number-state dynamics, observe the population converging to the target state $|0,3\rangle$ by the end of the operation. The population seems to increase in steps using the optimized controls.

```
transmon_pops = []
for step in simulation.output["states"]["value"]:
transmon_pops.append([])
for n in range(dimt):
transmn = np.zeros(dimt)
transmn[n] = 1.0
transmn_p = np.sum(np.abs(np.kron([transmn], np.eye(dimc)).dot(step) ** 2))
transmon_pops[-1].append(transmn_p)
fig, ax = plt.subplots(figsize=(8, 5))
for n in range(dimt):
ax.plot(np.array(transmon_pops).T[n], label=n)
plt.legend(title="Number state")
plt.ylabel("Population")
plt.xlabel("Time (ns)")
plt.title("Transmon dynamics")
plt.show()
```

In the above plot of the transmon number-state dynamics, the population oscillates between the $|0\rangle$ and $|1\rangle$ number states. This is correlated with the cavity dynamics above: the transmon excitation is induced by the drive away from the resonant interaction regime with the cavity, and the cavity population can subsequently increase when the transmon is tuned into the resonant regime. These dynamics are expected with the available controls, and are similar to those displayed by the analytical Fock state generation scheme in (Hofheinz et al., 2008). In the next section, we add another control degree of freedom to allow more flexible and faster operation dynamics.

### Creating optimized controls with an optimized-frequency transmon drive

Here we add an extra degree of control over the Fock state generation scheme by optimizing the frequency of the transmon drive. Again, the drive frequency is fixed, and this could optionally be made dynamic for additional control. As before, we specify the system and optimization parameters before optimizing the controls.

```
dimt = 4 # Transmon dimensions
dimc = 6 # Cavity dimensions
gate_duration = 100e-9 # s
max_rabi_rate_T = 2 * np.pi * 12e6
number_of_optimizer_vars = 64
number_of_segments = 128
sinc_cutoff_frequency_drives = 2 * np.pi * 20e7
# Set up system operators
a = np.kron(np.eye(dimt), np.diag(np.sqrt(np.arange(1, dimc)), k=1))
b = np.kron(np.diag(np.sqrt(np.arange(1, dimt)), k=1), np.eye(dimc))
H_sys = (
K / 2 * a.T @ a.T @ a @ a
+ alpha / 2 * b.T @ b.T @ b @ b
+ omega / 2 * (b.T @ a + b @ a.T)
)
detuning_T = b.T @ b
detuning_C = a.T @ a
# Define target operation
tr_g = np.zeros(dimt)
tr_g[0] = 1.0
cav_g = np.zeros(dimc)
cav_g[0] = 1.0
cav_3 = np.zeros(dimc)
cav_3[3] = 1.0
gstate = np.kron([tr_g], [cav_g]).astype(complex)
target = np.kron([tr_g], [cav_3]).astype(complex)
```

```
with qctrl.create_graph() as graph:
sinc_kernel = qctrl.operations.sinc_convolution_kernel(sinc_cutoff_frequency_drives)
# Set up drive frequency optimization
nu_T_var = qctrl.operations.optimization_variable(
count=1,
lower_bound=2 * np.pi * 6.25e9,
upper_bound=2 * np.pi * 6.57e9,
)
nu_T_signal = qctrl.operations.pwc_signal(
values=nu_T_var,
duration=gate_duration,
)
# Define cavity detuning from the drive frequency
delta_C = qctrl.operations.subtract(2 * np.pi * 6.57e9, nu_T_signal)
delta_C.name = "$\delta_C$"
# Set up transmon frequency optimization
delta_T_vars = qctrl.operations.optimization_variable(
count=number_of_optimizer_vars,
lower_bound=-0.05e9 * 2 * np.pi,
upper_bound=0.35e9 * 2 * np.pi,
)
delta_T_raw = qctrl.operations.pwc_signal(
values=delta_T_vars,
duration=gate_duration,
)
delta_T_filtered = qctrl.operations.convolve_pwc(
pwc=delta_T_raw,
kernel=sinc_kernel,
)
delta_T_signal = qctrl.operations.discretize_stf(
stf=delta_T_filtered,
duration=gate_duration,
segment_count=number_of_segments,
name="$\delta_T$",
)
# Set up transmon drive optimization
drive_iT_vars = qctrl.operations.optimization_variable(
count=number_of_optimizer_vars,
lower_bound=-max_rabi_rate_T,
upper_bound=max_rabi_rate_T,
)
drive_iT_raw = qctrl.operations.pwc_signal(
values=drive_iT_vars, duration=gate_duration
)
drive_iT_filtered = qctrl.operations.convolve_pwc(
pwc=drive_iT_raw, kernel=sinc_kernel
)
drive_iT_signal = qctrl.operations.discretize_stf(
stf=drive_iT_filtered,
duration=gate_duration,
segment_count=number_of_segments,
)
drive_qT_vars = qctrl.operations.optimization_variable(
count=number_of_optimizer_vars,
lower_bound=-max_rabi_rate_T,
upper_bound=max_rabi_rate_T,
)
drive_qT_raw = qctrl.operations.pwc_signal(
values=drive_qT_vars, duration=gate_duration
)
drive_qT_filtered = qctrl.operations.convolve_pwc(
pwc=drive_qT_raw, kernel=sinc_kernel
)
drive_qT_signal = qctrl.operations.discretize_stf(
stf=drive_qT_filtered,
duration=gate_duration,
segment_count=number_of_segments,
)
driveT_signal = drive_iT_signal + 1j * drive_qT_signal
driveT_signal.name = "$\gamma_T$"
# Construct the Hamiltonian
H_drive_T = qctrl.operations.pwc_operator_hermitian_part(driveT_signal * b.T)
H_detuning_T = delta_T_signal * detuning_T
H_C = delta_C * detuning_C
hamiltonian = H_sys + H_drive_T + H_detuning_T + H_C
noise_list = []
# Construct the state preparation cost
unitaries = qctrl.operations.time_evolution_operators_pwc(
hamiltonian=hamiltonian,
sample_times=np.array([gate_duration]),
)
states = unitaries @ gstate.T
states = states[:, :, 0]
infidelity = qctrl.operations.abs(
1
- qctrl.operations.abs(
qctrl.operations.sum(qctrl.operations.conjugate(target) * states[-1])
)
** 2,
name="infidelity",
)
start_time = time.time()
result = qctrl.functions.calculate_optimization(
graph=graph,
cost_node_name="infidelity",
output_node_names=["$\gamma_T$", "$\delta_T$", "$\delta_C$", "infidelity"],
optimization_count=10,
)
print("Run time (s):", time.time() - start_time)
print("Infidelity:", result.output["infidelity"]["value"])
```

```
print(
"Cavity detuning from drive (GHz):",
result.output["$\delta_C$"][0]["value"] / 2 / np.pi / 1e9,
)
controls = {
"$\delta_T$": result.output["$\delta_T$"],
"$\gamma_T$": result.output["$\gamma_T$"],
}
plot_controls(plt.figure(), controls, polar=False)
```

The plots display the transmon detuning from the drive frequency $\delta_T(t)$, as well as the optimised transmon drive $\gamma_{T}(t)$.

```
dimt = 5 # Higher transmon dimensions
dimc = 9 # Higher cavity dimensions
# System operators defined with higher dimensions
a = np.kron(np.eye(dimt), np.diag(np.sqrt(np.arange(1, dimc)), k=1))
b = np.kron(np.diag(np.sqrt(np.arange(1, dimt)), k=1), np.eye(dimc))
H_sys = (
K / 2 * a.T @ a.T @ a @ a
+ alpha / 2 * b.T @ b.T @ b @ b
+ omega / 2 * (b.T @ a + b @ a.T)
)
detuning_T = b.T @ b
detuning_C = a.T @ a
# Initial and target states
tr_g = np.zeros(dimt)
tr_g[0] = 1.0
cav_g = np.zeros(dimc)
cav_g[0] = 1.0
cav_3 = np.zeros(dimc)
cav_3[3] = 1.0
gstate = np.kron([tr_g], [cav_g]).astype(complex)
target = np.kron([tr_g], [cav_3]).astype(complex)
# Use optimized drives
drive_T_values = np.array([seg["value"] for seg in result.output["$\gamma_T$"]])
delta_T_values = np.array([seg["value"] for seg in result.output["$\delta_T$"]])
delta_C_values = np.array([seg["value"] for seg in result.output["$\delta_C$"]])
# Perform the simulation
with qctrl.create_graph() as graph:
# Build the time-dependent system Hamiltonian terms
driveT_signal = qctrl.operations.pwc_signal(
values=drive_T_values,
duration=gate_duration,
)
deltaT_signal = qctrl.operations.pwc_signal(
values=delta_T_values,
duration=gate_duration,
)
deltaC_signal = qctrl.operations.pwc_signal(
values=delta_C_values,
duration=gate_duration,
)
H_drive_T = qctrl.operations.pwc_operator_hermitian_part(driveT_signal * b.T)
H_detuning_T = deltaT_signal * detuning_T
H_detuning_C = deltaC_signal * detuning_C
# Construct the total Hamiltonian
hamiltonian = H_sys + H_drive_T + H_detuning_T + H_detuning_C
# Calculate states over time
unitaries = qctrl.operations.time_evolution_operators_pwc(
hamiltonian=hamiltonian,
sample_times=np.linspace(0, gate_duration, int(gate_duration * 1e9) + 1),
)
states = unitaries @ gstate.T
states.name = "states"
simulation = qctrl.functions.calculate_graph(
graph=graph,
output_node_names=["states"],
)
print(
"Infidelity:",
1 - np.abs(target @ simulation.output["states"]["value"][-1])[0][0] ** 2,
)
```

```
cavity_pops = []
for step in simulation.output["states"]["value"]:
cavity_pops.append([])
for n in range(dimc):
cavityn = np.zeros(dimc)
cavityn[n] = 1.0
cavityn_p = np.sum(np.abs(np.kron(np.eye(dimt), [cavityn]).dot(step)) ** 2)
cavity_pops[-1].append(cavityn_p)
fig, ax = plt.subplots(figsize=(12, 14))
im = ax.imshow(np.array(cavity_pops).T, aspect=4.5)
plt.ylabel("Number state")
plt.xlabel("Time (ns)")
cbar = plt.colorbar(im, shrink=0.35)
cbar.ax.set_ylabel("Population")
plt.title("Cavity dynamics")
plt.show()
```

The cavity number-state populations displayed in the figure exhibit faster accumulation with the optimized drive frequency. The dynamics are more complex than the step-wise increase observed for the off-resonant drive above.

```
transmon_pops = []
for step in simulation.output["states"]["value"]:
transmon_pops.append([])
for n in range(dimt):
transmn = np.zeros(dimt)
transmn[n] = 1.0
transmn_p = np.sum(np.abs(np.kron([transmn], np.eye(dimc)).dot(step) ** 2))
transmon_pops[-1].append(transmn_p)
fig, ax = plt.subplots(figsize=(8, 5))
for n in range(dimt):
ax.plot(np.array(transmon_pops).T[n], label=n)
plt.legend(title="Number state")
plt.ylabel("Population")
plt.xlabel("Time (ns)")
plt.title("Transmon dynamics")
plt.show()
```

The transmon dynamics are similarly more complex than for the off-resonant drive, and involve higher-level excitations that contribute to faster Fock state generation in the cavity. Here the optimized drive frequency permits the state generation in 100ns, in contrast to the 128ns scheme without this optimization. BOULDER OPAL optimization tools thus exploit the full system dynamics to generate optimal controls; further control freedoms could similarly be added to enhance the scheme speed or robustness to noise sources.