2.Stochastic quantum dynamics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Info: Nathan Shammah, RIKEN, nathan.shammah@gmail.com\n",
    "\n",
    "# Open Quantum Dynamics with QuTiP\n",
    " \n",
    " \n",
    " We use QuTiP's solvers to study the open dynamics of a quantum system evolving in time. \n",
    " \n",
    " \n",
    " - #### Lindblad Master Equation\n",
    " \n",
    "  - #### Lindblad Master Equation: Time evolution\n",
    " \n",
    "  - #### Lindblad Master Equation: steady-state solution\n",
    " \n",
    "- #### Stochastic Dynamics: Quantum Trajectories with Monte Carlo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from qutip import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the operators and Hamiltonian\n",
    "\n",
    "We consider the spin-boson system, which is a prototypical model of light-matter interaction in cavity quantum electrodynamics (cQED): a single two-level system coupled to a single mode of the photonic field. Its Hamiltonian is described by\n",
    "\n",
    "\\begin{eqnarray}\n",
    "H &=& \\omega_c a^\\dagger a + \\frac{\\omega_0}{2}\\sigma_z\n",
    "+\\frac{g}{2}\\sigma_x\\left(a+a^\\dagger\\right)+\\frac{\\omega_x}{2}\\sigma_x,\n",
    "\\end{eqnarray}\n",
    "\n",
    "where we added also the possibility of the system of a classical drive onto the system at a frequency $\\omega_x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# spins\n",
    "sx_reduced = sigmax()\n",
    "sy_reduced = sigmay()\n",
    "sz_reduced = sigmaz()\n",
    "sp_reduced = sigmap()\n",
    "sm_reduced = sigmam()\n",
    "\n",
    "# photons\n",
    "nph = 4\n",
    "a_reduced = destroy(nph)\n",
    "\n",
    "# tensor space\n",
    "sz = tensor(sz_reduced,qeye(nph))\n",
    "sx = tensor(sx_reduced,qeye(nph))\n",
    "sm = tensor(sm_reduced,qeye(nph))\n",
    "sp = sm.dag()\n",
    "a = tensor(qeye(2), a_reduced)\n",
    "\n",
    "# hamiltonians\n",
    "wc = 1\n",
    "w0 = 0.5*wc\n",
    "g = 0.1*w0\n",
    "wx = 0.2*w0\n",
    "Hcav = wc*a.dag()*a\n",
    "Hspin = w0*sz #+ wx*sx\n",
    "Hint =  g*sx*(a+a.dag())\n",
    "HintCheck =  g*tensor(sx_reduced,a_reduced+a_reduced.dag())\n",
    "H = Hcav + Hspin + Hint\n",
    "np.testing.assert_(Hint == HintCheck)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the initial state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial state of the system is given by \n",
    "\\begin{eqnarray}\n",
    "\\rho&=&\\rho_\\text{spin}\\otimes\\rho_\\text{phot}\n",
    "\\end{eqnarray}\n",
    "and in the case of a initially pure state, \n",
    "\\begin{eqnarray}\n",
    "\\rho&=&|\\psi\\rangle_\\text{spin}\\otimes|\\psi\\rangle_\\text{phot}\n",
    "\\end{eqnarray}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initial state\n",
    "psi0_spin = basis(2,0)\n",
    "psi0_phot = basis(nph,nph-int(nph/2))\n",
    "psi0 = tensor(psi0_spin,psi0_phot)\n",
    "rho0 = ket2dm(psi0)\n",
    "\n",
    "# times at which to calculate the variables\n",
    "tlist = np.linspace(0,50,2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lindblad master equation: $\\texttt{mesolve}$\n",
    "\n",
    "We now consider the time evolution of the open quantum system, in which $\\rho$ is dissipatively coupled to a spin and photonic bath,\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\frac{d}{dt}\\rho =-i\\lbrack H,\\rho\\rbrack+\\gamma\\mathcal{D}_{[\\sigma_-]}\\rho+\\kappa\\mathcal{D}_{[a]}\\rho.\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kappa = 0.3\n",
    "gamma = 0.3\n",
    "\n",
    "my_options = Options(average_states = True, store_states = True)\n",
    "results = mesolve(H, psi0, tlist, \n",
    "                  c_ops=[np.sqrt(kappa)*a,np.sqrt(gamma)*sz], \n",
    "                  e_ops=[a.dag()*a,sz], \n",
    "                  options=my_options, progress_bar=True)\n",
    "\n",
    "# store time evoluted variables\n",
    "nph_t = results.expect[0]\n",
    "sz_t = results.expect[1]\n",
    "rho_t = results.states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fs=20\n",
    "plt.figure(figsize=(14,4))\n",
    "plt.subplot(121)\n",
    "plt.plot(tlist, sz_t)\n",
    "plt.xlabel(r\"$t$\",fontsize=fs)\n",
    "plt.ylabel(r\"$\\langle \\sigma_z \\rangle$\",fontsize=fs)\n",
    "plt.subplot(122)\n",
    "plt.plot(tlist, nph_t/nph/0.5)\n",
    "plt.xlabel(r\"$t$\",fontsize=fs)\n",
    "plt.ylabel(r\"$\\langle a^\\dagger a \\rangle$\",fontsize=fs)\n",
    "plt.suptitle(\"Time evolution for spin and photonic excitations of a driven cavity with single two level system, with dissipation\")\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Steady state solver: $\\texttt{steadystate}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rhoss = steadystate(H,[np.sqrt(kappa)*a,np.sqrt(gamma)*sz])\n",
    "nph_ss= expect(a.dag()*a,rhoss)\n",
    "sz_ss= expect(sz,rhoss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(abs(rhoss.full()))\n",
    "plt.show()\n",
    "plt.close()\n",
    "plt.figure()\n",
    "plt.subplot(121)\n",
    "plt.imshow(abs(ptrace(rhoss,0).full()))\n",
    "plt.xlabel(\"atom\",fontsize=fs)\n",
    "plt.subplot(122)\n",
    "plt.imshow(abs(ptrace(rhoss,1).full()))\n",
    "plt.xlabel(\"photon\",fontsize=fs)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Liouvillian structure\n",
    "\n",
    "We can derive thermodynamical properties of the out-of-equilibrium system by studying the spectrum of the Liouvillian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = liouvillian(H,[np.sqrt(kappa)*a,np.sqrt(gamma)*sz])\n",
    "\n",
    "# represent the Liouvillian\n",
    "plt.figure(figsize=(14,6))\n",
    "plt.subplot(121)\n",
    "plt.imshow(np.real(L.full()))\n",
    "plt.title(\"Real values of the Liouvillian\",fontsize=fs)\n",
    "plt.subplot(122)\n",
    "plt.imshow(np.imag(L.full()))\n",
    "plt.title(\"Imaginary values of the Liouvillian\",fontsize=fs)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Liouvillian Spectrum, $\\texttt{eigenstates}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = liouvillian(H,[np.sqrt(kappa)*a,np.sqrt(gamma)*sz])\n",
    "\n",
    "# Plot the Liouvillian spectrum in the complex plane\n",
    "eigenvalues_L, eigenvectors_L = L.eigenstates()\n",
    "\n",
    "\n",
    "real_eigenvalues = np.real(eigenvalues_L)\n",
    "imag_eigenvalues = np.imag(eigenvalues_L)\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.plot(real_eigenvalues, imag_eigenvalues,\"o\")\n",
    "plt.title(\"Real values of the Liouvillian\",fontsize=fs)\n",
    "plt.ylabel(\"Im($\\lambda_i$)\",fontsize=fs)\n",
    "plt.xlabel(\"Re($\\lambda_i$)\",fontsize=fs)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us decipher the spectrum of the Liouvillian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,10))\n",
    "plt.plot(real_eigenvalues, imag_eigenvalues,\"o\")\n",
    "plt.title(\"Real values of the Liouvillian\")\n",
    "plt.axvline(x=0,color=\"red\")\n",
    "plt.axhline(y=0,color=\"black\",linestyle=\"dashed\")\n",
    "plt.ylabel(\"Im($\\lambda_i$)\",fontsize=fs)\n",
    "plt.xlabel(\"Re($\\lambda_i$)\",fontsize=fs)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We find symmetries in the Liouvillian spectrum with respect to the real and imaginary axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6,6))\n",
    "plt.plot(real_eigenvalues[-1],imag_eigenvalues[-1],\"o\",markersize=20,label=\"$\\lambda_0$\")\n",
    "plt.plot(real_eigenvalues[-2],imag_eigenvalues[-2],\"*\",markersize=20,label=\"$\\lambda_1$\")\n",
    "plt.plot(real_eigenvalues[-3],imag_eigenvalues[-3],\"s\",markersize=20,label=\"$\\lambda_2$\")\n",
    "plt.xlim([-2,2])\n",
    "plt.ylim([-2,2])\n",
    "plt.legend(fontsize=20)\n",
    "plt.ylabel(\"Im($\\lambda_i$)\",fontsize=fs)\n",
    "plt.xlabel(\"Re($\\lambda_i$)\",fontsize=fs)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Photon counting statistics: $\\texttt{mcsolve}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set dynamics options\n",
    "my_options = Options(average_states = False, store_states = True)\n",
    "# solve dynamics\n",
    "results_mc = mcsolve(H, psi0, tlist, \n",
    "                  c_ops=[np.sqrt(kappa)*a,np.sqrt(gamma)*sz], \n",
    "                  e_ops=[a.dag()*a,sz], \n",
    "                  options=my_options, progress_bar=True)\n",
    "\n",
    "# store time evoluted variables\n",
    "nph_mc_t = results_mc.expect[0]\n",
    "sz_mc_t = results_mc.expect[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default options for $\\texttt{mcsolve}$ imply solving the dynamics for 500 trajectories. \n",
    "\n",
    "This option can be controlled by setting a different number in $\\texttt{Options(ntraj=500})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_mc_t = results_mc.states\n",
    "len(rho_mc_t)\n",
    "#help(expect)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sz_stoch_t = []\n",
    "nph_stoch_t = []\n",
    "for i in range(len(rho_mc_t)):\n",
    "    sz_stoch_t.append(expect(sz,rho_mc_t[i]))\n",
    "    nph_stoch_t.append(expect(a.dag()*a,rho_mc_t[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,4))\n",
    "plt.subplot(121)\n",
    "plt.plot(tlist, sz_t, label=\"me solve\")\n",
    "plt.plot(tlist, sz_mc_t, label=\"mc solve (500 trajectories)\")\n",
    "plt.xlabel(r\"$t$\",fontsize=fs)\n",
    "plt.ylabel(r\"$\\langle \\sigma_z \\rangle$\",fontsize=fs)\n",
    "plt.legend(fontsize=15)\n",
    "plt.subplot(122)\n",
    "plt.plot(tlist, nph_t/nph/0.5, label=\"me solve\")\n",
    "plt.plot(tlist, nph_mc_t/nph/0.5, label=\"mc solve (500 trajectories)\")\n",
    "plt.xlabel(r\"$t$\",fontsize=fs)\n",
    "plt.ylabel(r\"$\\langle a^\\dagger a \\rangle$\",fontsize=fs)\n",
    "plt.legend(fontsize=15)\n",
    "plt.show()\n",
    "plt.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lindblad master equation as the limit of the stochastic evolution: $\\texttt{mesolve}$ and $\\texttt{mcsolve}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 50 trajectories\n",
    "my_options50 = Options(ntraj=50)\n",
    "results_mc50 = mcsolve(H, psi0, tlist, c_ops=[np.sqrt(kappa)*a,np.sqrt(gamma)*sz], \n",
    "                  e_ops=[a.dag()*a,sz], options=my_options50, progress_bar=True)\n",
    "nph_mc_t50 = results_mc50.expect[0]\n",
    "sz_mc_t50 = results_mc50.expect[1]\n",
    "# 100 trajectories\n",
    "my_options100 = Options(ntraj=100)\n",
    "results_mc100 = mcsolve(H, psi0, tlist, c_ops=[np.sqrt(kappa)*a,np.sqrt(gamma)*sz], \n",
    "                  e_ops=[a.dag()*a,sz], options=my_options100, progress_bar=True)\n",
    "nph_mc_t100 = results_mc100.expect[0]\n",
    "sz_mc_t100 = results_mc100.expect[1]\n",
    "# 200 trajectories\n",
    "my_options200 = Options(ntraj=200)\n",
    "results_mc200 = mcsolve(H, psi0, tlist, c_ops=[np.sqrt(kappa)*a,np.sqrt(gamma)*sz], \n",
    "                  e_ops=[a.dag()*a,sz], options=my_options200, progress_bar=True)\n",
    "nph_mc_t200 = results_mc200.expect[0]\n",
    "sz_mc_t200 = results_mc200.expect[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,6))\n",
    "# spin excitation\n",
    "plt.subplot(121)\n",
    "plt.plot(tlist, sz_mc_t50, label=\"mc solve (50 trajectories)\")\n",
    "plt.plot(tlist, sz_mc_t100, label=\"mc solve (100 trajectories)\")\n",
    "plt.plot(tlist, sz_mc_t200, label=\"mc solve (200 trajectories)\")\n",
    "plt.plot(tlist, sz_mc_t, label=\"mc solve (500 trajectories)\")\n",
    "plt.plot(tlist, sz_t, label=\"me solve\")\n",
    "plt.xlabel(r\"$t$\",fontsize=fs)\n",
    "plt.ylabel(r\"$\\langle \\sigma_z \\rangle$\",fontsize=fs)\n",
    "plt.legend(fontsize=15)\n",
    "# photonic excitation\n",
    "plt.subplot(122)\n",
    "plt.plot(tlist, nph_mc_t50/nph/0.5, label=\"mc solve (50 trajectories)\")\n",
    "plt.plot(tlist, nph_mc_t100/nph/0.5, label=\"mc solve (100 trajectories)\")\n",
    "plt.plot(tlist, nph_mc_t200/nph/0.5, label=\"mc solve (200 trajectories)\")\n",
    "plt.plot(tlist, nph_mc_t/nph/0.5, label=\"mc solve (500 trajectories)\")\n",
    "plt.plot(tlist, nph_t/nph/0.5, label=\"me solve\")\n",
    "plt.xlabel(r\"$t$\",fontsize=fs)\n",
    "plt.ylabel(r\"$\\langle a^\\dagger a \\rangle$\",fontsize=fs)\n",
    "plt.legend(fontsize=15)\n",
    "plt.show()\n",
    "plt.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe that progressively, as the number of trajectories is increased, the value given by the average over the trajectories gets closer to the master equation value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,4))\n",
    "plt.subplot(121)\n",
    "plt.plot(tlist, sz_t)\n",
    "plt.plot(tlist, sz_mc_t)\n",
    "plt.xlabel(r\"$t$\",fontsize=fs)\n",
    "plt.ylabel(r\"$\\langle s_z \\rangle$\",fontsize=fs)\n",
    "plt.subplot(122)\n",
    "plt.plot(tlist, nph_t/nph/0.5)\n",
    "plt.plot(tlist, nph_mc_t/nph/0.5)\n",
    "\n",
    "plt.xlabel(r\"$t$\",fontsize=fs)\n",
    "plt.ylabel(r\"$\\langle n_{ph} \\rangle$\",fontsize=fs)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting single quantum trajectories\n",
    "\n",
    "By setting:\n",
    "my_options = Options(average_states = False, store_states = True)\n",
    "\n",
    "we can inspect single quantum trajectories and their behaviour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,5))\n",
    "plt.subplot(121)\n",
    "for i in range(len(sz_stoch_t)):\n",
    "    plt.plot(tlist, sz_stoch_t[i],alpha=0.1)\n",
    "\n",
    "plt.plot(tlist, sz_t,\"b\",linewidth=3)\n",
    "plt.plot(tlist, sz_mc_t,\"y\",linewidth=2)\n",
    "plt.plot(tlist, sz_stoch_t[4]/nph/0.5,\"g\",linewidth=5)\n",
    "plt.plot(tlist, sz_stoch_t[6]/nph/0.5,\"r\",linewidth=5)\n",
    "plt.subplot(122)\n",
    "for i in range(len(sz_stoch_t)):\n",
    "    plt.plot(tlist, nph_stoch_t[i]/nph/0.5,alpha=0.1)\n",
    "plt.plot(tlist, nph_t/nph/0.5,\"b\",linewidth=3)\n",
    "plt.plot(tlist, nph_mc_t/nph/0.5,\"y\",linewidth=2)\n",
    "plt.plot(tlist, nph_stoch_t[4]/nph/0.5,\"g\",linewidth=5)\n",
    "plt.plot(tlist, nph_stoch_t[6]/nph/0.5,\"r\",linewidth=5)\n",
    "plt.xlabel(r\"$t$\",fontsize=fs)\n",
    "plt.ylabel(r\"$\\langle X \\rangle$\",fontsize=fs)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We obseverve that the single trajectories can provide information about the phase of the system. In presence of a phase transition the quantum trajectory jumps between possible choices of the degenerate steady state. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qutip.about()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}