{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# How to simulate open system dynamics\n",
"**Calculating the dynamics of a quantum system described by a GKS–Lindblad master equation**\n",
"\n",
"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)](https://doi.org/10.1063/1.522979) and [Lindblad](https://doi.org/10.1007/BF01608499):\n",
"\n",
"$$ \\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). $$\n",
"\n",
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Summary workflow\n",
"\n",
"### 1. Define Lindblad terms in computational graph\n",
"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](https://docs.q-ctrl.com/boulder-opal/user-guides/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.\n",
"\n",
"\n",
"### 2. Calculate density matrix evolution\n",
"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.\n",
"\n",
"**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](https://docs.q-ctrl.com/boulder-opal/user-guides/how-to-simulate-large-open-system-dynamics) covers efficient approaches for large systems.\n",
"\n",
"## Example: Simulating a single-qubit open system\n",
"\n",
"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$:\n",
"\n",
"$$ \\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]. $$\n",
"\n",
"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:\n",
"\n",
"$$ H_{\\rm s}(t) = \\frac{1}{2} \\left( \\Omega(t) \\sigma_- + \\Omega^* (t) \\sigma_+ \\right) + \\Delta(t) \\sigma_z, $$\n",
"\n",
"where $\\sigma_x$, $\\sigma_y$, and $\\sigma_z$ are the Pauli matrices and $\\sigma_\\pm = (\\sigma_x \\mp i \\sigma_y)/2$.\n",
"\n",
"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)$.\n",
"\n",
"\n",
"\n",
"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](https://docs.q-ctrl.com/open-controls/references/qctrl-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$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from qctrlopencontrols import new_primitive_control\n",
"from qctrlvisualizer import get_qctrl_style, plot_population_dynamics\n",
"\n",
"from qctrl import Qctrl\n",
"\n",
"plt.style.use(get_qctrl_style())\n",
"\n",
"# Start a Boulder Opal session.\n",
"qctrl = Qctrl()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Define control parameters.\n",
"omega_max = 2 * np.pi * 1e6 # Hz\n",
"total_rotation = np.pi\n",
"\n",
"# Obtain predefined pulse from Q-CTRL Open Controls.\n",
"predefined_pulse = new_primitive_control(\n",
" rabi_rotation=total_rotation, azimuthal_angle=0.0, maximum_rabi_rate=omega_max\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Defining the simulation graph\n",
" 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`. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/100 [00:00"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"closed_system_populations = (\n",
" np.abs(results.output[\"closed_system_states\"][\"value\"][:, :, 0]) ** 2\n",
")\n",
"open_system_populations = np.real_if_close(\n",
" np.diagonal(results.output[\"density_matrices\"][\"value\"], axis1=-1, axis2=-2)\n",
")\n",
"\n",
"plot_population_dynamics(\n",
" sample_times,\n",
" {\n",
" r\"$|{0}\\rangle$ (closed system)\": closed_system_populations[:, 0],\n",
" r\"$|{1}\\rangle$ (closed system)\": closed_system_populations[:, 1],\n",
" r\"$|{0}\\rangle$ (open system)\": open_system_populations[:, 0],\n",
" r\"$|{1}\\rangle$ (open system)\": open_system_populations[:, 1],\n",
" },\n",
")\n",
"plt.suptitle(\n",
" f\"Dynamics comparison between closed and open system ($T_2$ = {T2*1e6} μs)\"\n",
")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Summary\n",
"\n",
"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](https://docs.q-ctrl.com/boulder-opal/user-guides/how-to-simulate-quantum-dynamics-subject-to-noise-with-graphs) 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."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"| Package | Version |\n",
"| --------------------- | ------------ |\n",
"| Python | 3.10.8 |\n",
"| matplotlib | 3.6.3 |\n",
"| numpy | 1.23.1 |\n",
"| scipy | 1.10.0 |\n",
"| qctrl | 20.2.0 |\n",
"| qctrl-commons | 17.10.0 |\n",
"| qctrl-open-controls | 9.1.7 |\n",
"| boulder-opal-toolkits | 2.0.0-beta.4 |\n",
"| qctrl-visualizer | 4.6.0 |\n"
]
}
],
"source": [
"from qctrl.utils import print_environment_related_packages\n",
"\n",
"print_environment_related_packages()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.8"
}
},
"nbformat": 4,
"nbformat_minor": 4
}