# Set up quantum systems

**Building objects that describe your quantum system**

In order to use the Q-CTRL Python package to analyze and improve your quantum system, you must first translate your system into the language used by the Q-CTRL ecosystem. The Q-CTRL Python package uses two different representations of quantum systems:

*Graphs*, a general and flexible representation that can be used for any quantum system.*Drives, shifts, and drifts*, a more restrictive representation that forces the relationships between control signals and Hamiltonian terms to be linear.

Eventually all features will use graphs, but some features currently still require drives, shifts, and drifts.

This notebook demonstrates how to describe quantum systems using both of these representations.

```
import attr
import matplotlib.pyplot as plt
import numpy as np
from qctrlopencontrols import new_gaussian_control
from qctrlvisualizer import plot_controls, get_qctrl_style
from qctrl import Qctrl
# Starting a session with the API
qctrl = Qctrl()
plt.style.use(get_qctrl_style())
```

## Graphs

### Graph-based computation

Several functions in the Q-CTRL Python package use *graphs* to represent quantum systems. Graphs use collections of *nodes* and *edges* to represent computations that map inputs to outputs. Formally, these graphs are *directed* and *acyclic*, and are commonly known as data-flow graphs, but we refer to them simply as graphs. Each node in a graph describes a primitive computation performed on its input(s), for example the exponential of a number or the multiplication of two matrices. Different nodes are joined by directed edges, which describe how the output of one node is connected to the input of another node. By joining different types of nodes in different ways, you can build graphs representing computations as simple as basic arithmetic through to computations as complex as state propagation in high-dimensional quantum systems. The nodes provided by Q-CTRL are all available in the `operations`

namespace of the Q-CTRL Python package.

The graph representation offers three key advantages over other approaches to representing generic computations:

- Flexibility: no particular structure is enforced, so you can represent
*any*computation that can be expressed in terms of the provided nodes. - Efficiency: graphs can be evaluated extremely efficiently.
- Automatic differentiation: graphs can be automatically differentiated, enabling features like gradient-based optimization and calculation of Hessian matrices.

Before explaining how graphs are used to represent quantum systems, the next sections present a brief introduction to working with them in the abstract sense.

#### Creating graphs

The first step in creating a graph is to use the `qctrl.create_graph`

function to create the Python object representing the graph. Technically this object is a *Python context manager*, which you enter using the `with`

keyword. With the graph object created, you can then create graph nodes by calling functions in the `qctrl.operations`

namespace. All calls to `qctrl.operations`

functions must be made inside the context manager that you created, to ensure that nodes are added to the correct graph. Below we show how to create a simple graph that adds two numbers.

```
with qctrl.create_graph() as graph:
result = qctrl.operations.add(1.5, 2.5)
```

```
print(result)
```

The key point to understand is that the graph itself does not perform the computation; instead it is a *recipe for performing the computation remotely*. You can only get the computed values out of a graph after you have evaluated it, which you will see how to do later in this guide.

For now, we can explore some more consequences of the fact that graphs do not get evaluated immediately. Below, we add a node to the graph that represents an *optimization variable*. You'll see more of optimization variables later, but for now you only need to know that the `qctrl.operations.optimization_variable`

function returns an object representing `count`

(in this case 10) values.

```
with graph:
variables = qctrl.operations.optimization_variable(
count=10,
lower_bound=0,
upper_bound=1,
)
print(variables)
```

As with the addition in the earlier code block, the `variables`

object is not an array of numbers, but it is a *representation* of an array of numbers that will be computed in the cloud when you evaluate the graph. We call such an object a *tensor*. You can manipulate tensors in much the same way as NumPy arrays, although you need to use `qctrl.operations`

functions instead of NumPy functions. For certain basic arithmetic operations you can also use regular Python syntax, which is a convenient wrapper around the corresponding `qctrl.operations`

functions. Finally, when calling a `qctrl.operations`

function that accepts tensors, you can usually also pass NumPy arrays.

```
with graph:
first_variable = variables[0]
print(first_variable)
scaled_variables = variables * 5
print(scaled_variables)
added_variables = qctrl.operations.add(scaled_variables, np.linspace(0, 1, 10))
# This is equivalent to the following code:
# added_variables = scaled_variables + np.linspace(0, 1, 10)
print(added_variables)
multiplied_variables = scaled_variables * variables
print(multiplied_variables)
```

Tensors are not the only types that live in graphs. Some `qctrl.operations`

functions return other types of data too, representing (for example) piecewise-constant functions of time (PWCs) or sampleable tensor-valued functions of time (STFs). While these other types do not represent simple arrays of numbers, they are still similar to tensors in the sense that they represent the result of a remote computation that will be performed in the cloud. You will see examples of these types later in this guide, and you can also see the reference documentation for more details.

#### Evaluating graphs

There are currently two ways to evaluate graphs:

- Using
`calculate_graph`

, which simply calculates and returns the values of specific nodes in the graph. - Using
`calculate_optimization`

