seQuencing logo


sequencing is designed to make it easy to create and simulate time-dependent Hamiltonians.

This notebook is intended to show, at a high level, some of what sequencing can do. For a more in-depth look at the inner workings of sequencing and more advanced applications, check out the tutorials.


The general problem sequencing is designed to set up and solve is the time evolution of the state \(|\psi(t)\rangle\) or density matrix \(\rho(t)\) of a system composed of \(n\) oscillators or “modes” — each with its own nonlinearity, dimension, coherence properties, and interactions with other modes — acted upon by realistic time-dependent controls.

The typical time-dependent Hamiltonian constructed using sequencing has the following form (taking \(\hbar=1\)):

\[\begin{split} \begin{align} \hat{H}(t) &= \sum_{i=1}^n \hat{H}_{\mathrm{mode, }i} + \sum_{i\neq j}\hat{H}_{\mathrm{int, }ij} + \sum_{k}\hat{H}_{\mathrm{control, }k}(t)\\ &= \sum_{i=1}^n \delta_i\hat{a}_i^\dagger\hat{a}_i + \frac{K}{2}(\hat{a}_i^\dagger)^2(\hat{a}_i)^2\\ &+ \sum_{i\neq j}\chi_{ij}\hat{a}_i^\dagger\hat{a}_i\hat{a}_j^\dagger\hat{a}_j\\ &+ \sum_{\{\hat{A}_k\}}c_{\hat{A}_k}(t)\hat{A}_k \end{align}\end{split}\]
  • \(\hat{a}_i\) is the mode annihilation (lowering) operator.

  • \(\delta_i\) is the detuning of the mode \(i\) relative to the chosen rotating frame.

  • \(K_i\) is the Kerr nonlinearity (or self-Kerr) of mode \(i\).

  • \(\chi_{ij}\) is the cross-Kerr or dispersive shift between mode \(i\) and mode \(j\).

  • \(\{\hat{A}_k\}\) are a set of Hermitian control operators with complex time-dependent coefficients \(c_{\hat{A}_k}(t)\), where each \(\hat{A}_k\) may be composed of operators acting on one or more modes.


  • Although only the first-order Kerr nonlinearity \((\hat{a}_i^\dagger)^2(\hat{a}_i)^2\) is included by default, higher-order nonlinearities can easily be added.

  • Although the default interaction term is that of a dispersive coupling, \(\hat{H}_{\mathrm{int, }ij}=\chi_{ij}\hat{a}_i^\dagger\hat{a}_i\hat{a}_j^\dagger\hat{a}_j\), any type of multi-mode interaction can be included.

The finite coherence of each mode can also be included in the form of Lindblad collapse operators (see the Lindblad master equation section of the QuTiP documentation). The default collapse operators implemented in sequencing are:

  • \(\hat{C}_{\uparrow,i} = \sqrt{\gamma_{\uparrow,i}}\,\hat{a}_i^\dagger\), where \(\gamma_{\uparrow,i} = p_{\mathrm{therm, }i}T_{1,i}^{-1}\) is the mode’s energy excitation rate, computed as thermal population divided by \(T_1\).

  • \(\hat{C}_{\downarrow,i} = \sqrt{\gamma_{\downarrow,i}}\,\hat{a}_i\), where \(\gamma_{\downarrow,i} = T_{1,i}^{-1} - \gamma_{\uparrow,i}\) is the mode’s energy decay rate.

  • \(\hat{C}_{\phi,i} = \sqrt{2\gamma_{\phi,i}}\,\hat{a}_i^\dagger\hat{a}_i\), where \(\gamma_{\phi,i} = T_{\phi,i}^{-1} = T_{2,i}^{-1} - (2T_{1,i})^{-1}\) is the mode’s pure dephasing rate.

Only collapse operators with nonzero coefficients are included in simulations. By default, each mode is set to have ideal coherence properties (\(p_{\mathrm{therm, }i}=0\), \(T_{1,i}=T_{2,i}=\infty\)), meaning that no collapse operators are used and the evolution is unitary.

%config InlineBackend.figure_formats = ['svg']
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

from sequencing import Transmon, System, get_sequence, sync, delay

