Q-CTRL logo

Jupyter Get the notebook

How to define continuous analytical Hamiltonians

Use analytical expressions to construct your Hamiltonian

The Boulder Opal optimization and simulation engines can be applied to a wide range of quantum systems thanks to the flexible dataflow graph framework.

In particular, you can define analytical Hamiltonians by taking advantage of Bouldar Opal’s basic mathematical operations. In this notebook we show how to define signals from their analytical expressions and use them for simulation or optimization. Note that the Boulder Opal toolkit contains a library of predefined analytical signals as piecewise-constant (PWC) or sampleable (STF) functions that you can directly use without having to define them from the ground up. For more information, you can also read the How to create analytical signals for simulation and optimization user guide.

Summary workflow

1. Create an identity Stf node and define analytical signals

Once a graph has been set up to describe the computation, use the graph.identity_stf operation to create an STF representing the identity function $f(t) = t$.

Apply a sequence of mathematical functions to the identity function to create STF nodes representing your desired signals.

2. Construct the Hamiltonian

Create the Hamiltonian terms by multiplying the signals by the appropriate operators. You can directly use STF nodes created this way with operations such as graph.time_evolution_operators_stf or graph.infidelity_stf, or you can discretize them into PWCs with graph.discretize_stf. For more information on STFs and PWCs, see the Working with time-dependent functions in Boulder Opal topic.

3. Execute the graph for the simulation/optimization

With the graph object created, a simulation can be run using the qctrl.functions.calculate_graph function. The outputs and the graph must be provided. The function returns the results of the simulation.

Alternatively, an optimization can be run using the qctrl.functions.calculate_optimization function or the qctrl.functions.calculate_stochastic_optimization function. The cost, the outputs, and the graph must be provided. The function returns the results of the optimization.

Worked example: Simulation with analytical signals

In this example, we will simulate a Landau–Zener transition in a two-level system.

The Hamiltonian for the system is \begin{equation} H(t) = \Omega(t) \sigma_x + \Delta(t) \sigma_z , \end{equation} where $\Omega(t)$ is the Rabi coupling and $\Delta(t)$ is the detuning.

For the Rabi coupling we will implement a pulse of the form \begin{equation} \Omega(t) = \Omega_0 \sin^2 (\pi t / T) , \end{equation} and for the detuning a sweep of the form \begin{equation} \Delta(t) = - \Delta_0 \cos(\pi t / T) , \end{equation} where $\Omega_0$ and $\Delta_0$ are the maximum values of the Rabi coupling and detuning, and $T$ is the signal duration.

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

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

# Create the data flow graph describing the system.
graph = qctrl.create_graph()

# Define the duration of the signals.
duration = 2e-6  # s

# Define the number of samples for the discretizations.
sample_count = 512
sample_times = np.linspace(0.0, duration, sample_count)

# Define maximum values for the signals.
omega_max = 1.5e6  # rad/s
delta_max = 2.1e6  # rad/s

# Define the signals.
t = graph.identity_stf()
omega = omega_max * graph.sin(np.pi * t / duration) ** 2
delta = -delta_max * graph.cos(np.pi * t / duration)

# Construct the Hamiltonian.
hamiltonian = omega * graph.pauli_matrix("X") + delta * graph.pauli_matrix("Z")