, which optimizes special*optimization variable*nodes in order to minimize a specified*cost*node, and then returns the values of specific nodes in the graph. All graphs used for optimizations contain at least one of the special optimization variable nodes (and such nodes cannot be used with`calculate_graph`

).

You will see both of these functions used later in this guide.

### Quantum systems as 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).

To enable convenient and efficient computations, Q-CTRL offers a wide variety of nodes to help build graphs describing quantum systems, in addition to the standard mathematical operations common to other data-flow programming libraries. As with the basic nodes introduced earlier, the nodes for working with quantum systems are available in the `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 "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 input parameters.
- 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.

From this Hamiltonian you can then add nodes to perform further computations, for example calculating time evolution operators, operational infidelities, or filter function values.

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 Mølmer–Sørensen User guide for details).

Full descriptions of all the nodes provided by the Q-CTRL Python package are available in the reference documentation.

### 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.

#### Defining system parameters

First, define some basic parameters for the system.

```
# 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],
],
dtype=complex,
)
# 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

The next code blocks define a function that creates a graph object for the Hamiltonian of the tunable coupler system. The inputs to the function are the parameters that you may want to vary (which could be for many purposes, for example optimizable parameters in optimal control, input parameters for simulation, or unknown parameters for system identification). In this case, the function accepts a piecewise-constant (PWC) function of time for the coupling frequency $\omega_c$, and returns the corresponding PWC function of time for the Hamiltonian. Recall that a PWC object simply refers to a multi-dimensional array of real or complex numbers (and their associated time durations) that will be calculated remotely when you evaluate the graph.

In order to explain each step, each code block shows just one part of the function. In your own code, it would be simpler to merge these functions into one.

##### Using operations on scalar objects

As described earlier, you create graph nodes by calling functions in the `qctrl.operations`

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

functions. You can also use the shorthand `@`

for matrix multiplication, as long as the objects involved have at least two dimensions.

The following code block shows how to use these features to define the derived quantities $g_a^2/(2\Delta_a)$, $g_b^2/(2\Delta_b)$, and $\tilde g$, all of which depend on the piecewise-constant function of time (PWC) $\omega_c$. Note that operations between PWCs and NumPy arrays or Tensors correspond to applying the operation to every value that the PWC assumes in time, while operations between two PWCs correspond to the operation between the values that the objects assume in each time window. In this example, this means that all the returned quantities derived from the PWC $\omega_c$ are also PWCs.

```
def signals(omega_c):
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 = (
0.5
* (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)
return detuning_a, detuning_b, g_tilde
```

##### Creating signals and operators

The PWCs returned in the previous code block are scalar functions: they only contain one numerical value in each period of time, which means their shape is zero-dimensional. Scalar functions of time are also called *signals* in this context, and you can create them using the `pwc_signal`

operation.

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.

```
def operators(omega_c):
detuning_a, detuning_b, coupling = signals(omega_c)
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 = coupling * (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
)
return (
detuning_a_term,
detuning_b_term,
coupling_term,
fixed_frequency_a_term,
fixed_frequency_b_term,
)
```

##### 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).

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:

- 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. - 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.

```
def hamiltonian(omega_c):
(
detuning_a_term,
detuning_b_term,
coupling_term,
fixed_frequency_a_term,
fixed_frequency_b_term,
) = operators(omega_c)
return (
detuning_a_term
+ detuning_b_term
+ coupling_term
+ fixed_frequency_a_term
+ fixed_frequency_b_term
)
```

#### Building a graph for a specific computation

The `hamiltonian`

function defined above creates a graph object representing the Hamiltonian of the system. In order to do computations (for example simulations or optimizations), you embed this graph inside a larger graph that represents the overall computation. The next sections show how this procedure works for a simulation.

##### Creating the 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.

```
graph = qctrl.create_graph()
```

##### Creating the Hamiltonian from the input values

The first part of the graph sets up the necessary inputs to the Hamiltonian. For an optimization, these inputs typically depend on optimizable variables. For a simulation, the inputs can be known values for control pulses, or quantities derived from known values.

This example demonstrates the latter. 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 Optimize controls User guide for more details about incorporating filters into your graphs.

Note that you can give a name to any Tensor or PWC produced from a node, to make it easier to extract its value after the computation has been performed. The next code block gives names to the raw and filtered signals, so that you can later visualize the effect of the control lines.

```
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(
pwc=omega_c_raw,
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_filtered = qctrl.operations.discretize_stf(
stf=omega_c_smooth,
duration=gate_duration,
segment_count=100,
name="omega_c_filtered",
)
# Create the Hamiltonian using the filtered signal as the input
h = hamiltonian(omega_c_filtered)
```

##### Creating derived properties

The Q-CTRL Python package provides several functions for deriving properties of quantum systems from Hamiltonians. In the following example, the operation calculates the time evolution operators for the system (see the Simulation User guide for more details).

```
with graph:
time_evolution_operators = qctrl.operations.time_evolution_operators_pwc(
hamiltonian=h,
sample_times=np.linspace(0, gate_duration, 10),
name="time_evolution_operators",
)
```

