How to simulate large open system dynamics

Calculate the dynamics of a high-dimensional quantum system described by a GKS–Lindblad master equation

The Q-CTRL Python package enables you to simulate the evolution of open quantum systems that are interacting with their environment.

For larger systems, with Hilbert space dimensions more than roughly $10$, the exact evolution calculated using the techniques described above can be prohibitively expensive to compute. Instead, you can get better performance by using an approximate method to calculate the evolution. You can express the dynamics of these open quantum systems in terms of a master equation, such as the one described by Gorini, Kossakowski, and Sudarshan (GKS) and Lindblad. The tools described in this notebook allow you to solve this GKS–Lindblad master equation to obtain the time evolution for large systems using approximate methods.

Summary workflow

1. Select Krylov subspaces

The Q-CTRL Python package makes it easy to switch between the exact and approximate methods, simply by passing the krylov_subspace_dimension parameter to the qctrl.operations.density_matrix_evolution_pwc function. This parameter controls the precision of the method; a higher value leads to a more accurate solution, at the expense of computational efficiency.

2. Estimate the subspace dimension

The Q-CTRL Python package also provides a function qctrl.operations.estimated_krylov_subspace_dimension_arnoldi to estimate a value for this parameter based on properties of your system and your desired error tolerance. Note, however, that this estimate can sometimes be unnecessarily high; for performance-critical applications, for example when performing optimizations, you can manually fine-tune the parameter via preliminary simulations in order to find a smaller value that still yields sufficiently small errors.

In order to use this feature you will need to define a "sample" Hamiltonian, representing a typical value of your Hamiltonian. For time-dependent Hamiltonians, you should choose a value derived from order-of-magnitude estimates of the various time-dependent parameters in the system. Getting the values exactly right is usually unimportant, since the error bound used to calculate the estimate is quite loose. You also need to set the duration (total duration of the evolution you want to simulate) and maximum segment duration (longest duration of any single segment in that evolution).

3. Select sparse matrices (optional)

When you pass a krylov_subspace_dimension, you also have the option of setting up your system with sparse matrices. Using sparse matrices can lead to performance improvements if your Hamiltonian and collapse operators contain mostly zeros.

Worked example: Superconducting resonator coupled to transmon qubits

In this example, we consider a superconducting qubit coupled to a transmission line resonator (truncated to $10$ levels), with a drive on the resonator. We show how to simulate a measurement process that uses the quadratures of the resonator output to probe the qubit state.

Assuming that the qubit is far detuned from the resonator (the dispersive approximation), and that the resonator drive is off resonance from the qubit transition frequency, in the absence of losses the system can be described by the Hamiltonian (see Gambetta et al. and Blais et al.):

$$ H(t) = \frac{\omega_q+g^2/\Delta}{2} \sigma_z + (\omega_r-\omega_d) a^\dagger a + \frac{g^2}{\Delta} a^\dagger a \sigma_z + \Gamma(t) a^\dagger + \Gamma^*(t) a, $$

where $a$ is the annihilation operator for the resonator, $\sigma_z$ is the Pauli Z matrix for the qubit, $\Gamma(t)$ is the resonator drive, $\omega_r$ is the resonator frequency, $\omega_q$ is the qubit transition frequency, $\omega_d$ is the resonator drive frequency, $\Delta=\omega_q-\omega_r$ is the qubit detuning from the resonator, and $g$ is the coupling strength. Taking into account the photon loss $\kappa$ for the resonator and qubit decay $\gamma$, the system evolution follows the GKS–Lindblad equation:

$$ \frac{\mathrm{d} \rho(t)}{\mathrm{d} t} = - i \left[ H(t), \rho(t) \right] + \kappa \left[ a \rho(t) a^\dagger - \frac{1}{2} \rho(t) a^\dagger a - \frac{1}{2} a^\dagger a \rho(t) \right] + \gamma \left[ \sigma_- \rho(t) \sigma_+ - \frac{1}{2} \rho(t) \sigma_+ \sigma_- - \frac{1}{2} \sigma_+\sigma_- \rho(t) \right], $$

where $\sigma_{\pm} = (\sigma_x \mp i\sigma_y)/2$.

Below we show how you can simulate this system using BOULDER OPAL. Note that similar code can be used to perform optimizations, for example if you want to design custom measurement pulses optimized for your particular system.

import matplotlib.pyplot as plt
import numpy as np
from scipy.sparse import coo_matrix
from qctrlopencontrols import new_corpse_control, new_primitive_control
from qctrlvisualizer import get_qctrl_style, plot_controls

from qctrl import Qctrl

# Start a session with the API
qctrl = Qctrl()
# Define standard matrices.
identity_q = np.identity(2, 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)

