# How to represent quantum systems using graphs

Represent quantum systems for optimization, simulation, and other tasks using graphs

The typical purpose of graph objects in Boulder Opal is to represent quantum systems. The graph representation of such a system can be used to perform several tasks, including robust control (for calculating optimized control pulses), simulation (to understand the dynamics of the system in the presence of specific controls and noises), and system identification (to estimate the values of unknown system parameters based on measurements of the system).

Please refer to our topic Understanding graphs in Boulder Opal for context on what graphs are used for and why.

In what follows, we use graphs to define time-dependent Hamiltonians for simulation or optimization.

## Summary workflow

Here we outline the general procedure for representing your quantum systems and the computations performed on them in the form of graph objects.

### 1. Create graph inputs

For an optimization, these are optimizable variables that can be tuned by the optimizer to minimize a cost function. For a simulation, the inputs can be known values for control pulses, or quantities derived from known values.

### 2. Create signals

Signals are scalar-valued functions of time, which may be non-linear, have enforced temporal structure such as time symmetry, or more generally can depend arbitrarily on the inputs. They only contain one numerical value in each period of time, which means their shape is zero-dimensional. Signals represent the time-dependent envelope of the Hamiltonian. You create piecewise-constant (PWC) signals using the graph.pwc_signal operation.

### 3. Create Hamiltonian operators

You create operators or PWC operator-valued (2D) functions of time by multiplying constant matrices (for example Pauli matrices) with signals. Usually these operators represent individual terms in your Hamiltonian. You can also create constant operators to represent static terms in your Hamiltonian. Finally, you sum the individual operators created into a single Hamiltonian operator.

Once you have defined the nodes that describe your quantum system, you can add extra nodes representing the computations that you want to perform on it. For example, Boulder Opal offers time evolution operations, among many others.

Note that while this approach to constructing Hamiltonians is the most common, it is not a requirement. You can use graphs to perform a wide variety of other computations too. For example, Boulder Opal also provides specialized functions for working with trapped ions systems that take advantage of certain approximations to bypass Hamiltonian-level descriptions of the system (see the How to design error-robust Mølmer–Sørensen gates for trapped ions user guide for details).

## Worked example: Two-qubit system with a tunable coupling qubit

Consider a system consisting of three transmon qubits, in which the central qubit acts as a tunable coupler between the other two, as proposed by Yan et al. This system can be approximated by the following two-qubit Hamiltonian: $$\tilde H = \frac{\omega_a}{2} \sigma_{z, a} +\frac{\omega_b}{2} \sigma_{z, b} + \frac{g_a^2}{2 \Delta_a} \sigma_{z,a}+ \frac{g_b^2}{2 \Delta_b} \sigma_{z,b} + \tilde g \left(\sigma_{+, a} \sigma_{-,b} + \sigma_{-,a} \sigma_{+,b} \right) ,$$ where $\omega_k$ are the qubit frequencies, $\Delta_k=\omega_k-\omega_c$ are the detunings of each qubit from the coupling qubit, $g_k$ are the direct couplings between each qubit and the coupling qubit, and $\tilde g$ is the effective coupling between the two qubits.

The couplings $g_k,\tilde g$ are given by $$g_k \approx \frac{1}{2}\frac{C_{kc}}{\sqrt{C_kC_c}}\sqrt{\omega_k\omega_c},$$ $$\tilde g = \frac{1}{2}\left[\frac{\omega_c}{2\Delta}\eta - \frac{\omega_c}{2\Sigma}\eta + \eta + 1\right]\frac{C_{ab}}{\sqrt{C_aC_b}}\sqrt{\omega_a\omega_b} ,$$ where $C_k,C_c$ are the capacitances of the qubits, $C_{xy}$ are the qubit-qubit capacitances between each pair of qubits, and the derived quantities are given by $\eta=C_{ac}C_{bc}/C_{ab}C_c$, $1/\Sigma=(1/\Sigma_a+1/\Sigma_b)/2$, $\Sigma_k=\omega_k+\omega_c$, and $1/\Delta=(1/\Delta_a+1/\Delta_b)/2$. We can see that the effective coupling strength $\tilde g$ can be tuned by adjusting the frequency $\omega_c$ of the coupling qubit.

In this example, we will show how you can represent this system as a graph, and calculate its time evolution operators for a particular choice of the $\omega_c$ control. The chosen controls will produce a gate that is, up to single-qubit phases, close to an iSWAP.

We will consider a simple PWC pulse for $\omega_c$ which is then passed through imperfect control lines. These control lines distort the signal that actually reaches the system, and thus the relevant term in the Hamiltonian. See the How to add smoothing and band-limits to optimized controls user guide for more details about incorporating filters into your graphs.

From this filtered $\omega_c$ signal, we will also create signals representing the different terms in the Hamiltonian, such as $g_a^2/(2\Delta_a)$, $g_b^2/(2\Delta_b)$, and $\tilde g$. We will also add constant operators to the Hamiltonian, representing its static terms—in this case the fixed frequency terms for each qubit.

import matplotlib.pyplot as plt
import numpy as np
from qctrlvisualizer import plot_controls

from qctrl import Qctrl

# Start a Boulder Opal session.
qctrl = Qctrl()

