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

## Imports and initialization

All usage of the Q-CTRL Python package begins by importing the qctrl package and starting a session.

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

from qctrl import Qctrl

# Starting a session with the API
qctrl = Qctrl()


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


#### Working with graphs

Graphs provide a powerful and efficient framework for performing computations remotely, and have some important differences when compared to libraries like NumPy. You might expect that if we printed the result object, we would see a value of 4. However, this is not the case:

print(result)

TensorNodeData(operation=<qctrlcommons.node.wrapper.Operation object at 0x12f8aa5d0>, shape=())


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)

TensorNodeData(operation=<qctrlcommons.node.wrapper.Operation object at 0x12f8ca3d0>, shape=(10,))


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
print(first_variable)

scaled_variables = variables * 5
print(scaled_variables)

# This is equivalent to the following code:
# added_variables = scaled_variables + np.linspace(0, 1, 10)

multiplied_variables = scaled_variables * variables
print(multiplied_variables)

TensorNodeData(operation=<qctrlcommons.node.wrapper.Operation object at 0x12fa180d0>, shape=())
TensorNodeData(operation=<qctrlcommons.node.wrapper.Operation object at 0x12fa18a50>, shape=(10,))
TensorNodeData(operation=<qctrlcommons.node.wrapper.Operation object at 0x12fa1d150>, shape=(10,))
TensorNodeData(operation=<qctrlcommons.node.wrapper.Operation object at 0x12fa1d310>, shape=(10,))


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. This is true for most operations that accept arrays, Tensors, PWCs, or STFs as input, although the exact meaning of shape-compatibility can depend on the operation. For basic arithmetic operations, the shapes of the objects must be equal or capable of being broadcasted.

For these basic arithmetic operations, 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 one their dimensions has different sizes 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 2$ with an operator of shape $1\times 1$.

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_integral=qctrl.operations.sinc_integral_function(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,
segments_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"],
)

Your task calculate_graph has completed in 3s.



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


#### 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 = delta * sigma_z / 2

# Sum the operators
hamiltonian = detuning + drive + clock

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

Your task calculate_optimization has completed in 3s.
Infidelity:
4.140243703432134e-12 ### 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
)
mu = qctrl.operations.optimization_variable(
1, total_duration / 6, 5 * total_duration / 6
)
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,
segments_count=256,
name="carrier",
)
qctrl.operations.discretize_stf(
stf=envelope,
duration=total_duration,
segments_count=256,
name="envelope",
)
qctrl.operations.discretize_stf(
stf=resulting_signal,
duration=total_duration,
segments_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}")

Your task calculate_optimization has completed in 3s.
Infidelity: 4.1611158962950867e-13 ### 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.name, drives.name],
graph=graph,
optimization_count=1,
)

# Visualize the results
plot_controls(
plt.figure(),
{
"$\gamma_1(t)$": result.output[drives.name],
"$\gamma_2(t)$": result.output[drives.name],
},
)

print("Infidelity:")
print(result.cost)

Your task calculate_optimization has completed in 3s.
Infidelity:
1.1653278342294016e-10 ## 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

\begin{align*} H(t) = \sum_{j=1}^{N_\mathrm{drives}} \left(\gamma_j(t)C_j + \mathrm{H.c}\right) + \sum_{k=1}^{N_\mathrm{shifts}} \alpha_k(t)A_k + \sum_{l=1}^{N_\mathrm{drifts}} D_l, \end{align*}

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

\begin{align*} H(t) =& \sum_{j=1}^{N_\mathrm{drives}} \left(1+\beta_{\gamma_j}(t)\right)\left(\gamma_j(t)C_j + \mathrm{H.c}\right) + \sum_{k=1}^{N_\mathrm{shifts}} \left(1+\beta_{\alpha_k}(t)\right)\alpha_k(t)A_k \\ &+ \sum_{l=1}^{N^\mathrm{(noiseless)}_\mathrm{drifts}} D_l + \sum_{l=1}^{N^\mathrm{(noisy)}_\mathrm{drifts}} \beta_{D_l}(t) D_l, \end{align*}

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,
),
]


#### Visualizing the control terms

The piecewise-constant control terms can be visualized using the Q-CTRL Python Visualizer package, as shown below. Also note that the Python objects can be conveniently converted into dictionaries using the attr.asdict function.

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],
)

Your task calculate_filter_function has completed in 4s.



#### Extracting the results

The evaluated values can be extracted from the result object. The specific form of outputs is discussed in detail in the more specific user guides, so here we simply print the sample of the calculated filter function at the zero frequency.

print(filter_function_result.samples)

Sample(frequency=0.0, inverse_power=6.822551858178849e-12, inverse_power_uncertainty=0.0, frequency_domain_noise_operator=array([[ 2.27124650e-06+5.62506303e-23j, -3.09917070e-07-6.30887052e-07j,
6.75073290e-07+6.35608643e-07j,  1.92611530e-07-6.64553782e-07j],
[-3.09917070e-07+6.30887052e-07j,  2.28288397e-06-9.92665080e-24j,
7.13582458e-07-3.55570823e-07j, -4.60271699e-07-3.17672637e-07j],
[ 6.75073290e-07-6.35608643e-07j,  7.13582458e-07+3.55570823e-07j,
-2.17199063e-06+3.30734068e-24j, -2.50591803e-07-6.48277715e-07j],
[ 1.92611530e-07+6.64553782e-07j, -4.60271699e-07+3.17672637e-07j,
-2.50591803e-07+6.48277715e-07j, -2.38213984e-06-4.96313202e-23j]]))


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

Your task calculate_coherent_simulation has completed in 3s.
[[ 0.26959181-5.52229725e-01j  0.78889961+1.87892933e-18j]
[-0.78889961+1.12516347e-18j  0.26959181+5.52229725e-01j]]