Create a quantum Mode and System

qubit = Transmon('qubit', levels=3)
qubit.anharmonicity = -200e-3 # GHz
# A System can contain many interacting or non-interacting Modes,
# but here we have just one.
system = System('system', modes=[qubit])

Inspect and adjust Pulse parameters

print(list(qubit.pulses), end='\n\n')
print(qubit.gaussian_pulse, end='\n\n')
['smoothed_constant_pulse', 'gaussian_pulse']

GaussianPulse(name='gaussian_pulse', cls='sequencing.pulses.GaussianPulse', amp=1.0, detune=0.0, phase=0.0, dt=1, noise_sigma=0.0, noise_alpha=0.0, scale_noise=False, sigma=10, chop=4, drag=0.0)

  "amp": 1.0,
  "chop": 4,
  "cls": "sequencing.pulses.GaussianPulse",
  "detune": 0.0,
  "drag": 0.0,
  "dt": 1,
  "name": "gaussian_pulse",
  "noise_alpha": 0.0,
  "noise_sigma": 0.0,
  "phase": 0.0,
  "scale_noise": false,
  "sigma": 10

Plot the Pulse waveform

ax = qubit.gaussian_pulse.plot()
# Adjust the Gaussian width, sigma
qubit.gaussian_pulse.sigma = 15 # ns
ax = qubit.gaussian_pulse.plot(ax=ax)

The values plotted above represent time-dependent coefficients for the x and y operators defined by the Mode class:

Quantum object: dims = [[3], [3]], shape = (3, 3), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 1.0 & 0.0\\1.0 & 0.0 & 1.414\\0.0 & 1.414 & 0.0\\\end{array}\right)\end{equation*}
Quantum object: dims = [[3], [3]], shape = (3, 3), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & -1.0j & 0.0\\1.0j & 0.0 & -1.414j\\0.0 & 1.414j & 0.0\\\end{array}\right)\end{equation*}

If we temporarily set the number of levels in our Transmon to 2, we will see that x and y take the familiar form of the Pauli operators, \(\sigma_x\) and \(\sigma_y\).

with qubit.temporarily_set(levels=2):
Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 1.0\\1.0 & 0.0\\\end{array}\right)\end{equation*}
Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & -1.0j\\1.0j & 0.0\\\end{array}\right)\end{equation*}

Simulate a Gaussian \(\pi\) pulse

Define the initial state of the system, and operators for which you would like to compute expectation values.

init_state = qubit.fock(0) # |g> state
g_op = qubit.fock_dm(0) # |g><g| operator
e_op = qubit.fock_dm(1) # |e><e| operator

Generate the pulse sequence

And plot the time-dependent Hamiltonian coefficients.

# Define the pulse sequence
seq = get_sequence(system)

fig, ax = seq.plot_coefficients(subplots=False)
ax.set_title('Hamiltonian coefficients for a Gaussian pulse')
ax.set_xlabel('Time [ns]')
ax.set_ylabel('Hamiltonian coefficient [GHz]');

Simulate the pulse sequence

And plot expectation values vs. time.

# Run the master equation simulation
result =, e_ops=[g_op, e_op])

fig, ax = plt.subplots()
ax.plot(result.times, result.expect[0], label='$|g\\rangle$')
ax.plot(result.times, result.expect[1], label='$|e\\rangle$')
ax.set_xlabel('Time [ns]')
ax.set_title('Qubit population vs. time - Gaussian $\pi$ pulse')

Simulate a \(T_1\) measurement

def fit_exp_decay(xs, ys):
    slope, offset = np.polyfit(xs, np.log(ys), 1)
    amp = np.exp(offset)
    tau = -1 / slope
    return amp, tau

The pulse sequence for a \(T_1\) measurement is:

for delay_time in delay_times:
    Rx(pi) # excite the qubit from |g> -> |e>
    measure P(|e>)

Because there are no operations on the qubit after the delay time, we can simulate the whole measurement in a single master equation simulation.

Execute the sequence

qubit.t1 = 10e3 # 10e3 ns == 10 us

tmax = 40e3 # maximum delay time, 40 us

seq = get_sequence(system)
# T1 sequence

