How to simulate open system dynamics

Boulder Opal enables you to simulate not only the evolution of isolated quantum systems with and without noise, but also the evolution of systems that are interacting with their environment. 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:

$$ \frac{\mathrm{d}}{\mathrm{d} t} \rho(t) = - i \left[ H_{\rm s} (t), \rho(t) \right] + \gamma L \rho(t) L^\dagger - \frac{1}{2} \gamma \rho(t) L^\dagger L - \frac{1}{2} \gamma L^\dagger L \rho(t). $$

The tools described in this notebook allow you to solve this GKS–Lindblad master equation to obtain the time evolution for the system that you want.

Summary workflow

1. Define Lindblad terms in computational graph

To describe the dynamics of an open system with Boulder Opal, you start by setting up a graph object as described in the How to represent quantum systems using graphs user guide. The difference is that besides setting up a graph with a piecewise-constant system Hamiltonian $H_{\rm s} (t)$ for the controls, you must also define the Lindblad terms that appear in the GKS–Lindblad master equation above.

2. Calculate density matrix evolution

If you provide a tuple with the rate $\gamma$ and Lindblad operator $L$ as lindblad_terms to the graph operation graph.density_matrix_evolution_pwc, together with a piecewise constant hamiltonian and an initial_density_matrix, you obtain the value of the density matrix $\rho(t)$ for the sample_times that you requested in a call of the function qctrl.functions.calculate_graph. It also performs the same operation for a closed system (which takes into account the system Hamiltonian, but not the Lindbladian), so you can compare the two types of evolution.

Note that the density_matrix_evolution_pwc function can be inefficient when used with the default parameters if the dimension of your Hilbert space is more than roughly $10$. A separate user guide covers efficient approaches for large systems.

Example: Simulating a single-qubit open system

The master equation integration tool from Boulder Opal allows you to obtain the solution of a GKS–Lindblad equation. This equation describes the evolution of a system that obeys a system Hamiltonian $H_{\rm s}$ and also any number of extra Lindbladian terms consisting of a rate $\gamma_i$ and a Lindblad operator $L_i$:

$$ \frac{\mathrm{d} \rho(t)}{\mathrm{d} t} = - i \left[ H_{\rm s}(t), \rho(t) \right] + \sum_i \gamma_i \left[ L_i \rho(t) L_i^\dagger - \frac{1}{2} \rho(t) L_i^\dagger L_i - \frac{1}{2} L_i^\dagger L_i \rho(t) \right]. $$

For this example, consider a single-qubit system whose evolution you can control through a complex Rabi rate $\Omega(t)$ and a detuning $\Delta(t)$, so that the system Hamiltonian is:

$$ H_{\rm s}(t) = \frac{1}{2} \left( \Omega(t) \sigma_- + \Omega^* (t) \sigma_+ \right) + \Delta(t) \sigma_z, $$

where $\sigma_x$, $\sigma_y$, and $\sigma_z$ are the Pauli matrices and $\sigma_\pm = (\sigma_x \mp i \sigma_y)/2$.

Suppose this system is also subject to $T_2$ noise, which you can represent by the Lindblad operator $L= \sigma_z$ with associated rate $\gamma = 1/(2T_2)$. By using the tools from Boulder Opal to solve this master equation, you can obtain the evolution of the system in terms of $\rho(t)$.

The first step is to set up the Python object that represents the pulse, which you can define by yourself or by importing one of the predefined pulses in the Q-CTRL Open Controls package. In particular, this example illustrates the case where the target operation is an $X_{\pi}$ rotation applied to a system initially at the state $\left|0\right\rangle$, which corresponds to the initial density matrix $\rho(0) = \left| 0 \right\rangle \left\langle 0 \right|$. Such an operation should flip the qubit from the state $\left|0\right\rangle$ to the state $\left|1\right\rangle$.

