Simulation & optimal control of a two-qubit gate on nuclear spins in an NV quantum processor¶
In this tutorial, we’ll show you how to simulate and optimise a two-qubit controlled-Z (CZ) gate on nuclear spins in an NV quantum processor using the qruise-toolset. The system consists of an NV centre coupled to two nuclear spins: ${}^{13}\text{C}$ and ${}^{15}\text{N}$.
Note: We recommend the reader to first refer to the tutorial on Simulation & optimal control of a single-qubit gate in an NV quantum processor before continuing with this tutorial.
The tutorial is comprised of the following:
- Defining the Hamiltonian (theory)
- Defining the stationary Hamiltonian
- Defining the external drive
- Solving the Schrödinger equation
- Performing quantum optimal control
1. Defining the Hamiltonian (theory)¶
The general Hamiltonian for an NV centre coupled to $i$ nuclear spins (where $i$ is the index of each spin) is:
$$ H = \underbrace{D \left(S_z^2 - \tfrac{3}{2}\right)\vphantom{\sum}}_{\text{zero-field splitting}} + \underbrace{\gamma_e B_0 S_z \vphantom{\sum}}_{\text{NV Zeeman interaction}} + \underbrace{\gamma_e B_1^e(t) S_x \vphantom{\sum}}_{\text{NV drive interaction}} + \sum\limits_{i} \Big[\underbrace{\gamma_i B_0 J_{iz} \vphantom{\sum}}_{\text{nuclear Zeeman interaction}} + \underbrace{\gamma_i B_{1}^n(t) J_{ix} \vphantom{\sum}}_{\text{nuclear drive interaction}} + \underbrace{ S \cdot A_i \cdot J_i \vphantom{\sum}}_{\text{hyperfine interaction}} + \underbrace{Q_i \left(J_{iz}^2 - \tfrac{3}{2}\right) \vphantom{\sum}}_{\text{quadrupole interaction}} \Big]. \tag{1} $$
Here, $S_{x,y,z}$ are the spin-$1$ operators for the NV electronic spin, and $J_{i,x,y,z}$ are the spin-$\frac{1}{2}$ operators for the $i^\text{th}$ nuclear spin.
If you're not familiar with this Hamiltonian, please refer to the theory section at the start of our tutorial on Simulation & optimal control of a single-qubit gate in an NV quantum processor.
Single-qubit gates are performed within the computational subspace, which corresponds to the $|m_s = -1\rangle$ state (see diagram below), using spectrally selective radio-frequency pulses. To implement a CZ two-qubit gate, however, we also need to make use of the auxiliary subspace spanned by the $|m_s = 0\rangle$ state.
The CZ gate is implemented by applying a selective microwave $2\pi$ pulse to the NV electron spin, but only when both nuclear spins are in the $|1\rangle$ state. This pulse causes the system to briefly enter the auxiliary subspace ($|m_s = 0\rangle$) and then return to the computational subspace ($|m_s = -1\rangle$). The effect is to add a $\pi$ phase to the $|11\rangle$ state, leaving all other states unchanged. This conditional phase is what enables entanglement between the two qubits. For more details on this gate, see Waldherr et al..
Now let's consider the quantum system we're working with: one NV centre spin coupled to two nuclear spins, a ${}^{13}\text{C}$ and a ${}^{15}\text{N}$. We can simplify our Hamiltonian for a CZ gate on nuclei with the following simplifications & approximations:
Two-level approximation: The CZ gate only involves the $|m_s = 0\rangle$ and $|m_s = -1\rangle$ subspaces of the NV centre spin, so the $|m_s = +1\rangle$ state is excluded. As a result, the NV spin can be treated as an effective two-level system where the spin-1 operators are replaced by Pauli operators ($\sigma_{z,x}$).
No quadrupole interaction: Both ${}^{13}\text{C}$ and ${}^{15}\text{N}$ have spin-$\tfrac{1}{2}$, so the quadrupole term vanishes.
Secular approximation: Under a strong static magnetic field, $B_0$, aligned with the NV axis ($z$-axis), non-energy-conserving terms such as $S_x$, $S_y$, $J_{ix}$, and $J_{iy}$ can be neglected in the static part of the Hamiltonian. This reduces the hyperfine interaction to $A_{i,zz} S_z J_{iz}$.
Aligned nuclei: We consider only nuclei that are nearly aligned with the NV axis. This ensures that off-diagonal hyperfine components like $A_{i,xz}$ and $A_{i,yz}$ are small, so don’t appear explicitly in the Hamiltonian (though they still affect the energy level shifts).
Applying all these simplifications, the full Hamiltonian reduces to the effective lab-frame form used for nuclear CZ gate implementations:
$$ H = \frac{\Delta}{2} \sigma_z + \frac{\Omega}{2} B_1^e(t) \sigma_x + \sum_i (\alpha_i + \beta_i \sigma_z) J_{iz} + \gamma_i B_1^n(t) J_{ix} \tag{2} $$
with:
$\sigma_{z,x}$: Pauli operators acting on the NV spin
$\Delta= D - \gamma_e B_0 $: effective NV detuning
$\Omega$: effective NV drive coupling strength
$\alpha_i= -\frac{\gamma_i B_0}{2} - \frac{\omega_i}{2}$: effective nuclear Zeeman shift for nuclear spin $i$
$\beta_i= -\frac{\gamma_i B_0}{2} + \frac{\omega_i}{2}$: effective conditional coupling between the NV spin and nuclear spin $i$
$\omega_i= \sqrt{(A_{i,zz} + \gamma_i B_0)^2 + A_{i,xz}^2 + A_{i,yz}^2}$: effective nuclear transition frequency due to combined hyperfine and Zeeman interactions.
For more details on this derivation, please see Chen et al..
Tip: Make sure that your environment is correctly set to handle float64 precision by setting JAX_ENABLE_X64=True or add the following codeblock to your scripts' preamble:
# Enable float64 precision
import jax
jax.config.update("jax_enable_x64", True)
2. Defining the stationary Hamiltonian¶
First, we’ll set up some general parameters for the simulation of our qubit. These include the strength of the static magnetic field ($B_0$), the nuclear hyperfine coupling constants ($A_i$), and the gyromagnetic ratios ($\gamma_i$).
# NV constants
D = 2.87e9 # zero field splitting (Hz)
gamma_e = 28.024e9 # NV gyromagnetic ratio (Hz/T)
B0 = 0.1 # static magnetic field (T)
Delta = D - gamma_e * B0 # effective NV detuning (Hz)
# hyperfine coupling constants for 13C and 15N (Hz)
A_zz_C13 = 2.281e6
A_xz_C13 = 0.24e6
A_yz_C13 = 0.24e6
A_zz_N15 = 3.03e6
A_xz_N15 = 0.0 # zero due to alignment with NV axis
A_yz_N15 = 0.0 # zero due to alignment with NV axis
# gyromagnetic ratios for 13C and 15N (Hz/T)
gamma_C13 = 10.7084e6
gamma_N15 = -4.316e6
gamma_ratio = gamma_N15 / gamma_C13
We can then use these parameters to define $\omega_i$, $\alpha_i$, and $\beta_i$.
import jax.numpy as jnp
# effective nuclear transition frequencies (Hz)
omega_C13 = jnp.sqrt((A_zz_C13 + gamma_C13 * B0) ** 2 + A_xz_C13**2 + A_yz_C13**2)
omega_N15 = jnp.sqrt((A_zz_N15 + gamma_N15 * B0) ** 2 + A_xz_N15**2 + A_yz_N15**2)
# effective nuclear Zeeman shift (Hz)
alpha_C13 = -0.5 * (gamma_C13 * B0 + omega_C13)
alpha_N15 = -0.5 * (gamma_N15 * B0 + omega_N15)
# effective nuclear-NV coupling (Hz)
beta_C13 = -0.5 * (gamma_C13 * B0 - omega_C13)
beta_N15 = -0.5 * (gamma_N15 * B0 - omega_N15)
Before we can define our drive or our Hamiltonian, we need to define our operators. For this, we import the Pauli $\sigma$-operators and the identity operator from the QuTip library. In a multi-spin system, it's essential to ensure each operator acts only on the intended spin while leaving the other spins unaffected. To do this, we use tensor products, where the relevant $\sigma$-operator acts on the target spin and the identity acts on the other two spins.
from qutip import tensor, sigmax, sigmaz, qeye
import qutip as qt
# define identity operator using Hilbert space dimension
hilbert_dim = 2
Id = qeye(hilbert_dim)
# define Pauli operators for each spin
sigmax_NV = tensor(sigmax(), Id, Id)
sigmaz_NV = tensor(sigmaz(), Id, Id)
sigmax_C13 = tensor(Id, sigmax(), Id)
sigmaz_C13 = tensor(Id, sigmaz(), Id)
sigmax_N15 = tensor(Id, Id, sigmax())
sigmaz_N15 = tensor(Id, Id, sigmaz())
Now we can define the stationary part of our Hamiltonian according to Equation 2.
from qruise.toolset import Hamiltonian
H_zeeman_NV = jnp.pi * Delta * sigmaz_NV
H_zeeman_nuclei = jnp.pi * (alpha_C13 * sigmaz_C13 + alpha_N15 * sigmaz_N15)
H_coupling = (
jnp.pi * (beta_C13 * sigmaz_NV * sigmaz_C13 + beta_N15 * sigmaz_NV * sigmaz_N15) / 2
)
H = Hamiltonian(H_zeeman_NV + H_zeeman_nuclei + H_coupling)
3. Defining the external drive¶
For our external drive, we'll use a piecewise constant rectangular pulse with a finite number of time steps of duration $dt$, where each step has a constant amplitude set by the Rabi frequency, $\Omega$.
from qruise.toolset.control_stack import PWCPulse
# time parameters for simulation (s)
t0 = 0.0 # initial time
tfinal = 3e-6 # gate duration
grid_size = 1000 # number of time points
ts = jnp.linspace(t0, tfinal, grid_size)
dt = ts[1] - ts[0] # time step
# Rabi frequency for drive pulse (MHz)
rabi_freq = 1.0
pwc_pulse = PWCPulse(
"drive",
{
"env": (rabi_freq * jnp.ones_like(ts), True),
"t0": (t0, False),
"dt": (float(dt), False),
},
)
params = pwc_pulse.params
Note: The True and False values indicate whether the parameter is optimisable or fixed, respectively.
Recall that the CZ gate on nuclear spins is implemented by applying a frequency-selective $2\pi$ pulse to the NV spin, conditioned on the nuclear spins being in the $|11\rangle$ state. This pulse targets the $|111\rangle \leftrightarrow |011\rangle$ transition, returning the system to $|111\rangle$ with a phase difference of $\pi$ relative to all other computational states.
To get the frequency of this transition, we need to calculate the eigenvalues of the static Hamiltonian. For high-fidelity control, we'll also check that other transitions are well separated. If they're too close, we may need to drive the nuclear spins as well, so we'll also calculate their average transition frequencies.
# calculate eigenvalues of the static Hamiltonian
eigvals = jnp.linalg.eigh(H.H0)[0]
# NV spin transitions
transition_freqs_nv = [
eigvals[0] - eigvals[4], # 000 -> 100
eigvals[1] - eigvals[5], # 001 -> 101
eigvals[2] - eigvals[6], # 010 -> 110
eigvals[3] - eigvals[7], # 011 -> 111
]
# 13C nuclear spin transitions
transition_freqs_c13 = [
eigvals[0] - eigvals[2], # 000 -> 010
eigvals[1] - eigvals[3], # 001 -> 011
eigvals[4] - eigvals[6], # 100 -> 110
eigvals[5] - eigvals[7], # 101 -> 111
]
# 15N nuclear spin transitions
transition_freqs_n15 = [
eigvals[0] - eigvals[1], # 000 -> 001
eigvals[2] - eigvals[3], # 010 -> 011
eigvals[4] - eigvals[5], # 100 -> 101
eigvals[6] - eigvals[7], # 110 -> 111
]
lo_freq_nv = transition_freqs_nv[3] # 011 -> 111 transition frequency
lo_freq_c13 = jnp.mean(
jnp.array(transition_freqs_c13)
) # average of 13C transition frequencies
lo_freq_n15 = jnp.mean(
jnp.array(transition_freqs_n15)
) # average of 15N transition frequencies
We can now define separate drive functions to act on the NV centre and the ${}^{13}\text{C}$ and ${}^{15}\text{N}$ nuclear spins, allowing for easy switching between which one is driven. As an initial guess, each drive uses the uniform-amplitude piecewise constant pulse modulated by a sine wave at the corresponding target frequency.
# NV drive
def pwc_drive_nv(t, params):
"""Envelope: pwc pulse modulated by a sine"""
return 1e6 * pwc_pulse(t, params) * jnp.sin(lo_freq_nv * t)
# 13C drive
def pwc_drive_c13(t, params):
"""Envelope: pwc pulse modulated by a sine"""
return 1e6 * pwc_pulse(t, params) * jnp.sin(lo_freq_c13 * t)
# 15N drive
def pwc_drive_n15(t, params):
"""Envelope: pwc pulse modulated by a sine"""
return 1e6 * pwc_pulse(t, params) * jnp.sin(lo_freq_n15 * t)
It’s a good idea to check your drives behave as desired before proceeding, so let’s quickly plot them.
from qruise.toolset.plots import PlotUtil
ts1 = jnp.linspace(t0, tfinal, grid_size * 10) # finer grid for plotting
canvas = PlotUtil(x_axis_label="t [μs]", y_axis_label="Amplitude [Hz]", notebook=True)
canvas.plot(ts1 / 1e-6, pwc_drive_nv(ts1, params), labels=["NV drive"])
canvas.plot(ts1 / 1e-6, pwc_drive_c13(ts1, params), labels=["C13 drive"])
canvas.plot(ts1 / 1e-6, pwc_drive_n15(ts1, params), labels=["N15 drive"])
canvas.show_canvas()
Great! As expected, each drive exhibits a rectangular envelope modulated by sine functions at the defined drive frequencies.
Now that we our drives, we can add the time-dependent terms to our Hamiltonian.
H.add_term(sigmax_NV * jnp.pi, pwc_drive_nv)
H.add_term(sigmax_C13 * jnp.pi, pwc_drive_c13)
H.add_term(sigmax_N15 * jnp.pi, pwc_drive_n15)
4. Solving the Schrödinger equation¶
To solve the Schrödinger equation, we need to specify the equations and parameters that govern the system. In the qruise-toolset, we combine these using a Problem. To instantiate the Problem, we need:
- the Hamiltonian (
H) - the initial qubit state (
y0) - the pulse parameters (
params) - the time interval of the simulation (
t0totfinal).
Before we can define our initial qubit state, we need to define define and label our basis states. As we have three two-level systems, we have $2^3=8$ basis states: $|000\rangle, |001\rangle, |010\rangle, |011\rangle, |100\rangle, |101\rangle, |110\rangle, |111\rangle$. We'll set our initial qubit state as $|111\rangle$ state.
from qruise.toolset.utils import computational_basis
from qruise.toolset import Problem
# define basis states and label
num_subsystems = 3 # number of subsystems (spins)
basis, labels = computational_basis([hilbert_dim] * num_subsystems)
# define initial qubit state (|111> state)
y0 = basis[labels.index((1, 1, 1))]
# define problem
prob = Problem(H, y0, params, (t0, tfinal))
The user can then define the type of solver they want to use and the equation they want to solve (for example, the Schrödinger equation, the master equation, or the Lindblad equation).
In this case, we’ll choose to use the piece-wise constant solver (PWCSolver) to solve the Schrödinger equation (sepwc). We can then calculate the wavefunction of the system at each timestamp.
from qruise.toolset import PWCSolver
# uses PWCSolver to solve Schrödinger equation
solver = PWCSolver(n=grid_size, store=True)
solver.set_system(prob.sepwc())
times, res_init_state = solver.evolve(*prob.problem())
Note: Setting store=True in the PWCSolver ensures that the results are stored at each timestamp during the simulation.
Finally, we can view the results of our simulation by calculating and plotting the expectation value of the spin state populations against time.
from qruise.toolset import get_population
from qruise.toolset.plots import qruise_colour_palette
# Calculate populations
pop = get_population(res_init_state)
# plot state populations
canvas.init_canvas(
x_axis_label="t [s]",
y_axis_label="Population",
colour_palette=qruise_colour_palette,
)
canvas.plot(times, pop, labels=[str(label) for label in labels])
canvas.show_canvas()
Nice! You've just simulated a nuclear CZ gate in an NV quantum processor.
We observe a decrease in the population of the $|111\rangle$ basis state and an increase in the population of the other states. For a perfect CZ gate, however, we would expect the population of the $|111\rangle$ state to return to 1 and the population of the other states to return to zero at the end of the pulse. Clearly, the gate is imperfect.
To better visualise the gate's imperfections, we can plot the unitary matrix at the end of the simulation time and compare it to the ideal unitary. This allows us also to view the off-diagonal components, which are not visible in the population plot above.
import numpy as np
# calculate unitary from simulation
U0 = jnp.eye(8, dtype=jnp.complex128)
times, res_init_unitary = solver.evolve(U0, params, (t0, tfinal))
U_unopt = res_init_unitary[-1, :] # unitary at t=tfinal
# compute ideal CZ unitary
CZ_ideal = jnp.array(
[
[1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, -1],
],
dtype=np.complex128,
)
# convert unitaries into QuTip quantum objects
dims = [[hilbert_dim] * num_subsystems, [hilbert_dim] * num_subsystems]
Uideal_qopj = qt.Qobj(CZ_ideal, dims=dims) # ideal unitary
Usim_qobj = qt.Qobj(U_unopt, dims=dims) # simulated unitary
Now we can plot the real and imaginary parts of both the ideal and simulated unitary matrices using 3D bar charts. This visual comparison helps us identify discrepancies between the actual gate implemented by our control pulses and the perfect CZ gate, highlighting any unwanted off-diagonal elements or phase errors. We install and use the matplotlib library for this purpose.
%pip install -q matplotlib
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/pty.py:89: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock. pid, fd = os.forkpty()
Note: you may need to restart the kernel to use updated packages.
import matplotlib.pyplot as plt
# create 2×2 subplot layout
fig, axes = plt.subplots(2, 2, figsize=(10, 10), subplot_kw={"projection": "3d"})
# define titles and data
titles = [
"Ideal unitary (real)",
"Ideal unitary (imaginary)",
"Simulated unitary (real)",
"Simulated unitary (imaginary)",
]
qobjs = [Uideal_qopj, Uideal_qopj, Usim_qobj, Usim_qobj]
bar_styles = ["real", "img", "real", "img"]
# prepare each subplot
for ax, title, qobj, style in zip(axes.flat, titles, qobjs, bar_styles):
qt.matrix_histogram(
qobj,
x_basis=range(8),
y_basis=range(8),
bar_style=style,
ax=ax,
limits=[-1, 1],
colorbar=False,
)
ax.set_title(title, fontsize=10)
ax.collections[0].set_cmap("viridis") # Correct way to set colormap
ax.set_zticks([-1, 0, 1])
ax.set_xticklabels([f"{i:03b}" for i in range(8)], fontsize=4)
ax.set_yticklabels([f"{i:03b}" for i in range(8)], fontsize=4)
ax.tick_params(axis="both", labelsize=8)
# adjust layout and show plot
plt.subplots_adjust(wspace=0.3, hspace=0.3)
plt.show()
You can see that the match between the ideal unitary and the simulated unitary is not great. We can perform optimal control to try and improve this.
5. Performing quantum optimal control¶
To perform optimal control, we need a metric that quantifies how far we are from the target CZ gate. Here we define unitary overlap function as $$ \mathcal{L}=1-\mathcal{F}=1-\frac{|\text{Tr}[U(t=t_\text{final})^{\dagger} U_T]|^2 }{d^2}, \tag{7} $$ where $\mathcal{F}$ is the gate fidelity, $U(t=t_\text{final})$ is the achieved unitary at the end time of the pulse, $U_T$ is the target (desired) unitary, and $d = 2^n$, where $n$ is the number of qubits in the system. The objective is to optimise the parameters of the drive pulse by minimising $\mathcal{L}$.
def loss(x, y):
d = 2**num_subsystems
fid = (jnp.abs(jnp.trace(x.T.conj() @ y)) / d) ** 2
return 1 - fid
We are now ready to define a quantum optimal control problem (QOCProblem). Similar to the definition of a Problem for simulating the dynamics of a quantum system, QOCProblem further defines the optimisation protocol. As inputs, it takes
- the Hamiltonian (
H) - the initial unitary (
U0, identity) - the pulse parameters (
params) - the time interval of interest (
t0totfinal) - the desired unitary(
Ut, i.e.,CZ_ideal) - the loss function (
loss)
from qruise.toolset import QOCProblem
Ut = CZ_ideal
opt_prob = QOCProblem(H, U0, params, (t0, tfinal), Ut, loss)
The whole workflow now reduces to simulating the dynamics using the initial guess values for the parameters we defined earlier. At the end of the simulation, we get the new unitary and calculate the infidelity with respect to the target unitary. We then calculate the gradients of the loss function for the initial value of the parameters. Based on the gradient values, we update the parameter values iteratively until they converge and minimise the loss function. This procedure is called GRadient Ascent Pulse Engineering (GRAPE) algorithm.
The Optimiser class is provided for this specific task. As a user, you choose the solver you want to use to simulate the dynamics and which algorithm to use to carry out the optimisation. Here, we use a piecewise constant solver (PWCSolver) and the Broyden–Fletcher–Goldfarb–Shannon (BFGS) method from our OptimistixMinimiser library for optimisation.
from qruise.toolset import Optimiser, OptimistixMinimiser
import optimistix as optx
minimiser = OptimistixMinimiser(optx.BFGS, rtol=1e-4, atol=1e-4)
opt = Optimiser(minimiser, PWCSolver(n=grid_size, eq=prob.sepwc()))
opt.set_optimisation(opt_prob.loss)
result, opt_summary = opt.optimise(*opt_prob.problem())
We can now plot the optimised PWC pulses to see how they compare to the originals.
from qruise.toolset.plots import PlotUtil
canvas = PlotUtil(x_axis_label="t [s]", y_axis_label="Amplitude [Hz]", notebook=True)
canvas.plot(ts, pwc_drive_nv(ts, result), labels=["NV drive"])
canvas.plot(ts, pwc_drive_c13(ts, result), labels=["C13 drive"])
canvas.plot(ts, pwc_drive_n15(ts, result), labels=["N15 drive"])
canvas.show_canvas()
You can see that the optimised drive pulses are much more complex and irregular compared to the originals.
We can now simulate and plot the state populations using the optimised pulses to see if our optimisation worked. We'll also plot the unoptimised populations again so we can compare the results side-by-side.
# simulate with optimised parameters
times, res_opt_state = solver.evolve(y0, result, (t0, tfinal))
# get the state populations after optimisation
opt_pop = get_population(res_opt_state)
# plot unoptimised populations
canvas = PlotUtil(x_axis_label="Time (s)", y_axis_label="Population", notebook=True)
canvas.plot(ts, pop, labels=[str(label) for label in labels])
canvas.canvas.title.text = "Unoptimised populations"
canvas.show_canvas()
# plot optimised populations
canvas = PlotUtil(x_axis_label="Time (s)", y_axis_label="Population", notebook=True)
canvas.plot(times, opt_pop, labels=[str(label) for label in labels])
canvas.canvas.title.text = "Optimised populations"
canvas.show_canvas()
We can see here that the population of $|111\rangle$ returns back to 1 after a $2\pi$ rotation and the final populations of the other states has returned to zero. This is what we'd expect from a CZ gate; however, these plots don't tell us anything about the phase. We can again visualise this by plotting the unitaries.
# calculate unitary from simulation
times, res_opt_unitary = solver.evolve(U0, result, (t0, tfinal))
U_opt = res_opt_unitary[-1, :] # unitary at t=tfinal
Uopt_qobj = qt.Qobj(U_opt, dims=[[2, 2, 2], [2, 2, 2]]) # create QuTip quantum object
# create 2×2 subplot layout
fig, axes = plt.subplots(2, 2, figsize=(10, 10), subplot_kw={"projection": "3d"})
# define titles and data
titles = [
"Ideal unitary (real)",
"Ideal unitary (imaginary)",
"Optimised unitary (real)",
"Optimised unitary (imaginary)",
]
qobjs = [Uideal_qopj, Uideal_qopj, Uopt_qobj, Uopt_qobj]
bar_styles = ["real", "img", "real", "img"]
# prepare each subplot
for ax, title, qobj, style in zip(axes.flat, titles, qobjs, bar_styles):
qt.matrix_histogram(
qobj,
x_basis=range(8),
y_basis=range(8),
bar_style=style,
ax=ax,
limits=[-1, 1],
colorbar=False,
)
ax.set_title(title, fontsize=10)
ax.collections[0].set_cmap("viridis") # Correct way to set colormap
ax.set_zticks([-1, 0, 1])
ax.set_xticklabels([f"{i:03b}" for i in range(8)], fontsize=4)
ax.set_yticklabels([f"{i:03b}" for i in range(8)], fontsize=4)
ax.tick_params(axis="both", labelsize=8)
# adjust layout and show plot
plt.subplots_adjust(wspace=0.3, hspace=0.3)
plt.show()
As you can see, our final unitary shows a much better match with the ideal CZ gate, up to a global phase. Our optimisation was successful!
We can quantify this success by printing the fidelity before and after optimisation. (Remember: we defined our loss function as the infidelity = 1 - fidelity).
# Calculate final fidelities
fidelity_initial = 1 - loss(U0, Ut)
fidelity_optimised = 1 - loss(U_opt, Ut)
print(f"Initial fidelity: {fidelity_initial*100:.5f}%")
print(f"Optimised fidelity: {fidelity_optimised*100:.5f}%")
Initial fidelity: 56.25000% Optimised fidelity: 99.99995%
We can clearly see that our optimisation has significantly enhanced the fidelity!
Congratulations, you’ve just performed optimal control on a two-qubit CZ gate! You can now use these methods to start optimising your QPU performance.