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 the Q-CTRL Python package is to represent quantum systems. The graph representation of such a system can be used to perform several tasks, including optimal 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).

Summary workflow for graph-based representations of quantum systems

1. Create a graph object

The first step is to use the qctrl.create_graph function to create the Python object representing the graph. This graph will contain all the nodes that you create when describing the quantum system. The graph object needs to be active (by using Python with syntax) whenever you create new nodes.

2. Define graph nodes representing your quantum system

The next step is to add nodes to the graph representing your quantum system. To enable convenient and efficient computations, Q-CTRL offers a wide variety of nodes specific to quantum systems, in addition to the standard mathematical operations common to other data-flow programming libraries. All nodes are available in the qctrl.operations namespace of the Q-CTRL Python package.

Typically you will first create nodes that describe the Hamiltonian of your quantum system. A standard flow is:

  • Create the inputs to your graph. 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.
  • Create "signals", or 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.
  • Create "operators", or operator-valued functions of time, by multiplying constant matrices (for example Pauli matrices) by signals.
  • Sum the operators into a single Hamiltonian operator.

Note that while this approach of 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, the Q-CTRL Python package 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).

General mathematical rules for assembling graph nodes

Creating nodes and applying operations

You create graph nodes by calling functions in the qctrl.operations namespace. There are three basic types of nodes:

  • Tensors, representing multidimensional arrays of data,
  • PWCs, representing piecewise-constant functions of time, and
  • STFs, representing sampleable time-dependent functions.

For basic arithmetic (for example addition and multiplication), you can use standard Python syntax, which is a convenient short-hand for the corresponding qctrl.operations nodes.

Creating signals and operators

Signals are scalar functions: they only contain one numerical value in each period of time, which means their shape is zero-dimensional. You can create a PWC signal using qctrl.operations.pwc_signal operation.

You can combine signals with constant matrices to produce operators, which represent piecewise-constant operator-valued (2D) functions of time. Usually these operators represent individual terms in your Hamiltonian. You can also create constant operators, representing static terms in your Hamiltonian.

Using operations on multi-dimensional objects

You can directly sum PWC operators to yield the overall system Hamiltonian, as long as all the operators have compatible shapes. In general, you can apply basic elementwise arithmetic operations between arrays, Tensors, PWCs, or STFs such as +, -, *, /, // (floor division), and ^ (exponentiation).

Note that operations between PWCs/STFs and NumPy arrays or Tensors correspond to applying the operation to every value that the PWC/STF assumes in time, while operations between two PWCs/STFs correspond to the operation between the values that the objects assume in each time window. Operations between a PWC and an STF are not allowed.

If the two objects don't have the same shape, the Q-CTRL Python package attempts to broadcast them following the same broadcasting rules of the NumPy package. They consist of these two steps:

  1. If one of the objects has fewer dimensions than the other, more dimensions with size 1 are added to the front until both have the same number of dimensions.
  2. If one of the dimensions has size 1 for one of the objects but not for the other, the length of this dimension is increased to the largest value between the two.

Given these rules, it is not possible to broadcast two objects if they have different sizes in one of their dimensions and neither of those sizes is 1. For example, it is not possible to broadcast an operator of shape $2\times 2$ with an operator of shape $3\times 3$ (they have different sizes in both dimensions but none of them is 1), but it is possible to broadcast an operator of shape $2\times 1$ with an operator of shape $1\times 2$.

Besides basic arithmetic, you can also calculate other operations between objects, such as matrix multiplication (using the @ operator or qctrl.operations.matmul), the Kronecker product (with qctrl.operations.kron), or create a complex object from its real and imaginary parts (via qctrl.operations.complex_value). However, in these cases the shape compatibility depends on the particular operation.

The Q-CTRL Python package also includes many operations that you can apply to Tensors, PWCs, or STFs. These include basic mathematic functions (like sqrt, exp, log), trigonometric functions (sin, cosh, arctan), functions for complex objects (real, angle, conjugate), or even matrix functions (trace, adjoint). You can apply them by just calling the relevant node, for instance, qctrl.operations.sqrt. You can search the full list of available operations to learn more about them.

3. Define graph nodes representing computations on your quantum system

Once you have defined nodes that describe your quantum system, you can add nodes representing the computation to perform on it. The Q-CTRL Python package provides several functions for performing computations on quantum systems, for example for calculating the time evolution operators for the system.

4. Add names to the output nodes for your computation

You can give names to any nodes whose values you want to extract from the graph. For example, for a simulation you likely want to give names to the nodes representing the time evolution operators or evolved states of your system, while for an optimization you likely want to give names to the nodes representing the control pulses that create the optimized gates.

5. Run your computation

The final step is to use a function in the qctrl.functions namespace to run your computation and get the values of the requested output nodes.

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 $\omega_k,\omega_c$ are the frequencies of the qubits, $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$, $1/\Delta=(1/\Delta_a+1/\Delta_b)/2$.

By adjusting the frequency $\omega_c$ of the coupling qubit, the effective coupling strength $\tilde g$ can be tuned.

