{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "94786846",
      "metadata": {
        "id": "94786846"
      },
      "source": [
        "# Solving the Lindblad Master equation with QuTiP"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Theory"
      ],
      "metadata": {
        "id": "O29HWycvQQ7s"
      },
      "id": "O29HWycvQQ7s"
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Dynamics: \n",
        "\n",
        "$$\n",
        "\\hat{H}_{tot} = \\hat{H}_{sys} + \\hat{H}_{env} + \\hat{H}_{int}\n",
        "$$\n",
        "\n",
        "\n",
        "---\n",
        "\n",
        "### Liouville - von Neumann equation:\n",
        "\n",
        "$$\n",
        "\\dot{\\rho}_{tot}(t) = -\\frac{i}{\\hbar} [\\hat{H}_{tot}(t), \\rho_{tot}(t)]\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\rho_{sys}(t) = \\text{Tr}_{env}[\\rho_{tot}(t)]\n",
        "$$\n",
        "\n",
        "\n",
        "---\n",
        "\n",
        "\n",
        "### Assumptions:\n",
        "\n",
        "##### 1. Environment is in thermal equilibrium at temperature $T.$ \n",
        "\n",
        "##### 2. At initial time $t_{0},$ system and environment are weakly coupled i.e \n",
        "\n",
        "$$\n",
        "\\rho_{tot}(0) = \\rho_{sys}(0) \\otimes \\rho_{env}\n",
        "$$\n",
        "\n",
        "##### 3. Born approximation: The system and environment remains separable i.e\n",
        "\n",
        "$$\n",
        "\\rho_{tot}(t) \\approx \\rho_{sys}(t) \\otimes \\rho_{env}\n",
        "$$\n",
        "\n",
        "##### 4. Markov approximation: The state of the system at $t$ only depends the state of the system at $t + \\Delta t.$\n",
        "\n",
        "##### 5. Secular or rotating wave approximation.\n",
        "\n",
        "\n",
        "---\n",
        "\n",
        "\n",
        "### Gorini - Kossakowski - Sudarshan - Lindblad Master equation (written in QuTiP notation):\n",
        "$$\n",
        "\\dot{\\rho}_{sys}(t) = -\\frac{i}{\\hbar} [ \\hat{H}_{sys}(t), \\rho_{sys}(t) ] + \\sum_{n} \\left( \\hat{C}_{n} \\rho_{sys}(t) \\hat{C}^{\\dagger}_{n} - \\frac{1}{2} \\biggl\\{ \\hat{C}^{\\dagger}_{n} \\hat{C}_{n}, \\rho_{sys}(t) \\biggl\\} \\right)\n",
        "$$\n",
        "\n",
        "where \n",
        "$$\n",
        "\\hat{C}_{n} = \\sqrt{\\gamma_{n}} \\hat{A}_{n}\n",
        "$$\n",
        "\n",
        "$\\gamma_{n}:$ corresponding rate. \n",
        "\n",
        "$\\hat{C}_{n}:$ collaspe operator (QuTiP terminology).\n",
        "\n",
        "$\\hat{A}_{n}:$ dissipation operator. "
      ],
      "metadata": {
        "id": "QL0Xk_Q9Ero8"
      },
      "id": "QL0Xk_Q9Ero8"
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Example: Two-level system decay and driven by an oscillating field"
      ],
      "metadata": {
        "id": "B9AZa8BaaV0U"
      },
      "id": "B9AZa8BaaV0U"
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Pen & paper derivation"
      ],
      "metadata": {
        "id": "4HIgDxo2b90V"
      },
      "id": "4HIgDxo2b90V"
    },
    {
      "cell_type": "markdown",
      "source": [
        "Density matrix:\n",
        "$$\n",
        "| 0 \\rangle \\equiv \n",
        "\\begin{bmatrix} \n",
        "1 \\\\ 0 \n",
        "\\end{bmatrix}, \\\\\n",
        "| 1 \\rangle \\equiv \n",
        "\\begin{bmatrix} \n",
        "0 \\\\ 1 \n",
        "\\end{bmatrix}, \\\\\n",
        "\\rho_{sys}(t) = p_{00}(t) |0 \\rangle \\langle 0| + p_{01}(t) |0 \\rangle \\langle 1| + p_{10}(t) |1 \\rangle \\langle 0| + p_{11}(t) |1 \\rangle \\langle 1|, \\\\\n",
        "\\rho_{sys}(t) = \n",
        "\\begin{bmatrix} \n",
        "p_{00}(t) & p_{01}(t) \\\\ p_{10}(t) & p_{11}(t)\n",
        "\\end{bmatrix}.\n",
        "$$\n",
        "\n",
        "Driving frequency: $\\Omega.$\n",
        "\n",
        "System Hamiltonian: \n",
        "$$\n",
        "\\hat{H}_{sys}(t) = E_{0} |0 \\rangle \\langle 0| + E_{1} |1 \\rangle \\langle 1| + \\hbar \\Omega |0 \\rangle \\langle 1| + \\hbar \\Omega |1 \\rangle \\langle 0|.\n",
        "$$\n",
        "\n",
        "Change energy reference to $E_{0} = 0$ and let $E_{1} = E.$ Also, for simplicity, assume $\\hbar = 1.$ \n",
        "\n",
        "Then,\n",
        "$$\n",
        "\\hat{H}_{sys}(t) = E |1 \\rangle \\langle 1| + \\Omega \\big[ |0 \\rangle \\langle 1| + |1 \\rangle \\langle 0| \\big], \\\\\n",
        "\\hat{H}_{sys}(t) = E \n",
        "\\begin{bmatrix} \n",
        "0 & 0 \\\\ 0 & 1\n",
        "\\end{bmatrix} + \\Omega \n",
        "\\begin{bmatrix} \n",
        "0 & 1 \\\\ 1 & 0\n",
        "\\end{bmatrix}.\n",
        "$$\n",
        "\n",
        "\n",
        "---\n",
        "\n",
        "\n",
        "Decay rate: $\\beta.$\n",
        "\n",
        "Dissipation operator: \n",
        "$$\n",
        "\\hat{\\sigma}^{-} = | 0 \\rangle \\langle 1 | = \n",
        "\\begin{bmatrix} \n",
        "1 \\\\ 0 \n",
        "\\end{bmatrix} [0 \\quad 1] = \n",
        "\\begin{bmatrix} \n",
        "0 & 1 \\\\ 0 & 0 \n",
        "\\end{bmatrix}, \\\\\n",
        "(\\hat{\\sigma}^{-})^{\\dagger} = \\hat{\\sigma}^{+} = \n",
        "\\begin{bmatrix} \n",
        "0 & 0 \\\\ 1 & 0 \n",
        "\\end{bmatrix}.\n",
        "$$\n",
        "\n",
        "\n",
        "---\n",
        "\n",
        "\n",
        "Recall the Lindblad Master equation (in QuTiP notation): \n",
        "$$\n",
        "\\dot{\\rho}_{sys}(t) = -\\frac{i}{\\hbar} \\biggl[ \\hat{H}_{sys}(t), \\rho_{sys}(t) \\biggl] + \\sum_{n} \\left( \\hat{C}_{n} \\rho_{sys}(t) \\hat{C}^{\\dagger}_{n} - \\frac{1}{2} \\biggl\\{ \\hat{C}^{\\dagger}_{n} \\hat{C}_{n}, \\rho_{sys}(t) \\biggl\\} \\right), \\\\\n",
        "\\hbar = 1, \\\\\n",
        "\\dot{\\rho}_{sys}(t) = -i \\biggl[ \\hat{H}_{sys}(t), \\rho_{sys}(t) \\biggl] + \\sum_{n} \\left( \\hat{C}_{n} \\rho_{sys}(t) \\hat{C}^{\\dagger}_{n} - \\frac{1}{2} \\biggl\\{ \\hat{C}^{\\dagger}_{n} \\hat{C}_{n}, \\rho_{sys}(t) \\biggl\\} \\right).\n",
        "$$\n",
        "\n",
        "Put together terms:\n",
        "1. System Hamiltonian and density matrix:\n",
        "$$\n",
        "\\hat{H}_{sys}(t) = E \n",
        "\\begin{bmatrix} \n",
        "0 & 0 \\\\ 0 & 1\n",
        "\\end{bmatrix} + \\Omega \n",
        "\\begin{bmatrix} \n",
        "0 & 1 \\\\ 1 & 0\n",
        "\\end{bmatrix}, \\\\\n",
        "\\rho_{sys}(t) = \n",
        "\\begin{bmatrix} \n",
        "p_{00} & p_{01} \\\\ p_{10} & p_{11}\n",
        "\\end{bmatrix}.\n",
        "$$\n",
        "\n",
        "2. Commutator $-i \\biggl[ \\hat{H}_{sys}(t), \\rho_{sys}(t) \\biggl]:$\n",
        "$$\n",
        "-i \\biggl[ \\hat{H}_{sys}(t), \\rho_{sys}(t) \\biggl] = \n",
        "-i \\begin{bmatrix} \n",
        "\\Omega p_{10} - (E + \\Omega) p_{01} & \\Omega p_{11} - \\Omega p_{00} \\\\ (E + \\Omega) p_{00} - (E + \\Omega) p_{11} & (E + \\Omega) p_{01} - \\Omega p_{10} \n",
        "\\end{bmatrix}.\n",
        "$$\n",
        "\n",
        "3. Collaspe operator:\n",
        "$$\n",
        "\\hat{C}_{n} = \\sqrt{\\beta} \\hat{\\sigma}^{-} = \\sqrt{\\beta} \n",
        "\\begin{bmatrix} \n",
        "0 & 1 \\\\ 0 & 0 \n",
        "\\end{bmatrix}, \\\\\n",
        "\\hat{C}^{\\dagger}_{n} = \\sqrt{\\beta} (\\hat{\\sigma}^{-})^{\\dagger} = \\hat{\\sigma}^{+} = \\sqrt{\\beta} \n",
        "\\begin{bmatrix} \n",
        "0 & 0 \\\\ 1 & 0 \n",
        "\\end{bmatrix}.\n",
        "$$\n",
        "\n",
        "4. Term $\\hat{C}_{n} \\rho_{sys}(t) \\hat{C}^{\\dagger}_{n}$ and anti-commutator $-\\frac{1}{2} \\biggl\\{ \\hat{C}^{\\dagger}_{n} \\hat{C}_{n}, \\rho_{sys}(t) \\biggl\\}:$\n",
        "$$\n",
        "\\hat{C}_{n} \\rho_{sys}(t) \\hat{C}^{\\dagger}_{n} = \\beta\n",
        "\\begin{bmatrix} \n",
        "p_{11} & 0 \\\\ 0 & 0 \n",
        "\\end{bmatrix}, \\\\\n",
        "-\\frac{1}{2} \\biggl\\{ \\hat{C}^{\\dagger}_{n} \\hat{C}_{n}, \\rho_{sys}(t) \\biggl\\} = -\\frac{\\beta}{2} \n",
        "\\begin{bmatrix} \n",
        "0 & p_{01} \\\\ p_{10} & 2 p_{11} \n",
        "\\end{bmatrix}.\n",
        "$$\n",
        "\n",
        "\n",
        "---\n",
        "\n",
        "Lindblad Master equation of the system:\n",
        "\n",
        "$$\n",
        "\\dot{\\rho}_{sys}(t) = -i \n",
        "\\begin{bmatrix} \n",
        "\\Omega p_{10} - (E + \\Omega) p_{01} & \\Omega p_{11} - \\Omega p_{00} \\\\ (E + \\Omega) p_{00} - (E + \\Omega) p_{11} & (E + \\Omega) p_{01} - \\Omega p_{10}. \n",
        "\\end{bmatrix} -\\frac{\\beta}{2} \n",
        "\\begin{bmatrix} \n",
        "-2 p_{11} & p_{01} \\\\ p_{10} & 2 p_{11} \n",
        "\\end{bmatrix}, \\\\\n",
        "\\dot{\\rho}_{sys}(t) = \n",
        "\\begin{bmatrix} \n",
        "\\dot{p}_{00} & \\dot{p}_{01} \\\\ \\dot{p}_{10} & \\dot{p}_{11}\n",
        "\\end{bmatrix}. \\\\\n",
        "$$\n",
        "\n",
        "Set of four ODEs:\n",
        "$$\n",
        "\\dot{p}_{00} = -i \\Omega p_{10} + i (E + \\Omega) p_{01} + \\beta p_{11}, \\\\\n",
        "\\dot{p}_{01} = -i \\Omega p_{11} + i \\Omega p_{00} - \\frac{\\beta}{2} p_{01}, \\\\\n",
        "\\dot{p}_{10} = -i (E + \\Omega) p_{00} + i (E + \\Omega) p_{11} - \\frac{\\beta}{2} p_{01}, \\\\\n",
        "\\dot{p}_{11} = -i (E + \\Omega) p_{01} + i \\Omega p_{10} - \\beta p_{11}.\n",
        "$$"
      ],
      "metadata": {
        "id": "jBJkXDbeNLYG"
      },
      "id": "jBJkXDbeNLYG"
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Using QuTiP "
      ],
      "metadata": {
        "id": "JiZ76mEnb2hC"
      },
      "id": "JiZ76mEnb2hC"
    },
    {
      "cell_type": "code",
      "source": [
        "# Install QuTiP with pip3 \n",
        "!pip3 install qutip"
      ],
      "metadata": {
        "id": "ZcNGOyvzq_lI"
      },
      "id": "ZcNGOyvzq_lI",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Import required packages\n",
        "import qutip\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt"
      ],
      "metadata": {
        "id": "GH5KEmk0q-hr"
      },
      "id": "GH5KEmk0q-hr",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# System parameters\n",
        "excited_state_energy = 1\n",
        "driving_freq = 0.1\n",
        "decay_rate = 0.3\n",
        "\n",
        "# Starting or initial conditions. Probability corresponds to each pure state \n",
        "p_00 = 0\n",
        "p_01 = 0\n",
        "p_10 = 0\n",
        "p_11 = 1\n",
        "\n",
        "# Define system Hamiltonian\n",
        "Hamiltonian_sys = excited_state_energy * qutip.sigmam() + driving_freq * qutip.sigmax()\n",
        "\n",
        "# Define initial density matrix \n",
        "state_0 = qutip.basis(2,0)\n",
        "state_1 = qutip.basis(2,1)\n",
        "init_density_matrix = p_00 * state_0 * state_0.dag() + p_11 * state_1 * state_1.dag() + p_01 * state_0 * state_1.dag() + p_10 * state_1 * state_0.dag()\n",
        "\n",
        "# Define collaspe operator, equals sqrt(rate) * interacting_operator\n",
        "collapse_op = np.sqrt(decay_rate) * qutip.sigmap()\n",
        "\n",
        "# Define operators to compute population of Fock state 0 and Fock state 1\n",
        "population_state_0 = state_0 * state_0.dag() \n",
        "population_state_1 = state_1 * state_1.dag() "
      ],
      "metadata": {
        "id": "_a_MJKI37r8h"
      },
      "id": "_a_MJKI37r8h",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Define discrete time steps\n",
        "times = np.linspace(0.0, 50.0, 1000)"
      ],
      "metadata": {
        "id": "uBAeNp32gAWy"
      },
      "id": "uBAeNp32gAWy",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Compute expectation value or density matrix at each time step\n",
        "results = qutip.mesolve(Hamiltonian_sys, init_density_matrix, times, collapse_op, [population_state_0, population_state_1])"
      ],
      "metadata": {
        "id": "nz0LWtfedqTu"
      },
      "id": "nz0LWtfedqTu",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plt.figure(figsize=(15, 10), dpi=80)\n",
        "plt.plot(times, results.expect[0])\n",
        "plt.plot(times, results.expect[1])\n",
        "\n",
        "plt.ylabel(\"Expectation value\")\n",
        "plt.xlabel(\"Time\")\n",
        "plt.legend((\"Population at ground state\", \"Population at excited state\"))\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "sGlx10X-g65N"
      },
      "id": "sGlx10X-g65N",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Reference"
      ],
      "metadata": {
        "id": "9k8oDehpo-JB"
      },
      "id": "9k8oDehpo-JB"
    },
    {
      "cell_type": "markdown",
      "source": [
        "D. Manzano, \"A short introduction to the Lindblad master equation,\" AIP Advances 10, 025106 (2020), DOI: 10.1063/1.5115323\n",
        "\n",
        "S. Maniscalco, \"Open quantum systems: Opportunities & challenges,\" Kavli Institute for Theoretical Physics, https://www.youtube.com/watch?v=BVuhHn3-PAM\n",
        "\n",
        "T. Osborne, \"Theory of quantum noise and decoherence, Lecture 7\", University of Hannover, https://www.youtube.com/watch?v=LTj5UL89-ro&list=PLDfPUNusx1EotNvZr1mbu3-0QThQYilFD&index=7\n",
        "\n",
        "J. R. Johansson, P. D. Nation, and F. Nori, \"QuTiP 2: A Python framework for the dynamics of open quantum systems\" Comp. Phys. Comm. 184, 1234 (2013), DOI: 10.1016/j.cpc.2012.11.019\n",
        "\n",
        "J. R. Johansson, P. D. Nation, and F. Nori, \"QuTiP: An open-source Python framework for the dynamics of open quantum systems,\" Comp. Phys. Comm. 183, 1760–1772 (2012), DOI: 10.1016/j.cpc.2012.02.021\n"
      ],
      "metadata": {
        "id": "zGJ6TzV_pB8L"
      },
      "id": "zGJ6TzV_pB8L"
    }
  ],
  "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.6"
    },
    "colab": {
      "provenance": [],
      "collapsed_sections": [
        "94786846",
        "O29HWycvQQ7s",
        "B9AZa8BaaV0U",
        "4HIgDxo2b90V",
        "JiZ76mEnb2hC",
        "9k8oDehpo-JB"
      ]
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}