#### Calling a function with the graph

The graph constructed represents both the quantum system and the desired computation on that quantum system. In this case the computation is a simple simulation, but you can create similar graphs for both optimization and system identification (see the Optimization User guide for examples). To run the simulation, use the `calculate_graph`

function, which evaluates 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(
graph=graph,
output_node_names=["time_evolution_operators", "omega_c_raw", "omega_c_filtered"],
)
```

#### 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 print 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.

```
plot_controls(
plt.figure(),
{
"$\omega_c^\\mathrm{raw}$": result.output["omega_c_raw"],
"$\omega_c^\\mathrm{filtered}$": result.output["omega_c_filtered"],
},
)
plt.show()
print("Final time evolution operator:")
print(np.round(result.output["time_evolution_operators"]["value"][-1], 3))
```

#### Summary

This example showed the general procedure for describing a quantum system, and a computation performed on that quantum system, as a graph object. This approach of defining the Hamiltonian inside a function and then plugging this function into the graph creation is simple and powerful, but not the only possible approach. For example, for simple systems you might find it easier to define your entire graph within a single block, while for for complex systems it might be helpful to split creation of the graph across several different functions. 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.

### Example: One-qubit system for an optimization

Next we present a full example of a simple one-qubit system set up for optimization. We consider the Hamiltonian:

\begin{align*} H(t) = & \frac{\nu}{2} \sigma_z + \frac{\Omega(t)}{2} \sigma_- + \frac{\Omega^*(t)}{2} \sigma_+ + \frac{\Delta(t)}{2} \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, and $\sigma_k$ are the Pauli matrices.

We will use optimization to find controls $\Omega(t)$ and $\Delta(t)$ that produce a NOT gate. As above, we use a graph to represent the computation. This time, the inputs to the Hamiltonian are *optimization variables*, that represent quantities that can be controlled and optimized. Once the Hamiltonian is constructed, it is used as the input to an *infidelity* node, which represents the cost function minimized by the optimizer in order to find the optimal controls. See the Optimization User guide for details and more examples.

```
# Define standard matrices
sigma_x = np.array([[0, 1], [1, 0]])
sigma_z = np.array([[1, 0], [0, -1]])
sigma_m = np.array([[0, 1], [0, 0]])
# Define physical constants
nu = 2 * np.pi * 0.5 * 1e6 # Hz
omega_max = 2 * np.pi * 0.5e6 # Hz
delta_max = 2 * np.pi * 0.5e6 # Hz
segment_count = 50
duration = 10e-6 # s
# Create the graph describing the system
with qctrl.create_graph() as graph:
# Create the optimizable variables for the controls
omega_moduli = qctrl.operations.optimization_variable(
count=segment_count,
lower_bound=0,
upper_bound=omega_max,
)
omega_phases = qctrl.operations.optimization_variable(
count=segment_count,
lower_bound=0,
upper_bound=2 * np.pi,
is_lower_unbounded=True,
is_upper_unbounded=True,
)
delta_values = qctrl.operations.optimization_variable(
count=segment_count,
lower_bound=-delta_max,
upper_bound=delta_max,
)
# Create the Hamiltonian
# Create the signals
omega = qctrl.operations.complex_pwc_signal(
moduli=omega_moduli,
phases=omega_phases,
duration=duration,
name="omega",
)
delta = qctrl.operations.pwc_signal(
values=delta_values,
duration=duration,
name="delta",
)
# Create the operators
detuning = nu * sigma_z / 2
drive = qctrl.operations.pwc_operator_hermitian_part(omega * sigma_m)
clock_shift = delta * sigma_z / 2
# Sum the operators
hamiltonian = detuning + drive + clock_shift
# Calculate the infidelity, derived from the Hamiltonian
infidelity = qctrl.operations.infidelity_pwc(
hamiltonian=hamiltonian,
target_operator=qctrl.operations.target(sigma_x),
name="infidelity",
)
# Call the function
result = qctrl.functions.calculate_optimization(
graph=graph,
cost_node_name="infidelity",
output_node_names=["omega", "delta"],
optimization_count=1,
)
# Visualize the results
plot_controls(
plt.figure(),
{"$\Omega(t)$": result.output["omega"], "$\Delta(t)$": result.output["delta"]},
)
print("Infidelity:")
print(result.cost)
```

```
```

### Example: Detuning scan of a one-qubit system

Next we demonstrate to how to set up a one-qubit system for parallel simulations with different parameters, that is, a batch of simulations. We consider the Hamiltonian:

\begin{align*} H(t) = & \frac{\nu}{2} \sigma_z + \frac{\Omega(t)}{2} \sigma_- + \frac{\Omega^*(t)}{2} \sigma_+, \end{align*}where $\nu$ is the qubit detuning, $\Omega(t)$ is a time-dependent Rabi rate, and $\sigma_k$ are the Pauli matrices. We will use a set of different qubit detuning values to create a batch of Hamiltonians and a Gaussian shaped time-dependent Rabi rate that produces a NOT gate. As in the previous example, we will use a graph to represent the computation, which allows the use of batches in most cases and performs their computation in a single operation.

Batches are created by adding leading dimensions to constants, signals and operators in your system. In this example, the constant detuning is extended to an array of shape (101,) to represent a batch of 101 detuning values. In order to disambiguate the one dimensional batch of constants from an array of signal values you have to specify the number of batch dimensions when creating a PWC or an STF object.

General broadcasting rules apply when working with batches, that is, two batches of shape ($b_1$, 1) and ($b_2$,) lead to a batch of shape ($b_1$, $b_2$). In this example, you add a batch of Hamiltonian terms (the detuning) to a single Hamiltonian term (the drive), and the result is a batch of Hamiltonian terms. Once the Hamiltonian is constructed, it is used as the input to an *infidelity* node, which represents a batch of cost functions, with the same batch dimension as the Hamiltonian.

```
# Define standard matrices
sigma_x = np.array([[0, 1], [1, 0]])
sigma_z = np.array([[1, 0], [0, -1]])
sigma_m = np.array([[0, 1], [0, 0]])
# Define physical constants
duration = 10e-6 # s
segment_count = 50
# Define the detuning values
detuning_values = 2 * np.pi * np.linspace(-5, 5, 101) * 1e4
# The values of the normalized Gaussian
omega_values = new_gaussian_control(
rabi_rotation=np.pi,
duration=duration,
segment_count=segment_count,
width=duration / 2,
).rabi_rates
# Create the graph describing the system
with qctrl.create_graph() as graph:
# Create the Hamiltonian
# Create the signals
# The detuning is a one dimensional batch of constants
delta = qctrl.operations.constant_pwc(
constant=detuning_values, duration=duration, batch_dimension_count=1
)
omega = qctrl.operations.pwc_signal(values=omega_values, duration=duration)
# Create the operators
# The detuning is a batch of shape (101,) of PWC operators of shape (2, 2)
detuning = delta * sigma_z / 2
# The drive is a PWC operator of shape (2, 2)
drive = qctrl.operations.pwc_operator_hermitian_part(omega * sigma_m)
# The sum is a batch of shape (101,) of PWC operators of shape (2, 2)
hamiltonian = detuning + drive
# Calculate the infidelities, derived from the Hamiltonian
infidelities = qctrl.operations.infidelity_pwc(
hamiltonian=hamiltonian,
target_operator=qctrl.operations.target(sigma_x),
name="infidelities",
)
# Call the function
detuning_result = qctrl.functions.calculate_graph(
graph=graph, output_node_names=["infidelities"]
)
# Visualize the results
plt.plot(
detuning_values / (2 * np.pi * 1e6),
detuning_result.output["infidelities"]["value"],
)
plt.xlabel("Frequency drift (MHz)")
plt.ylabel("Infidelity")
plt.title("Detuning scan")
plt.show()
```

### Example: Optimizing a Gaussian envelope function applied to fast oscillating controls

You can also optimize signals whose components relate to each other in more complicated ways. For example, consider that you want to add additional restrictions to the optimal controls $\Omega(t)$ that you create when optimizing an X gate for a single-qubit system whose Hamiltonian has the form:

$$ H(t) = \frac{\Omega(t)}{2} \sigma_x. $$In this case, suppose that you want the controls $\Omega(t)$ to be composed of an oscillating carrier signal $c(t)$ multiplied by a slowly-varying Gaussian envelope $E(t)$:

$$ \Omega(t) = c(t) E(t). $$The smooth carrier signal $c(t)$ is the superposition of a number of sine functions with frequencies $\{ f_n \}$ and optimizable amplitudes $\{ c_n \}$ and phases $\{ \phi_n \}$, and the envelope $E(t)$ is a Gaussian function whose standard deviation $\sigma$ and mean $\mu$ you want to optimize:

$$ c(t) = \sum_n c_n \sin (2\pi f_n t + \phi_n), $$$$ E(t) = \exp \left\{ - \frac{(t- \mu)^2}{2\sigma^2} \right\}. $$You can separately create the two components of your controls (the fast-varying carrier signal and the slowly-varying envelope) and then multiply their sampleable functions of time (STF) to obtain the resulting controls.

The following example shows how you can obtain the result of an optimization where you apply an operation between two signals, and show how you can plot both the original signals and the resulting final signal. Note that to fetch the results you must first discretize the STF signals.

```
# Define standard matrices.
sigma_x = np.array([[0, 1], [1, 0]], dtype=complex)
# Define optimization parameters.
signal_frequencies = np.arange(0, 5) * 5e6 # Hz
total_duration = 1e-6 # s
# Set up the optimization graph.
with qctrl.create_graph() as graph:
t = qctrl.operations.identity_stf()
# Define the fast-varying carrier component of the controls.
amplitudes = qctrl.operations.optimization_variable(
len(signal_frequencies), -1e7, 1e7
)
phases = qctrl.operations.optimization_variable(
len(signal_frequencies), 0, 2 * np.pi
)
carrier = sum(
[
amplitudes[i] * qctrl.operations.sin(2 * np.pi * f * t + phases[i])
for i, f in enumerate(signal_frequencies)
]
)
# Define the slowly-varying Gaussian envelope.
sigma = qctrl.operations.optimization_variable(
1, total_duration / 600, total_duration / 6
)[0]
mu = qctrl.operations.optimization_variable(
1, total_duration / 6, 5 * total_duration / 6
)[0]
envelope = qctrl.operations.exp(-((t - mu) ** 2) / (2 * sigma ** 2))
# Multiply the two STFs.
resulting_signal = carrier * envelope
# Create a Hamiltonian operator from the controls.
hamiltonian = resulting_signal * sigma_x / 2
# Create infidelity as the cost function.
qctrl.operations.infidelity_stf(
hamiltonian=hamiltonian,
target_operator=qctrl.operations.target(operator=sigma_x),
sample_times=np.linspace(0, total_duration, 500),
name="infidelity",
)
# Discretize STFs to fetch them.
qctrl.operations.discretize_stf(
stf=carrier, duration=total_duration, segment_count=256, name="carrier"
)
qctrl.operations.discretize_stf(
stf=envelope, duration=total_duration, segment_count=256, name="envelope"
)
qctrl.operations.discretize_stf(
stf=resulting_signal,
duration=total_duration,
segment_count=256,
name="resulting_signal",
)
# Calculate optimization.
result = qctrl.functions.calculate_optimization(
optimization_count=3,
graph=graph,
cost_node_name="infidelity",
output_node_names=["carrier", "envelope", "resulting_signal"],
)
# Plot controls.
plot_controls(plt.figure(), result.output)
print(f"Infidelity: {result.cost}")
```

### Example: Mølmer–Sørensen system for an optimization

Next we present a full example of an ion chain system under a Mølmer–Sørensen interaction set up for optimization. The Hamiltonian of such a system is:

$$ H(t) = i\frac{\hbar}{2} \sum_{j=1}^N \sigma_{x, j} \sum_{p=1}^{3N} \eta_{pj} \left( \gamma_{j}(t) a_{p}^\dagger - \gamma^\ast_{j}(t) a_p \right), $$where $N$ is the number of ions, $\sigma_x$ is the Pauli $X$ operator, $a_p$ is the annihilation operator for the mode $p$, $\eta_{pj}$ is the Lamb–Dicke parameter of mode $p$ for the ion $j$, and $\gamma_j(t)$ is the drive applied ion the ion $j$.

Unlike the examples above, where we first modeled the Hamiltonian and then calculated the dynamics of the system based on the Hamiltonian, for an ion chain system we can use built-in nodes to calculate the dynamics of the system directly. Using information about the properties of the ions and the drive fields, the `ms_phases`

and `ms_displacements`

nodes compute, respectively, the relative phases for all ion pairs and the phase space displacements of the collective modes of motion. The `ms_infidelity`

node can then be used to calculate the infidelity of the realized system dynamics relative to a target gate (characterized by zero residual displacement and relative phases matching specified target values).

In this example we consider a chain of two ${}^{171}{\rm Yb}^{+}$ ions, with the goal of finding a pair of drives that produce zero residual displacement and a relative phase of $\pi/4$. You can visit the Optimizing Mølmer–Sørensen gates user guide for more details and examples, including information about how to calculate the ion chain properties that we hardcode in this example.

```
# Define ion chain properties
relative_detunings = np.array(
[
[-33076.63544983, -4700.0],
[-135006.15433009, -104700.0],
[-1304700.0, -1085084.75772933],
],
)
lamb_dicke_parameters = np.array(
[
[[-0.05519423, 0.05519423], [0.0547026, 0.0547026]],
[[-0.05707613, 0.05707613], [0.0564966, 0.0564966]],
[[0.0, 0.0], [0.0, 0.0]],
],
)
# Define physical constants
gamma_max = 2 * np.pi * 100e3 # Hz
segment_count = 100
duration = 2e-4 # s
# Define target phases as a lower triangular matrix
target_phases = np.array([[0, 0], [np.pi / 4, 0]])
# Create the graph describing the system
with qctrl.create_graph() as graph:
# Create the drive signals based on optimizable variables
drives = [
qctrl.operations.complex_pwc_signal(
moduli=qctrl.operations.optimization_variable(
count=segment_count, lower_bound=0, upper_bound=gamma_max
),
phases=qctrl.operations.optimization_variable(
count=segment_count,
lower_bound=0,
upper_bound=2 * np.pi,
is_upper_unbounded=True,
is_lower_unbounded=True,
),
duration=duration,
)
for _ in range(2)
]
# Calculate the system dynamics (relative phases and displacements),
# derived from the drive signals
relative_phases = qctrl.operations.ms_phases(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings,
)
displacements = qctrl.operations.ms_displacements(
drives=drives,
lamb_dicke_parameters=lamb_dicke_parameters,
relative_detunings=relative_detunings,
)
# Calculate the infidelity
infidelity = qctrl.operations.ms_infidelity(
phases=relative_phases,
displacements=displacements,
target_phases=target_phases,
name="infidelity",
)
# Call the function
result = qctrl.functions.calculate_optimization(
cost_node_name="infidelity",
output_node_names=[drives[0].name, drives[1].name],
graph=graph,
optimization_count=1,
)
# Visualize the results
plot_controls(
plt.figure(),
{
"$\gamma_1(t)$": result.output[drives[0].name],
"$\gamma_2(t)$": result.output[drives[1].name],
},
)
print("Infidelity:")
print(result.cost)
```

## Drives, shifts, and drifts

Some functions use a simpler, but less flexible, representation than graphs for describing quantum systems. Specifically, we consider Hamiltonians that can be decomposed into a sum of *drive*, *shift*, and *drift* terms. A drive term is a constant non-Hermitian operator multiplied by a complex scalar-valued function of time, plus its Hermitian conjugate. A shift term is a constant Hermitian operator multiplied by a real scalar-valued function of time. A drift term is a constant Hermitian operator. We can write such a Hamiltonian in the form

where $\gamma_j(t)$ are the complex scalar-valued functions of time modulating the non-Hermitian drive operators $C_j$, $\alpha_k(t)$ are the real scalar-valued functions of time modulating the Hermitian shift operators $A_k$, and $D_l$ are the static Hermitian drift operators.

In order to describe noise processes, each drive, shift, and drift term can also have an associated noise amplitude. Noise on drive and shift terms is used to represent *amplitude noise*, which is relative noise on the modulations $\gamma_j(t)$ or $\alpha_k(t)$ (corresponding to laser intensity fluctuations, for example). Noise on drift terms is used to represent *additive noise*, which is noise that directly modulates the drift operators $D_l$ (corresponding to dephasing noise, for example). Formally, when noises are included the Hamiltonian can be written

where $\beta_x(t)$ are noise amplitudes. Different functions support different forms for the noise amplitudes.

### Worked example: Two-qubit system for a filter function calculation

Consider a two-qubit system represented by the following Hamiltonian:

\begin{align*} H(t) = & \frac{\nu_a}{2} \sigma_{z,a} + \frac{\nu_b}{2} \sigma_{z,b} + \frac{\Omega_a(t)}{2} \sigma_{-,a} + \frac{\Omega^*_a(t)}{2} \sigma_{+,a} + \frac{\Omega_b(t)}{2} \sigma_{-,b} + \frac{\Omega^*_b(t)}{2} \sigma_{+,b}\\ & + \frac{\Lambda(t)}{2} \sigma_{-,a} \sigma_{+,b} + \frac{\Lambda^*(t)}{2} \sigma_{+,a} \sigma_{-,b} + \frac{\Delta_a(t)}{2} \sigma_{z,a} + \frac{\Delta_b(t)}{2} \sigma_{z,b} + \eta_a(t) \sigma_{z,a} + \eta_b(t)\sigma_{z,b}, \end{align*}where $\nu_k$ are the qubit detunings, $\Omega_k(t)$ are time-dependent Rabi rates (for example created by a microwave drive), $\Lambda(t)$ is a time-dependent coupling between the qubits, $\Delta_k(t)$ are time-dependent clock shifts, $\eta_k(t)$ are small, slowly-varying, stochastic dephasing noise processes and $\sigma_{\{+,-,z\}, k}$ are the Pauli matrices.

We consider custom pulses defined by the piecewise-constant functions:

\begin{align*} \Omega_a(t) &= \begin{cases} 2\pi\times 1.5e^{0.5\pi i}\times 10^{6} & 0\leq t\leq 2\times 10^{-6}\\ 0 & 2\times 10^{-6}\leq t\leq 4\times 10^{-6}\\ 2\pi\times 0.5e^{-0.5\pi i}\times 10^{6} & 4\times10^{-6}\leq t\leq 6\times 10^{-6} \end{cases}\\ \Omega_b(t) &= \begin{cases} 2\pi\times 1.0e^{0.25\pi i}\times 10^{6} & 0\leq t\leq 2\times 10^{-6}\\ 0 & 2\times 10^{-6}\leq t\leq 4\times 10^{-6}\\ 2\pi\times 0.5e^{1.25\pi i}\times 10^{6} & 4\times10^{-6}\leq t\leq 6\times 10^{-6} \end{cases}\\ \Lambda(t) &= \begin{cases} 0 & 0\leq t\leq 2\times 10^{-6}\\ 2\pi\times 1.0e^{\pi i}\times 10^{6} & 2\times 10^{-6}\leq t\leq 4\times 10^{-6}\\ 0 & 4\times10^{-6}\leq t\leq 6\times 10^{-6} \end{cases}\\ \Delta_a(t) &= \begin{cases} 2\pi\times 10.0\times 10^{6} & 0\leq t\leq 2\times 10^{-6}\\ 0 & 2\times 10^{-6}\leq t\leq 4\times 10^{-6}\\ 2\pi\times -10.0\times 10^{6} & 4\times10^{-6}\leq t\leq 6\times 10^{-6} \end{cases}\\ \Delta_b(t) &= \begin{cases} 2\pi\times -1.0\times 10^{6} & 0\leq t\leq 2\times 10^{-6}\\ 0 & 2\times 10^{-6}\leq t\leq 4\times 10^{-6}\\ 2\pi\times 1.5\times 10^{6} & 4\times10^{-6}\leq t\leq 6\times 10^{-6} \end{cases} \end{align*}Such a setup would be useful for performing simulation or characterizing susceptibility to noise via filter functions. In this case we create Python objects suitable for calculating a filter function, but the same procedure applies to other computations too.

#### Decomposing the quantum Hamiltonian into drives, shifts, and drifts

For functions that use the drives, shifts, and drifts representation of quantum systems, the first step in analyzing your system is to split the Hamiltonian of your quantum system into the appropriate set of drives, shifts, and drifts. In this case, the Hamiltonian is represented by eight terms:

- a
`drift`

term with`operator`

$\nu_a \sigma_{z,a}/2 + \nu_b \sigma_{z,b}/2$, - a
`drive`

term with`operator`

$\sigma_{-,a}/2$ and complex`control`

$\Omega_a(t)$, - a
`drive`

term with`operator`

$\sigma_{-,b}/2$ and complex`control`

$\Omega_b(t)$, - a
`drive`

term with`operator`

$\sigma_{-,a} \sigma_{+,b}/2$ and complex`control`

$\Lambda(t)$, - a
`shift`

term with`operator`

$\sigma_{z,a}/2$ and real`control`

$\Delta_a(t)$, - a
`shift`

term with`operator`

$\sigma_{z,b}/2$ and real`control`

$\Delta_b(t)$, - a
`drift`

term with`operator`

$\sigma_{z,a}$ and`noise`

amplitude $\eta_a(t)$, - a
`drift`

term with`operator`

$\sigma_{z,b}$ and`noise`

amplitude $\eta_b(t)$.

With this decomposition in hand, we may proceed to set up the appropriate Python objects.

#### Creating the control pulses

Each drive and shift control term in your system must have an associated control pulse, which describes the piecewise-constant scalar-valued function of time that modulates the control operator. A piecewise-constant function of time is represented as a list of segment objects, each describing the duration and value of a constant segment of the function. The function is obtained by concatenating these segments in time. The `qctrl.types.RealSegmentInput`

and `qctrl.types.ComplexSegmentInput`

objects represent such segments for real and complex piecewise-constant functions of time, respectively.

```
# Define physical constants
segment_duration = 1 * 1e-6 # s
# Define control pulses
drive_A_control = [
qctrl.types.ComplexSegmentInput(
duration=segment_duration,
value=2 * np.pi * 1.5 * np.exp(0.5j * np.pi) * 1e6,
),
qctrl.types.ComplexSegmentInput(duration=segment_duration, value=0),
qctrl.types.ComplexSegmentInput(
duration=segment_duration,
value=2 * np.pi * 0.5 * np.exp(-0.5 * np.pi) * 1e6,
),
]
drive_B_control = [
qctrl.types.ComplexSegmentInput(
duration=segment_duration,
value=2 * np.pi * 1.0 * np.exp(0.25j * np.pi) * 1e6,
),
qctrl.types.ComplexSegmentInput(duration=segment_duration, value=0),
qctrl.types.ComplexSegmentInput(
duration=segment_duration,
value=2 * np.pi * 1.0 * np.exp(1.25j * np.pi) * 1e6,
),
]
coupling_control = [
qctrl.types.ComplexSegmentInput(duration=segment_duration, value=0),
qctrl.types.ComplexSegmentInput(
duration=segment_duration,
value=2 * np.pi * 1.0 * np.exp(1.0j * np.pi) * 1e6,
),
qctrl.types.ComplexSegmentInput(duration=segment_duration, value=0),
]
clock_A_control = [
qctrl.types.RealSegmentInput(
duration=segment_duration,
value=2 * np.pi * 10.0 * 1e6,
),
qctrl.types.RealSegmentInput(duration=segment_duration, value=0),
qctrl.types.RealSegmentInput(
duration=segment_duration,
value=2 * np.pi * -10.0 * 1e6,
),
]
clock_B_control = [
qctrl.types.RealSegmentInput(
duration=segment_duration,
value=2 * np.pi * -1.0 * 1e6,
),
qctrl.types.RealSegmentInput(duration=segment_duration, value=0),
qctrl.types.RealSegmentInput(
duration=segment_duration,
value=2 * np.pi * 1.5 * 1e6,
),
]
```

```
plot_controls(
plt.figure(),
{
"$\Omega_a(t)$": map(attr.asdict, drive_A_control),
"$\Delta_a(t)$": map(attr.asdict, clock_A_control),
},
)
```

#### Creating the control terms

Each control term in your quantum system is represented as a Python object specific to the type of computation you are performing. In the case of filter functions, the relevant objects are `qctrl.types.filter_function.Drive`

, `qctrl.types.filter_function.Shift`

, and `qctrl.types.filter_function.Drift`

. These objects wrap an operator and, in the case of drives and shifts, a control.

```
# Define standard matrices
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_m = np.array([[0.0, 1.0], [0.0, 0.0]], dtype=complex)
# Define physical constants
nu_a = 2 * np.pi * 0.5 * 1e6 # Hz
nu_b = 2 * np.pi * 0.5 * 1e6 # Hz
# Define control objects
drift = qctrl.types.filter_function.Drift(
operator=nu_a * np.kron(sigma_z, identity) / 2
+ nu_b * np.kron(sigma_z, identity) / 2,
)
drive_A = qctrl.types.filter_function.Drive(
operator=np.kron(sigma_m, identity) / 2,
control=drive_A_control,
)
drive_B = qctrl.types.filter_function.Drive(
operator=np.kron(identity, sigma_m) / 2,
control=drive_B_control,
)
coupling = qctrl.types.filter_function.Drive(
operator=np.kron(sigma_m, sigma_m.T) / 2,
control=coupling_control,
)
clock_A = qctrl.types.filter_function.Shift(
operator=np.kron(sigma_z, identity) / 2,
control=clock_A_control,
)
clock_B = qctrl.types.filter_function.Shift(
operator=np.kron(identity, sigma_z) / 2,
control=clock_B_control,
)
```

#### Creating the noises

Each noise in your quantum system is represented as a drive, shift, or drift object with the `noise`

field set. In this example we have only drift noises, which are created as new terms in addition to the control terms defined above. In other cases you might need to consider control noises, which are represented by setting the `noise`

field of the control terms.

The data in the `noise`

field depends on the specific computation you are performing. In the case of filter functions, the `noise`

is simply a boolean indicating that the term is perturbed by noise (see the reference documentation for the exact meaning of this boolean).

```
dephasing_A = qctrl.types.filter_function.Drift(
operator=np.kron(sigma_z, identity),
noise=True,
)
dephasing_B = qctrl.types.filter_function.Drift(
operator=np.kron(identity, sigma_z),
noise=True,
)
```

#### Calling a function with the drives, shifts, and drifts

With the Python objects representing the Hamiltonian set up, you can perform computations using the functions in the Q-CTRL Python package. In this case we call the `qctrl.functions.calculate_filter_function`

function.

```
filter_function_result = qctrl.functions.calculate_filter_function(
duration=3 * 1e-6,
frequencies=np.linspace(0, 1e6, 100),
drives=[drive_A, drive_B, coupling],
shifts=[clock_A, clock_B],
drifts=[drift, dephasing_A],
)
```

```
print(filter_function_result.samples[0])
```

#### Summary

We have shown the general procedure for describing a quantum system as a set of drive, shift, and drift objects, and for performing computations on such systems. You can visit the other User guides and Application notes to see more examples of using this decomposition to perform computations, and the reference documentation to see more details about all functions and types in the Q-CTRL Python package.

### Example: One-qubit system for a simulation

Next we present a full example of a one-qubit system set up for simulation. We consider the Hamiltonian:

\begin{align*} H(t) = & \frac{\nu}{2} \sigma_z + \frac{\Omega(t)}{2} \sigma_- + \frac{\Omega^*(t)}{2} \sigma_+ + \frac{\Delta(t)}{2} \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, and $\sigma_k$ are the Pauli matrices.

In this case we use primitive control pulses $\Omega(t) = 2\pi\times 1.0e^{0.5\pi i} \times 10^6$ and $\Delta(t) = 2\pi\times 0.2\times 10^6$ (for $0\leq t\leq 2\times 10^{-6}$).

Below we show how to set up the drive, shift, and drift terms for a simulation of this system via the `qctrl.functions.coherent_simulation`

function. In this example we also show how to create dictionaries instead of class instances for the control pulses. Any inputs to callables in the `qctrl.functions`

or `qctrl.types`

namespaces can be provided as either classes or dictionaries; while using classes is safer and can benefit from autocompletion functionality, using dictionaries can lead to easier syntax in simple situations.

```
# Define physical constants
nu = 2 * np.pi * 0.5 * 1e6 # Hz
duration = 2 * 10e-6 # s
# Define control pulses
drive_control = [
{"duration": duration, "value": 2 * np.pi * 1.0 * np.exp(0.5j * np.pi) * 1e6},
]
shift_control = [
{"duration": duration, "value": 2 * np.pi * 0.2 * 1e6},
]
# Define control objects
drift = qctrl.types.coherent_simulation.Drift(operator=nu * sigma_z / 2)
drive = qctrl.types.coherent_simulation.Drive(
operator=sigma_m / 2,
control=drive_control,
)
clock = qctrl.types.coherent_simulation.Shift(
operator=sigma_z / 2,
control=shift_control,
)
# Call the simulation function
simulation_result = qctrl.functions.calculate_coherent_simulation(
duration=duration,
drives=[drive],
shifts=[clock],
drifts=[drift],
)
# Extract data from the result
print(simulation_result.samples[-1].evolution_operator)
```