# Define system parameters.

# Physical parameters.
C_a = 70  # fF
C_b = 72  # fF
C_ac = 4  # fF
C_bc = 4.2  # fF
C_c = 200  # fF
C_ab = 0.1  # fF
omega_a = 2 * np.pi * 4e9  # Hz
omega_b = 2 * np.pi * 4e9  # Hz
eta = C_ac * C_bc / (C_ab * C_c)
gate_duration = 100e-9  # s

# Create graph.
graph = qctrl.create_graph()

# Standard qubit operators.
iswap = np.array(
[
[1.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 1j, 0.0],
[0.0, 1j, 0.0, 0.0],
[0.0, 0.0, 0.0, 1.0],
]
)

# 2-qubit operators.
sigma_z_a = graph.pauli_kronecker_product([("Z", 0)], 2)
sigma_z_b = graph.pauli_kronecker_product([("Z", 1)], 2)
sigma_plus_minus = graph.pauli_kronecker_product([("P", 0), ("M", 1)], 2)

# Create basic piecewise-constant signal.
omega_c_raw_values = (
(5.4 - 0.34 * np.array([0, 1, 2, 3, 2, 1, 0])) * 2 * np.pi * 1e9
)  # Hz
omega_c_raw = graph.pwc_signal(
duration=gate_duration, values=omega_c_raw_values, name="omega_c_raw"
)

# Apply a sinc filter (and rediscretize) to simulate
# the effect of control line imperfections.
omega_c = graph.utils.filter_and_resample_pwc(
pwc=omega_c_raw,
cutoff_frequency=2 * np.pi * 0.1e9,
segment_count=100,
name="omega_c",
)

# Calculate derived couplings in the Hamiltonian terms.
Delta_a = omega_a - omega_c
Delta_b = omega_b - omega_c
Delta = 2 / (1 / Delta_a + 1 / Delta_b)
Sigma_a = omega_a + omega_c
Sigma_b = omega_b + omega_c
Sigma = 2 / (1 / Sigma_a + 1 / Sigma_b)

g_a = 0.5 * C_ac * graph.sqrt(omega_a * omega_c / (C_a * C_c))
g_b = 0.5 * C_bc * graph.sqrt(omega_b * omega_c / (C_b * C_c))
g_tilde = (
0.5
* (omega_c * eta / (2 * Delta) - omega_c * eta / (2 * Sigma) + eta + 1)
* C_ab
* graph.sqrt(omega_a * omega_b / (C_a * C_b))
)

detuning_a = g_a**2 / (2 * Delta_a)
detuning_b = g_b**2 / (2 * Delta_b)

# Build the system's Hamiltonian from its terms.

detuning_a_term = detuning_a * sigma_z_a
# This is equivalent to using graph.pwc_operator:
# detuning_a_term = graph.pwc_operator(detuning_a, sigma_z_a).

detuning_b_term = detuning_b * sigma_z_b
coupling_term = g_tilde * 2.0 * graph.hermitian_part(sigma_plus_minus)

# Tensors and arrays are considered constant operators.
fixed_frequency_a_term = omega_a * sigma_z_a / 2

# We can also explicity declare constant PWC operators.
fixed_frequency_b_term = graph.constant_pwc_operator(
operator=omega_b * sigma_z_b / 2, duration=gate_duration
)

# Add all operators into full Hamiltonian.
hamiltonian = (
detuning_a_term
+ detuning_b_term
+ coupling_term
+ fixed_frequency_a_term
+ fixed_frequency_b_term
)

# Compute the time evolution operators for the system.
# We provide a name to extract them after you execute the graph.
time_evolution_operators = graph.time_evolution_operators_pwc(
hamiltonian=hamiltonian,
sample_times=np.linspace(0, gate_duration, 10),
name="time_evolution_operators",
)

# Evaluate the time evolution operators for the system, and the raw and filtered input signals defined earlier.
result = qctrl.functions.calculate_graph(
graph=graph,
output_node_names=["time_evolution_operators", "omega_c_raw", "omega_c"],
)

Your task calculate_graph (action_id="1164458") has completed.

plot_controls(
{
"$\omega_c^\mathrm{raw}$": result.output["omega_c_raw"],
"$\omega_c^\mathrm{filtered}$": result.output["omega_c"],
}
)

print("Final time evolution operator:")
print(np.round(result.output["time_evolution_operators"]["value"][-1], 3))


Final time evolution operator:
[[-0.126-0.992j  0.   +0.j     0.   +0.j     0.   +0.j   ]
[ 0.   +0.j    -0.169-0.058j  0.   +0.984j  0.   +0.j   ]
[ 0.   +0.j    -0.   +0.984j -0.169+0.058j  0.   +0.j   ]
[ 0.   +0.j     0.   +0.j     0.   +0.j    -0.126+0.992j]]


This notebook was run using the following package versions. It should also be compatible with newer versions of the Q-CTRL Python package.

Package version
Python 3.9.12
matplotlib 3.5.1
numpy 1.21.5
scipy 1.7.3
qctrl 19.1.0
qctrlcommons 17.1.1
qctrltoolkit 1.5.0
qctrlvisualizer 3.2.1