This example shows how you can represent this system as a graph, and calculate its time evolution operators for a particular choice of the $\omega_c$ control.

Defining system parameters

First, define some basic parameters for the system.

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

from qctrl import Qctrl

# Starting a session with the API
qctrl = Qctrl()
# Standard qubit operators
identity = np.array([[1.0, 0.0], [0.0, 1.0]], dtype=complex)
sigma_z = np.array([[1.0, 0.0], [0.0, -1.0]], dtype=complex)
sigma_minus = np.array([[0.0, 1.0], [0.0, 0.0]], dtype=complex)
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 = np.kron(sigma_z, identity)
sigma_z_b = np.kron(identity, sigma_z)
sigma_plus_minus = np.kron(sigma_minus.T, sigma_minus)

# 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

Building the Hamiltonian

Creating the graph object

The graph object obtained from the qctrl.create_graph function contains all the nodes you create. The object is a Python context manager, which you can activate using the with keyword, and which needs to be active any time you create a node using the qctrl.operations functions.

graph = qctrl.create_graph()

Creating the inputs

This example sets up a simulation with inputs that are known values for control pulses. Consider a simple piecewise-constant pulse for $\omega_c$ which is then passed through imperfect control lines. The effect of these control lines is to distort the signal that actually reaches the system, and the values from this distorted signal are the input to function that creates 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.

# Enter the graph context, to ensure that nodes are added to the correct graph
with graph:

    # 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 = qctrl.operations.pwc_signal(
        duration=gate_duration, values=omega_c_raw_values, name="omega_c_raw"

    # Apply a sinc filter to simulate the effect of control line imperfections
    omega_c_smooth = qctrl.operations.convolve_pwc(
        kernel=qctrl.operations.sinc_convolution_kernel(2 * np.pi * 0.1e9),

    # Discretize the resulting smooth signal to obtain a finely-sampled
    # piecewise-constant approximation
    omega_c = qctrl.operations.discretize_stf(

Creating signals

The omega_c object in the previous code block is a piecewise-constant (PWC) scalar function of time, also called a PWC signal in this context. The following code block shows how to manipulate this signal, which represents $\omega_c$, to create PWC signals representing the derived quantities $g_a^2/(2\Delta_a)$, $g_b^2/(2\Delta_b)$, and $\tilde g$.

with graph:
    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 * qctrl.operations.sqrt(omega_a * omega_c / (C_a * C_c))
    g_b = 0.5 * C_bc * qctrl.operations.sqrt(omega_b * omega_c / (C_b * C_c))
    g_tilde = (
        * (omega_c * eta / (2 * Delta) - omega_c * eta / (2 * Sigma) + eta + 1)
        * C_ab
        * qctrl.operations.sqrt(omega_a * omega_b / (C_a * C_b))

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

Creating operators

You can combine the signals that you created in the previous section with constant matrices to produce operators, which represent PWC operator-valued (2D) functions of time. Usually these operators represent individual terms in your Hamiltonian. You can also create constant operators, representing static terms in your Hamiltonian—in this case the fixed frequency terms for each qubit.

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

    detuning_b_term = detuning_b * sigma_z_b
    coupling_term = g_tilde * (sigma_plus_minus + sigma_plus_minus.T)

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

    # You can also explicity declare constant PWC operators.
    fixed_frequency_b_term = qctrl.operations.constant_pwc_operator(
        operator=omega_b * sigma_z_b / 2, duration=gate_duration

    # You can add operators.
    hamiltonian = (
        + detuning_b_term
        + coupling_term
        + fixed_frequency_a_term
        + fixed_frequency_b_term

Describing the computation to perform on the system

You can then add nodes that perform computations on your Hamiltonian. The code below shows how to compute the time evolution operators for the system, and gives them a name so that they can be fetched from the computation.

with graph:
    time_evolution_operators = qctrl.operations.time_evolution_operators_pwc(
        sample_times=np.linspace(0, gate_duration, 10),

Running the computation

Now that the graph is defined, you can use the qctrl.functions.calculate_graph function to evaluate the time evolution operators for the system as well as the raw and filtered input signals that you defined earlier (and, more generally, can be used to calculate the values of any named nodes in the graph).

result = qctrl.functions.calculate_graph(
    output_node_names=["time_evolution_operators", "omega_c_raw", "omega_c"],
Your task calculate_graph has completed.

Extracting the results

You can extract the evaluated values from the result object. The specific form of outputs is discussed in detail in the more specific user guides, so this example simply displays a visualization of the raw and filtered input signals (using the Q-CTRL Visualizer Python package), and prints the final time evolution operator. You can see that the chosen controls produce a gate that is, up to single-qubit phases, close to an iSWAP.

        "$\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 example showed the general procedure for describing a quantum system, and a computation performed on that quantum system, as a graph object. You can visit the other User guides and Application notes to see several more approaches to using graphs for describing different types of systems and solving different types of problems, and the reference documentation to see more details about all functions and types in the Q-CTRL Python package.