result =, e_ops=[g_op, e_op])

ts = result.times / 1e3 # [us]
g_pops = result.expect[0] # P(|g>)
e_pops = result.expect[1] # P(|e>)

Analyze the results

Fit to the model for energy relaxation, \(P_e(t) = e^{-\frac{t}{T_1}}\), by fitting \(\log(P_e)\) vs. \(t\) to a line.

# Only fit P(|e>) starting after the pi pulse is finished
t0 = qubit.gaussian_pulse.sigma * qubit.gaussian_pulse.chop
fit_ts = ts[t0:]
fit_amp, fit_t1 = fit_exp_decay(fit_ts, e_pops[t0:])

Plot the results and the fit.

fig, (ax, bx) = plt.subplots(1, 2, sharey=True, figsize=(9,4))
# plot population vs. time during the pi pulse
ax.plot(ts[:t0], g_pops[:t0], lw=3, label='$P(|g\\rangle)$')
ax.plot(ts[:t0], e_pops[:t0], lw=3, label='$P(|e\\rangle)$')
ax.set_title('$\pi$ pulse')

# plot population vs. time for the whole sequence
bx.plot(ts, g_pops, lw=3, label='$P(|g\\rangle)$')
bx.plot(ts, e_pops, lw=3, label='$P(|e\\rangle)$')
bx.plot(fit_ts, np.exp(-fit_ts / fit_t1), 'k--', label=f'Fit, $T_1$ = {fit_t1:.3f} $\mu$s')
bx.set_title('Full evolution')
for a in (ax, bx):
    a.set_xlabel('Time [$\mu$s]')
fig.suptitle('$T_1$ sequence');

Simulate a \(T_2\) echo measurement

The pulse sequence for a \(T_2\) echo measurement is:

for delay_time in delay_times:
    Ry(pi / 2) # generate a superposition state |g> + e^(i*phi)|e>
    delay(delay_time / 2)
    Rx(pi) # perform the echo
    delay(delay_time / 2)
    Ry(pi / 2) # rotate back to one of the poles of the Bloch sphere
    measure P(|e>)

In this case, unlike for the \(T_1\) measurement, we need to perform a separate master equation simulation for each delay time.

Visualize the \(T_2\) echo sequence

dt = 200 # ns

seq = get_sequence(system)


fig, ax = seq.plot_coefficients(subplots=False)
ax.set_title(f'$T_2$ echo sequence at dt = {dt} ns');

Execute the sequence

qubit.t2 = 10e3 # 10 us
qubit.t1 = np.inf # Pure dephasing only

delay_times = 1e3 * np.linspace(0, 30, 11) # sweep delay from 0 to 30 us

g_pops = []
e_pops = []

print('Running dt [us] =', end=' ')
for dt in delay_times:
    print(dt/1e3, end=' ')

    seq = get_sequence(system)

    # T2 echo sequence

    result =, e_ops=[g_op, e_op])
Running dt [us] = 0.0 3.0 6.0 9.0 12.0 15.0 18.0 21.0 24.0 27.0 30.0

Analyze the results

ts = delay_times / 1e3 # [us]
# rescale |e> population to decay from 1 to 0 for fitting
ys = 2 * np.array(e_pops) - 1
fit_amp, fit_t2 = fit_exp_decay(ts, ys)
# invert the rescaling for plotting
fit_ys = (1 + np.exp(-ts / fit_t1)) / 2
fig, ax = plt.subplots()
ax.plot(ts, g_pops, 'o', label='$P(|g\\rangle)$')
ax.plot(ts, e_pops, 'o', label='$P(|e\\rangle)$')
ax.plot(ts, fit_ys, 'k--', label=f'Fit, $T_2$ = {fit_t2:.3f} $\mu$s')
ax.set_xlabel('Delay time [$\mu$s]')
ax.set_title('$T_2$ echo sequence')

Learn more

For more in depth tutorials and examples, see Tutorials.

from qutip.ipynbtools import version_table
Number of CPUs2
Python3.8.6 (default, Oct 19 2020, 15:10:29) [GCC 7.5.0]
OSposix [linux]
Tue Aug 30 19:19:17 2022 UTC
[ ]: