Simulating the dynamics of quantum systems in the presence of noise

The Q-CTRL Python package enables you to simulate the dynamics of quantum systems that are affected by various noise processes. Such simulations can provide useful insights into the expected real-world performance of candidate control solutions. In this notebook we show how to use the Q-CTRL Python package to perform simulations of systems affected by different types of noise.

Imports and initialization

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

# Essential imports
import numpy as np
from qctrl import Qctrl

# Predefined pulse imports
from qctrlopencontrols import new_predefined_driven_control

# Plotting imports
import matplotlib.pyplot as plt

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

Worked example: single qubit subject to dephasing noise

In this example we will show how to simulate the dynamics of a single qubit experiencing stochastic dephasing noise, driven by a $\pi/2$ CORPSE pulse. The Hamiltonian of the quantum system is:

\begin{align*} H(t) = & \frac{\Omega(t)}{2} \sigma_- + \frac{\Omega^*(t)}{2} \sigma_+ + \frac{\eta(t)}{2} \sigma_z, \end{align*}

where $\Omega(t)$ is a time-dependent Rabi rate, $\eta(t)$ is a stochastic dephasing noise process and $\sigma_k$ are the Pauli matrices.

We assume that $\eta(t)$ represents pink noise, as defined below.

def pink(frequencies, frequency_cutoff):
    return 1/(frequencies + frequency_cutoff)

frequencies = np.linspace(0, 2*1e6, 1000)

power_densities = 4e9*pink(frequencies=frequencies, frequency_cutoff=0.05*1e6)

plt.plot(frequencies/1e6, power_densities*1e6, c="#FB00A5")
plt.fill_between(frequencies/1e6, 0, power_densities*1e6, alpha=0.25, color="#FB00A5")
plt.xlabel("Frequency (MHz)")
plt.ylabel("Power density (1/MHz)")
plt.title("Dephasing noise spectral density")

We now proceed to show how the system, driven by stochastic noise drawn from this power spectral density, may be simulated using the Q-CTRL Python package.

Creating the system, control, pulse and noise

As described in the Setup feature guide, we first set up Python objects representing the system, control, pulse and noise.

# Define standard matrices
identity = np.array([[1., 0.],[0., 1.]], dtype=np.complex)
sigma_x = np.array([[0., 1.],[1., 0.]], dtype=np.complex)
sigma_z = np.array([[1., 0.],[0., -1.]], dtype=np.complex)
sigma_m = np.array([[0., 0.],[1., 0.]], dtype=np.complex)
square_root_sigma_x = 0.5*np.array([[1+1j, 1-1j], [1-1j, 1+1j]], dtype=np.complex)

# Define control parameters
omega_max = 2*np.pi * 1e6 #Hz
total_rotation = np.pi/2

# Define system object
system =
    name='Single qubit with dephasing',