n_r = 10  # 10 levels in the resonator
identity_r = np.identity(n_r, dtype=complex)
a = np.diag(np.sqrt(np.arange(1, n_r)), k=1)

# Define system parameters.
omega_r = 2 * np.pi * 10 * 1e9  # Hz
omega_q = 2 * np.pi * 8 * 1e9  # Hz
omega_d = 2 * np.pi * 10 * 1e9  # Hz
delta = omega_q - omega_r
g = 2 * np.pi * 0.05 * 1e9  # Hz
kappa = 2 * np.pi * 1e6  # Hz
gamma = 2 * np.pi * 1e3  # Hz
duration = 100 * 1e-9  # s

# Define drive parameters.
drive_amplitude = 2 * np.pi * 0.02 * 1e9  # Hz
drive_rise_duration = 20 * 1e-9  # s
segment_count = 1000

# Set up Hamiltonian. We use the `coo_matrix` function from SciPy to
# convert the matrix to a sparse matrix, which can lead to faster
# computations for large systems.
base_hamiltonian_value = coo_matrix(
    (omega_q + g ** 2 / delta) / 2 * np.kron(identity_r, sigma_z)
    + (omega_r - omega_d) * np.kron(a.T.conj() @ a, identity_q)
    + g ** 2 / delta * np.kron(a.T.conj() @ a, sigma_z)
drive_hamiltonian_value = coo_matrix(
    drive_amplitude * np.kron(a.T.conj() + a, identity_q)

# Set up Lindblad terms.
lindblad_terms = [
    (kappa, coo_matrix(np.kron(a, identity_q))),
    (gamma, coo_matrix(np.kron(identity_r, sigma_m))),

# Set up density matrices for the qubit in |0> and |1> states.
psi_0 = np.kron([1] + [0] * (n_r - 1), [1, 0])
rho_0 = np.outer(psi_0, psi_0)
psi_1 = np.kron([1] + [0] * (n_r - 1), [0, 1])
rho_1 = np.outer(psi_1, psi_1)

# Set up observables for the quadratures of the resonator field.
I = np.kron(a.T.conj() + a, identity_q)
Q = np.kron(-1j * (a - a.T.conj()), identity_q)

# Calculate an appropriate Krylov subspace dimension for the system and
# desired error tolerance.
# You need to define a "sample" Hamiltonian, representing a typical
# value of your Hamiltonian. In this case we use
# the value when the system is maximally-driven.
# You also need to set the duration (total duration of the evolution
# you want to simulate) and maximum segment duration (longest duration
# of any single segment in that evolution).
with qctrl.create_graph() as graph:
        hamiltonian_sample=coo_matrix(base_hamiltonian_value + drive_hamiltonian_value),
        maximum_segment_duration=duration / segment_count,
krylov_subspace_dimension = qctrl.functions.calculate_graph(

print(f"Estimated Krylov subspace dimension: {krylov_subspace_dimension}")

# Run the simulation.
with qctrl.create_graph() as graph:
    base_hamiltonian = qctrl.operations.constant_sparse_pwc_operator(

    # Prepare a tanh envelope for the drive signal.
    drive_envelope_times = np.linspace(0, duration, segment_count)
    drive_envelope_values = np.tanh(drive_envelope_times / drive_rise_duration)
    drive_envelope = qctrl.operations.pwc_signal(

    drive_hamiltonian = qctrl.operations.sparse_pwc_operator(

    density_matrices = qctrl.operations.density_matrix_evolution_pwc(
        initial_density_matrix=np.array([rho_0, rho_1]),
            [base_hamiltonian, drive_hamiltonian]
        sample_times=np.linspace(0, duration, 500),

    I_values_0 = qctrl.operations.trace(density_matrices[0] @ I, name="I_0")
    Q_values_0 = qctrl.operations.trace(density_matrices[0] @ Q, name="Q_0")
    I_values_1 = qctrl.operations.trace(density_matrices[1] @ I, name="I_1")
    Q_values_1 = qctrl.operations.trace(density_matrices[1] @ Q, name="Q_1")

result = qctrl.functions.calculate_graph(
    output_node_names=["I_0", "Q_0", "I_1", "Q_1"],

# Plot the results. Note how the I quadrature can be used to distinguish
# between the qubit states.
plt.title("I quadrature")
    np.linspace(0, duration, 500) * 1e9,
    np.linspace(0, duration, 500) * 1e9,
plt.xlabel("Time (ns)")
plt.ylabel("I (arb. units)")

plt.title("Q quadrature")
    np.linspace(0, duration, 500) * 1e9,
    np.linspace(0, duration, 500) * 1e9,
plt.xlabel("Time (ns)")
plt.ylabel("Q (arb. units)")
Your task calculate_graph has completed.
Estimated Krylov subspace dimension: 22
Your task calculate_graph has started.
Your task calculate_graph has completed.