Simulate open system dynamics

Calculating the dynamics of a quantum system described by a GKS–Lindblad master equation

The Q-CTRL Python package enables you to simulate not only the evolution of isolated quantum systems (see the Simulate quantum dynamics user guide), 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. The tools described in this user guide allow you to solve this GKS–Lindblad master equation to obtain the time evolution for the system that you want.

Import 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 controls
from qctrlopencontrols import (
    new_corpse_control,
    new_primitive_control,
)

# Plotting imports
import matplotlib.pyplot as plt
from qctrlvisualizer import (
    get_qctrl_style,
    plot_controls,
)

plt.style.use(get_qctrl_style())

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

Worked example: Simulating a single-qubit open system

The master equation integration tool from the Q-CTRL Python package 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} \Omega(t) \sigma_- + \frac{1}{2} \Omega^* (t) \sigma_+ + \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 the Q-CTRL Python package to solve this master equation, you can obtain the evolution of the system in terms of $\rho(t)$.

Creating the pulse

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$.

# Define standard matrices.
identity = np.identity(2, dtype=np.complex)
sigma_x = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex)
sigma_z = np.array([[1.0, 0.0], [0.0, -1.0]], dtype=np.complex)
sigma_m = np.array([[0.0, 1.0], [0.0, 0.0]], dtype=np.complex)

# Define control parameters.
omega_max = 2 * np.pi * 1e6  # Hz
total_rotation = np.pi
initial_state = np.array([[1], [0]], dtype=np.complex)
initial_density_matrix = np.kron(initial_state, initial_state.T.conj())

# 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,
)

Creating the graph

To describe the dynamics of an open system with the Q-CTRL Python package, you start by setting up a graph object as described in the Set up quantum systems 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:

$$ \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). $$

If you provide a tuple with the rate $\gamma$ and Lindblad operator $L$ as lindblad_terms to the operation qctrl.operations.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. The code block below illustrates how you can do that 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.

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

# Set up simulation graph.
with qctrl.create_graph() as graph:
    # Define the controls using the predefined pulse.
    omega = qctrl.operations.tensor_pwc(
        values=predefined_pulse.rabi_rates
        * np.exp(1j * predefined_pulse.azimuthal_angles),
        durations=predefined_pulse.durations,
    )
    delta = qctrl.operations.tensor_pwc(
        values=predefined_pulse.detunings, durations=predefined_pulse.durations
    )

    # Define terms of the Hamiltonian.
    sigma_m_term = qctrl.operations.pwc_operator(signal=omega, operator=sigma_m)
    sigma_z_term = qctrl.operations.pwc_operator(signal=delta, operator=sigma_z)
    sigma_x_term = qctrl.operations.pwc_operator_hermitian_part(operator=sigma_m_term)
    hamiltonian = qctrl.operations.pwc_sum([sigma_x_term, sigma_z_term])

    # Define Lindblad term of the master equation.
    lindblad_operator = sigma_z
    T2 = 1e-6  # s

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

    # Calculate vector state evolution for a closed system.
    qctrl.operations.matmul(
        qctrl.operations.time_evolution_operators_pwc(
            hamiltonian=hamiltonian, sample_times=sample_times,
        ),
        initial_state,
        name="closed_system_states",
    )

# Run simulation.
results = qctrl.functions.calculate_graph(
    graph=graph, output_node_names=["closed_system_states", "open_system_states"],
)
Your task calculate_graph has completed in 3s.

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.

plt.figure()
plt.xlabel(f"Time ($\mu$s)")
plt.ylabel(r"Population of $| 0 \rangle$")
plt.plot(
    sample_times * 1e6,
    np.abs([density_matrix[0][0] for density_matrix in results.output["open_system_states"]["value"]]),
    label=f"Open system with $T_2$ = {T2*1e6} μs",
)
plt.plot(
    sample_times * 1e6,
    np.abs([state_vector[0][0] for state_vector in results.output["closed_system_states"]["value"]])
    ** 2,
    label="Closed system",
)
plt.legend()

plt.figure()
plt.xlabel(f"Time ($\mu$s)")
plt.ylabel(r"Population of $| 1 \rangle$")
plt.plot(
    sample_times * 1e6,
    np.abs([density_matrix[1][1] for density_matrix in results.output["open_system_states"]["value"]]),
    label=f"Open system with $T_2$ = {T2*1e6} μs",
)
plt.plot(
    sample_times * 1e6,
    np.abs([state_vector[1][0] for state_vector in results.output["closed_system_states"]["value"]])
    ** 2,
    label="Closed system",
)
plt.legend()