# Define control object
drive = qctrl.factories.drive_controls.create(
    name='Rabi rate',

# Define pulse object using pulses from Q-CTRL Open Controls
pulse = new_predefined_driven_control(

drive_pulse = qctrl.factories.custom_pulses.create(
    segments=[{'duration': d, 'value': v} 
              for d, v in zip(pulse.durations,
                              pulse.rabi_rates * np.exp(1j*pulse.azimuthal_angles))],

# Define noise object
dephasing_noise = qctrl.factories.additive_noises.create(

Creating the noise spectral density

In order to perform a realistic stochastic simulation, we must provide the noise spectral density describing the behavior of the dephasing noise process, which is represented as a piecewise_linear_noise_spectral_density object. Note that we create the sampled_points using the frequencies and power densities defined above.

dephasing_noise_spectral_density = qctrl.factories.piecewise_linear_noise_spectral_densities.create(
    sampled_points=[{'frequency': f,
                     'power': p,
                     'power_uncertainty': 0.,
                     'weight': 0.} for f, p in zip(frequencies, power_densities)],

Creating the target (optional)

Simulations can optionally be provided with a target time evolution operator, which will cause operational infidelities to be calculated throughout the simulated time period. If provided, the target is represented as a target object, which requires a system, unitary operator and projection operator (which describes the subspace of interest, in this case the full Hilbert space).

target = qctrl.factories.targets.create(

Creating the simulation

Computation of a simulation with colored noise is represented as a colored_noise_simulation object, which contains options for the simulation:

  • point_times, which is a list of sample times at which the dynamics should be calculated (the final time must match the duration of the control pulse),
  • trajectory_count, which is the number of independent simulations to run using different random realizations of the noise processes,
  • initial_state_vector, which gives an initial state vector to be propagated according to the calculated dynamics.
simulation = qctrl.factories.colored_noise_simulations.create(
    point_times=np.linspace(0, pulse.duration, 100),
    initial_state_vector=np.array([1., 0.]),

Running the simulation

With the colored_noise_simulation object prepared, all that remains is to perform the actual simulation. This may be accomplished using the colored_noise_simulations service, as shown below.

simulation_result =
100%|██████████| 100/100 [00:18<00:00,  5.47it/s]

Extracting the average infidelities

A colored noise simulation produces many different types of data. The first we consider is the average infidelity, which gives the operational infidelity of the time evolution operator relative to the specified target, averaged over the ensemble of trajectories. All "average" quantities are available in the average_frames field of the simulation.

# Extract infidelities
average_infidelities = [frame['infidelity'] for frame in simulation_result.simulations[0].average_frames]

# Plot infidelities
plt.plot(simulation_result.simulations[0].point_times*1e6, average_infidelities)
plt.xlabel("Time ($\mu$s)")
plt.title("Average infidelity")

Extracting the average state matrices

The object describing an ensemble of state vectors is a density matrix. The density matrix at each time step can also be extracted from the average frame, and used (for example) to visualize trajectories of observables.

# Extract density matrices
density_matrices = [frame['state_matrix'] for frame in simulation_result.simulations[0].average_frames]

# Calculate and plot average sigma_x expectation
          for density_matrix in density_matrices])
plt.axhline(y=0., c='k', lw=0.5)
plt.xlabel("Time ($\mu$s)")
plt.title(r"Average $\sigma_x$ expectation")

Extracting the individual noise trajectories

In addition to the ensemble information available in the average frames, we can extract information about the individual simulation runs. Here we show how to extract the random realizations of the dephasing noise process.

# Extract noise trajectories
noise_samples_list = [trajectory.noise_samples[0]['sampled_points']
                      for trajectory in simulation_result.simulations[0].trajectories]

# Plot noise trajectories
for noise_samples in noise_samples_list:
    plt.plot([noise_sample['time']*1e6 for noise_sample in noise_samples],
             [noise_sample['value']/1e6 for noise_sample in noise_samples])
plt.xlabel("Time ($\mu$s)")
plt.ylabel("Noise value (MHz)")
plt.title("Noise trajectories")

Extracting the individual time evolution trajectories

Next we show how to extract the trajectories of the unitary time evolution operators, for the individual simulation runs.

# Extract first two time evolution operator trajectories
time_evolution_operators_list = [[frame['unitary_operator'] for frame in trajectory.frames]
                                 for trajectory in simulation_result.simulations[0].trajectories[0:2]]

# Print initial and final time evolution operators, for both trajectories
print("Time evolution operators at start:")
print("Time evolution operators at end:")
Time evolution operators at start:
[[1.+0.j 0.+0.j]
 [0.+0.j 1.+0.j]]
[[1.+0.j 0.+0.j]
 [0.+0.j 1.+0.j]]

Time evolution operators at end:
[[ 0.70543832+0.02355827j  0.00322172-0.70837236j]
 [-0.00322172-0.70837236j  0.70543832-0.02355827j]]
[[ 0.70726153+0.00854805j  0.03093181-0.70622325j]
 [-0.03093181-0.70622325j  0.70726153-0.00854805j]]


We have shown how to use the Q-CTRL Python package to simulate the dynamics of a system in the presence of a stochastic noise process. We have also shown how to extract several quantities of interest from the simulation result, including both ensemble objects describing the average system behavior and objects associated with specific realizations of the stochastic noise process.

Example: single qubit with leakage

In this example we show how to simulate the dynamics of a single qubit subject to leakage to a third energy level (i.e. a qutrit), driven by a primitive $\pi$ pulse. The qutrit is treated as an oscillator (truncated to three levels) with an anharmonicity of $\chi$, described by the Hamiltonian:

\begin{align*} H(t) = & \frac{\chi}{2} (a^\dagger)^2 a^2 + \frac{\Omega(t)}{2} a + \frac{\Omega^*(t)}{2} a^\dagger, \end{align*}

where $a = \left|0 \right\rangle \left\langle 1 \right| + \sqrt{2} \left|1 \right\rangle \left\langle 2 \right|$ is the lowering operator and $\Omega(t)$ is a time-dependent Rabi rate.

Below we show how to use the Q-CTRL Python package to simulate the time evolution of this coherent system.

# Define standard matrices
a = np.array([[0., 1., 0.],
              [0., 0., np.sqrt(2)],
              [0., 0., 0.]], dtype=np.complex)
ad2a2 =

# Define system parameters
chi = 2*np.pi * 3 * 1e6 #Hz
omega_max = 2*np.pi * 1e6 #Hz
total_rotation = np.pi

# Define system object
system =
    name='Single qubit with leakage',

# Define control object
drive = qctrl.factories.drive_controls.create(
    name='Rabi rate',

# Define drift object
anharmonicity = qctrl.factories.drift_controls.create(

# Define pulse object using pulses from Q-CTRL Open Controls
pulse = new_predefined_driven_control(

drive_pulse = qctrl.factories.custom_pulses.create(
    segments=[{'duration': d, 'value': v} 
              for d, v in zip(pulse.durations,
                              pulse.rabi_rates * np.exp(1j*pulse.azimuthal_angles))],

# Define simulation object
simulation = qctrl.factories.coherent_simulations.create(
    point_times=np.linspace(0, pulse.duration, 100),
    initial_state_vector=np.array([1., 0., 0.]),

# Run simulation
simulation_result =

# Extract and print final time evolution operator
print("Time evolution operator at end:")

# Extract and plot state populations
times = np.array([frame['time']
                  for frame in simulation_result.simulations[0].trajectories[0].frames])
state_vectors = np.array([frame['state_vector']
                          for frame in simulation_result.simulations[0].trajectories[0].frames])
for state in range(3):
    plt.plot(times*1e6, np.abs(state_vectors[:, state])**2, label=f"P{state}")
plt.xlabel("Time ($\mu$s)")
plt.title("State populations")
100%|██████████| 100/100 [00:06<00:00, 14.94it/s]
Time evolution operator at end:
[[ 0.05730981-0.14539786j  0.23505273-0.92588994j -0.08764586+0.23531105j]
 [ 0.23505273-0.92588994j -0.06664016+0.18738223j -0.19346042+0.10246022j]
 [-0.08764586+0.23531105j -0.19346042+0.10246022j -0.82544823+0.45569409j]]


Comprehensive knowledge base of quantum control theory