# Compute the time-evolution unitaries.
unitaries = graph.time_evolution_operators_stf(
    hamiltonian=hamiltonian, sample_times=sample_times

# Define the initial state.
initial_state = graph.fock_state(2, 0)[:, None]

# Evolve the initial state with the unitary to get the final state.
evolved_states = unitaries @ initial_state
evolved_states.name = "states"

# Discretize the signals as STFs can't be directly exported.
graph.discretize_stf(omega, duration, sample_count, name=r"$\Omega$")
graph.discretize_stf(delta, duration, sample_count, name=r"$\Delta$")

# Run simulation.
result = qctrl.functions.calculate_graph(
    graph=graph, output_node_names=["states", r"$\Omega$", r"$\Delta$"]

# Plot the analytical signals used.
    {name: result.output[name] for name in [r"$\Omega$", r"$\Delta$"]}, smooth=True
plt.suptitle("Landau–Zener signals")

# Plot the population dynamics.
populations = np.abs(result.output["states"]["value"].squeeze()) ** 2
plot_populations(sample_times, {rf"$|{k}\rangle$": populations[:, k] for k in range(2)})
plt.suptitle("Landau–Zener transition dynamics")
Your task calculate_graph (action_id="1404015") has completed.



Worked example: Optimization with analytical signals

In this example, we will implement an X gate in a noisy single-qubit system, by optimizing the parameters of two signals.

The Hamiltonian for the system is \begin{equation} H(t) = \alpha(t) \sigma_{z} + \frac{1}{2}\left(\gamma(t)\sigma_{-} + \gamma^*(t)\sigma_{+}\right) + \delta \sigma_{z} + \zeta(t) \sigma_{z} , \end{equation} where $\alpha(t)$ and $\gamma(t)$ are signals with optimizable parameters, $\delta$ is the qubit detuning, and $\zeta(t)$ is a dephasing noise process.

The $\alpha(t)$ signals acting on $\sigma_{z}$ is an arctangent ramp, defined by \begin{equation} \alpha(t) = \alpha_0 \arctan\left( \frac{t - t_0}{t_\mathrm{ramp}} \right) , \end{equation} where $\pm a_0$ are the asymptotic values of the signal, $t_0$ is the center time, and $t_\mathrm{ramp}$ is the characteristic duration of the ramp.

The $\gamma(t)$ signal acting on $\sigma_{-}$ is a Gaussian signal, defined by \begin{equation} \gamma(t) = \gamma_0 \exp \left(- \frac{(t-t_0)^2}{2\sigma^2} \right) , \end{equation} where $\gamma_0$ is the amplitude, $\sigma$ is the width of the Gaussian, and $t_0$ is the center of the signal. As the Boulder Opal signal library includes a node to create a Gaussian pulse, graph.signals.gaussian_pulse_stf, we will create $\gamma(t)$ with it rather than defining it from scratch.

The Boulder Opal Toolkits are currently in beta phase of development. Breaking changes may be introduced.

We will fix the characteristic ramp time $t_\mathrm{ramp}$ and Gaussian width $\sigma$, and aim to find optimal values for the rest of signal parameters to implement the target quantum gate.

In this case, we will discretize the smooth analytical signals into piecewise-constant functions to simulate a discrete interface with an experiment.

graph = qctrl.create_graph()

# Define the duration of the signals.
duration = 3e-6  # s
sample_count = 64

# Define α(t) as an arctan ramp.
t = graph.identity_stf()
alpha_ramp = duration / 8
alpha_zero = graph.optimization_variable(
    count=1, lower_bound=-3e7, upper_bound=3e7, name="alpha_amplitude"
alpha_center = graph.optimization_variable(
    count=1, lower_bound=0.0, upper_bound=duration, name="alpha_center"
alpha_stf = alpha_zero[0] * graph.arctan((t - alpha_center[0]) / alpha_ramp)

# Define γ(t) as a Gaussian pulse.
gamma_max = 2.0 * np.pi * 1e6
gamma_width = duration / 8
gamma_modulus = graph.optimization_variable(
    count=1, lower_bound=-gamma_max, upper_bound=gamma_max
gamma_phase = graph.optimization_variable(
    upper_bound=2.0 * np.pi,
gamma_amplitude = gamma_modulus * graph.exp(1j * gamma_phase)
gamma_amplitude.name = "gamma_amplitude"

gamma_center_time = graph.optimization_variable(
    count=1, lower_bound=0.0, upper_bound=duration, name="gamma_center"
gamma_stf = graph.signals.gaussian_pulse_stf(
    gamma_amplitude, gamma_width, gamma_center_time

# Discretize the signals.
alpha = graph.discretize_stf(alpha_stf, duration, sample_count, name=r"$\alpha$")
gamma = graph.discretize_stf(gamma_stf, duration, sample_count, name=r"$\gamma$")

# Create Hamiltonian.
hamiltonian = alpha * graph.pauli_matrix("Z") + graph.hermitian_part(
    gamma * graph.pauli_matrix("M")

# Add detuning term.
delta = 2 * np.pi * 0.25e6  # rad/s
hamiltonian += delta * graph.pauli_matrix("Z")

# Define target gate.
target = graph.target(operator=graph.pauli_matrix("X"))

# Add dephasing noise term.
zeta = 2 * np.pi * 20e3  # rad/s
dephasing = zeta * graph.pauli_matrix("Z")

# Robust infidelity.
robust_infidelity = graph.infidelity_pwc(

result = qctrl.functions.calculate_optimization(

print(f"Optimized cost: {result.cost:.3e}")

# Print out the converged signal parameters.
print("Optimal signal parameters")
print(f"\tα(t) amplitude:\t{result.output['alpha_amplitude']['value'][0]:.3e} rad/s")
print(f"\tα(t) center:\t{result.output['alpha_center']['value'][0]:.3e} s")
print(f"\tγ(t) amplitude:\t{result.output['gamma_amplitude']['value'][0]:.3e} rad/s")
print(f"\tγ(t) center:\t{result.output['gamma_center']['value'][0]:.3e} s")

# Plot the optimal analytical signals found.
plot_controls({name: result.output[name] for name in [r"$\alpha$", r"$\gamma$"]})
plt.suptitle("Optimal signals")
Your task calculate_optimization (action_id="1404017") has completed.
Optimized cost: 6.748e-04
Optimal signal parameters
	α(t) amplitude:	3.943e+06 rad/s
	α(t) center:	1.683e-06 s
	γ(t) amplitude:	1.729e+06-6.041e+06j rad/s
	γ(t) center:	1.569e-06 s


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

Package Version
Python 3.9.12
matplotlib 3.5.1
numpy 1.21.6
scipy 1.7.3
qctrl 19.7.2
qctrl-commons 17.6.0
boulder-opal-toolkits 2.0.0-beta
qctrl-visualizer 4.3.0