plt.show()

Summary

The open-system tools from the Q-CTRL Python package allow you to use the graph framework to calculate more general types of system evolution. While the tools described in the Simulate quantum dynamics 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.

Example: Quasi-static scan of a system exposed to $T_1$ noise

The quasi-static scans described in the Quasi-static error susceptibility user guide are useful to compare the performance of pulses against noise of different strengths. That user guide explains how you can calculate the quasi-static scan for a closed quantum system that evolves according to the Schrödinger equation. To create quasi-static scans for a system that evolves according to a master equation, you need to employ the graph framework to obtain the points that compose the quasi-static scan.

For this example, consider a system under the following system Hamiltonian:

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

where $\sigma_x$, $\sigma_y$, and $\sigma_z$ are the Pauli matrices and $\sigma_\pm = (\sigma_x \mp i \sigma_y)/2$. The signal $\Omega(t)$ represents your complex-valued controls, while $\eta$ is the noise coefficient. Assuming this noise changes sufficiently slowly that it is constant in each run of the gate, each point of the quasi-static scan represents the performance of the pulse under different strengths of $\eta$.

To quantify the performance of the pulse, you can use the state infidelity as a metric. Using the fact that the target state is a pure state $\rho_{\rm target}$, you can find the infidelity of a pulse that evolves the system to a state $\rho(t)$ through the formula:

$$ \mathcal{I} (t) = 1 - \left( \mathrm{Tr} \left\{ \sqrt{ \sqrt{\rho_{\rm target}} \rho(t) \sqrt{\rho_{\rm target}} } \right\} \right)^2 = 1 - \mathrm{Tr} \left\{ \rho_{\rm target} \rho(t) \right\}. $$

Consider as well that the system can decay at a certain rate $T_1$, so that the complete GKS–Lindblad master equation that describes its evolution is:

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

In this master equation, $\sigma_-$ is the Lindblad operator with an associated rate $1/T_1$.

To see how a set of pulses performs against additive $\sigma_z$ noise, when also subject to $T_1$ noise, you can use the open system tools to calculate the evolution of the master equation for different values of quasi-static noise $\eta$. This example shows how you can create an optimized pulse under these conditions by setting the sum of the infidelities under different quasi-static strengths as the cost function of the optimization. The code block below compares the performance of the predefined pulse and the optimized pulse, plotting the shape of the controls and the quasi-static scans for both of them.

# Define standard matrices.
identity = np.identity(2, dtype=np.complex)
sigma_x = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex)
sigma_z = np.array([[1.0, 0.0], [0.0, -1.0]], dtype=np.complex)
sigma_m = np.array([[0.0, 1.0], [0.0, 0.0]], dtype=np.complex)

# Define noise and control parameters.
T1 = 100e-6  # s
omega_max = 2 * np.pi * 1e6  # Hz
total_rotation = np.pi / 2
target_operator = (identity - 1j * sigma_x) / np.sqrt(2)
initial_density_matrix = np.array([[1, 0], [0, 0]], dtype=np.complex)
target_density_matrix = (
    target_operator @ initial_density_matrix @ target_operator.T.conj()
)

# Define scan parameters
scan_point_count = 20
scan_range = np.linspace(-omega_max / 2, omega_max / 2, scan_point_count)

# Obtain predefined pulse from Q-CTRL Open Controls.
corpse_pulse = new_corpse_control(
    rabi_rotation=total_rotation, azimuthal_angle=0.0, maximum_rabi_rate=omega_max,
)
duration = sum(corpse_pulse.durations)


def state_infidelity(density_matrix, name):
    """
    Calculates the state infidelity of the `density_matrix` with respect to the
    `target_density_matrix`, which represents a pure state.
    """
    return qctrl.operations.abs(
        1
        - qctrl.operations.trace(
            qctrl.operations.matmul(target_density_matrix, density_matrix)
        ),
        name=name,
    )