import matplotlib.pyplot as plt
import numpy as np
from qctrlopencontrols import new_primitive_control
from qctrlvisualizer import get_qctrl_style, plot_population_dynamics

from qctrl import Qctrl

plt.style.use(get_qctrl_style())

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

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

# Obtain predefined pulse from Q-CTRL Open Controls.
predefined_pulse = new_primitive_control(
    rabi_rotation=total_rotation, azimuthal_angle=0.0, maximum_rabi_rate=omega_max
)

Defining the simulation graph

The code block below illustrates how to obtain the value of the density matrix $\rho(t)$ for the sample_times that you requested in a call of the function qctrl.functions.calculate_graph.

# Define time parameters of the simulation.
duration = sum(predefined_pulse.durations)
sample_times = np.linspace(0, duration, 100)

# Create simulation graph.
graph = qctrl.create_graph()

initial_state = graph.fock_state(2, 0)
initial_density_matrix = graph.outer_product(initial_state, initial_state)

# Define the controls using the predefined pulse.
omega = graph.pwc(
    values=predefined_pulse.rabi_rates * np.exp(1j * predefined_pulse.azimuthal_angles),
    durations=predefined_pulse.durations,
)
delta = graph.pwc(
    values=predefined_pulse.detunings, durations=predefined_pulse.durations
)

# Define terms of the Hamiltonian.
sigma_x_term = graph.hermitian_part(omega * graph.pauli_matrix("M"))
sigma_z_term = delta * graph.pauli_matrix("Z")
hamiltonian = sigma_x_term + sigma_z_term

# Calculate vector state evolution for a closed system.
graph.matmul(
    graph.time_evolution_operators_pwc(hamiltonian, sample_times),
    initial_state[:, None],
    name="closed_system_states",
)

# Define Lindblad term of the master equation.
lindblad_operator = graph.pauli_matrix("Z")
T2 = 1e-6  # s

# Calculate density matrix evolution according to the master equation.
graph.density_matrix_evolution_pwc(
    initial_density_matrix=initial_density_matrix,
    hamiltonian=hamiltonian,
    lindblad_terms=[(1 / (2 * T2), lindblad_operator)],
    sample_times=sample_times,
    name="density_matrices",
)

# Run simulation.
results = qctrl.functions.calculate_graph(
    graph=graph, output_node_names=["closed_system_states", "density_matrices"]
)

Your task calculate_graph (action_id="1618605") has started.
Your task calculate_graph (action_id="1618605") has completed.

Plotting the system evolution

After you obtain the evolution of the system with and without $T_2$ noise, you can plot the results at each of the sample_times that you requested to see how the system goes from the initial state $\left|0\right\rangle$ to the target state $\left| 1 \right\rangle$. While the evolution for a closed system leads perfectly to the desired final state, the open system evolution is disrupted by the presence of $T_2$ noise.

closed_system_populations = (
    np.abs(results.output["closed_system_states"]["value"][:, :, 0]) ** 2
)
open_system_populations = np.real_if_close(
    np.diagonal(results.output["density_matrices"]["value"], axis1=-1, axis2=-2)
)

plot_population_dynamics(
    sample_times,
    {
        r"$|{0}\rangle$ (closed system)": closed_system_populations[:, 0],
        r"$|{1}\rangle$ (closed system)": closed_system_populations[:, 1],
        r"$|{0}\rangle$ (open system)": open_system_populations[:, 0],
        r"$|{1}\rangle$ (open system)": open_system_populations[:, 1],
    },
)
plt.suptitle(
    f"Dynamics comparison between closed and open system ($T_2$ = {T2*1e6} μs)"
)
plt.show()

Summary

The open-system tools from Boulder Opal allow you to use the graph framework to calculate more general types of system evolution. While the tools described in this user guide show how you can obtain the evolution of a system that follows the Schrödinger equation, the open system tools allow you to solve a more general master equation. You can also use these tools in optimization graphs, if you have a cost function that represents what you want to optimize.

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