# Set up optimization graph.
with qctrl.create_graph() as graph:
    # Define the controls for the predefined pulse.
    corpse_omega = qctrl.operations.tensor_pwc(
        values=corpse_pulse.rabi_rates
        * np.exp(1j * corpse_pulse.azimuthal_angles),
        durations=corpse_pulse.durations,
        name="corpse_pulse",
    )

    # Define controls for the optimized pulse.
    optimized_omega = qctrl.operations.complex_pwc_signal(
        moduli=qctrl.operations.bounded_optimization_variable(10, 0, omega_max),
        phases=qctrl.operations.unbounded_optimization_variable(10, 0, 2 * np.pi),
        duration=duration,
        name="optimized_pulse",
    )

    # Create control Hamiltonians for the predefined pulse and the optimized pulse.
    noiseless_corpse_hamiltonian = qctrl.operations.pwc_operator_hermitian_part(
        operator=qctrl.operations.pwc_operator(
            signal=corpse_omega, operator=sigma_m
        )
    )
    noiseless_optimized_hamiltonian = qctrl.operations.pwc_operator_hermitian_part(
        operator=qctrl.operations.pwc_operator(signal=optimized_omega, operator=sigma_m)
    )

    # Create list of detuning terms for the quasi-static scan.
    scan_terms = [
        qctrl.operations.constant_pwc_operator(
            operator=amplitude * sigma_z, duration=duration
        )
        for amplitude in scan_range
    ]

    # Add detuning terms to the control Hamiltonians.
    corpse_hamiltonians = [
        qctrl.operations.pwc_sum([noiseless_corpse_hamiltonian, scan_term])
        for scan_term in scan_terms
    ]
    optimized_hamiltonians = [
        qctrl.operations.pwc_sum([noiseless_optimized_hamiltonian, scan_term])
        for scan_term in scan_terms
    ]

    # Define decay term of the master equation.
    lindblad_operator = sigma_m

    # Calculate the final density matrices from the master equation.
    corpse_density_matrices = [
        qctrl.operations.density_matrix_evolution_pwc(
            initial_density_matrix=initial_density_matrix,
            hamiltonian=corpse_hamiltonian,
            lindblad_terms=[(1 / T1, lindblad_operator)],
            sample_times=np.array([duration]),
        )[-1]
        for corpse_hamiltonian in corpse_hamiltonians
    ]
    optimized_density_matrices = [
        qctrl.operations.density_matrix_evolution_pwc(
            initial_density_matrix=initial_density_matrix,
            hamiltonian=optimized_hamiltonian,
            lindblad_terms=[(1 / T1, lindblad_operator)],
            sample_times=np.array([duration]),
        )[-1]
        for optimized_hamiltonian in optimized_hamiltonians
    ]

    # Calculate infidelities with respect to the target.
    corpse_infidelities = [
        state_infidelity(
            density_matrix=density_matrix, name=f"corpse_infidelity{number}",
        )
        for number, density_matrix in enumerate(corpse_density_matrices)
    ]
    optimized_infidelities = [
        state_infidelity(
            density_matrix=density_matrix, name=f"optimized_infidelity{number}",
        )
        for number, density_matrix in enumerate(optimized_density_matrices)
    ]

    # Create cost as the sum of all the infidelities of the optimized pulse.
    qctrl.operations.sum(optimized_infidelities, name="cost")

# Run optimization.
results = qctrl.functions.calculate_optimization(
    graph=graph,
    cost_node_name="cost",
    output_node_names=["corpse_pulse", "optimized_pulse"]
    + [f"corpse_infidelity{number}" for number in range(scan_point_count)]
    + [f"optimized_infidelity{number}" for number in range(scan_point_count)],
)

# Plot pulses.
fig = plt.figure()
plt.title("CORPSE pulse")
plt.xticks(([0]))
plt.yticks(([0]))
plot_controls(fig, {r"$\Omega$": results.output[f"corpse_pulse"]})
fig2 = plt.figure()
plt.title("Optimized pulse")
plt.xticks(([0]))
plt.yticks(([0]))
plot_controls(fig2, {r"$\Omega$": results.output[f"optimized_pulse"]})

# Plot quasi-static scan.
fig3 = plt.figure()
plt.xlabel(r"Detuning $\eta$/(2$\pi$) (MHz)")
plt.ylabel("Infidelity")
plt.plot(
    scan_range * 1e-6 / (2 * np.pi),
    np.array(
        [
            results.output[f"corpse_infidelity{number}"]["value"]
            for number in range(scan_point_count)
        ]
    ),
    label="CORPSE pulse",
)
plt.plot(
    scan_range * 1e-6 / (2 * np.pi),
    [
        results.output[f"optimized_infidelity{number}"]["value"]
        for number in range(scan_point_count)
    ],
    label="Optimized pulse",
)
plt.legend()
plt.show()
Your task calculate_optimization has started.
Your task calculate_optimization has completed in 199s.