{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ad5f297b",
   "metadata": {},
   "source": [
    "# Bose-Einstein Condensates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9dee761",
   "metadata": {},
   "outputs": [],
   "source": [
    "using KadanoffBaym\n",
    "using LinearAlgebra\n",
    "\n",
    "using PyPlot\n",
    "# PyPlot.plt.style.use(\"./paper.mplstyle\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ee0896c",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac99780f",
   "metadata": {},
   "source": [
    "In this notebook we will use `KadanoffBaym.jl` to study _dephasing_ in Bose-Einstein condensates (see Chp. 3 [here](https://bonndoc.ulb.uni-bonn.de/xmlui/handle/20.500.11811/8961)). To do this, we will need to solve _one-time_ differential equations for the condensate amplitude $\\varphi(t)$ and the equal-time Keldysh Green function $G^K(t, t)$.\n",
    "\n",
    "The Lindblad master equation describing this systems reads\n",
    "\n",
    "\\begin{align}\n",
    "    \\begin{split}\n",
    "    \t\\dot{\\hat{\\rho}} &= -i\\omega_0[a^\\dagger a, \\hat{\\rho}] +\\frac{\\lambda}{2}\\left\\{ 2a\\hat{\\rho} a^{\\dagger} - \\left( a^{\\dagger}a\\hat{\\rho} + \\hat{\\rho} a^{\\dagger}a \\right)\\right\\} + \\frac{\\gamma}{2}\\left\\{ 2a^{\\dagger}\\hat{\\rho} a - \\left( aa^{\\dagger}\\hat{\\rho} + \\hat{\\rho} aa^{\\dagger} \\right)\\right\\} \\\\\n",
    "        &+ D\\left\\{ 2a^{\\dagger}a\\hat{\\rho} a^{\\dagger}a - \\left( a^{\\dagger}aa^{\\dagger}a\\hat{\\rho} + \\hat{\\rho} a^{\\dagger}aa^{\\dagger}a \\right)\\right\\},    \n",
    "    \\end{split}\n",
    "\\end{align}\n",
    "\n",
    "where $\\lambda > 0$ is the loss parameter, $\\gamma > 0$ represents the corresponding gain, and $D > 0$ is the constant that introduces dephasing. The derivation for the equations of motion for $\\varphi(t)$ and $G^K(t, t)$ is again given [here](https://bonndoc.ulb.uni-bonn.de/xmlui/handle/20.500.11811/8961) and leads to\n",
    "\n",
    "\\begin{align}\n",
    "    \\begin{split}\n",
    "    \\dot{\\varphi}(t) &=  -i\\omega_0\\varphi(t) -\\frac{1}{2}{(\\lambda - \\gamma + {2} D)}\\varphi(t), \\\\\n",
    "    \\dot{G}^K(t, t) &= -{(\\lambda - \\gamma)}G^K(t, t) - i{\\left(\\lambda + \\gamma + {2} D |\\varphi(t)|^2\\right)}.\n",
    "    \\end{split}\n",
    "\\end{align}\n",
    "\n",
    "To make these expressions more transparent, we set $\\varphi(t) = \\sqrt{2N(t)}\\mathrm{e}^{i \\theta(t)}$ and $G^K(t, t) = -i{(2\\delta N(t) + 1)}$, where $N$ and $\\delta N$ are the condensate and non-condensate occupation, respectively. For these quantities, we obtain\n",
    "\n",
    "\\begin{align}\n",
    "    \\begin{split}\n",
    "    \\dot{N} &=  {(\\gamma - \\lambda -{2} D)}N, \\\\\n",
    "    \\delta \\dot{N} &=  \\gamma{(\\delta N + 1)} - \\lambda\\delta N + {2}DN.\n",
    "    \\end{split}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c74db597",
   "metadata": {},
   "source": [
    "## Defining the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b16fea3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "ω₀ = 1.0\n",
    "λ = 0.0\n",
    "γ = 0.0\n",
    "D = 1.0 \n",
    "\n",
    "# initial occupations\n",
    "N = 1.0\n",
    "δN = 0.0\n",
    "\n",
    "# One-time function for the condensate\n",
    "φ = GreenFunction(zeros(ComplexF64, 1), OnePoint)\n",
    "\n",
    "# Allocate the initial Green functions (time arguments at the end)\n",
    "GK = GreenFunction(zeros(ComplexF64, 1, 1), SkewHermitian)\n",
    "\n",
    "# Initial conditions\n",
    "GK[1, 1] = -im * (2δN + 1)\n",
    "φ[1] = sqrt(2N);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4410fa04",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we leave the vertical equation empty since we can solve for GK in equal-time only\n",
    "function fv!(out, ts, h1, h2, t, t′)\n",
    "    out[1] = zero(out[1])\n",
    "end\n",
    "\n",
    "# diagonal equation for GK\n",
    "function fd!(out, ts, h1, h2, t, _)\n",
    "    out[1] = -(λ - γ) * GK[t, t] - im * (λ + γ + 2D * abs2(φ[t]))\n",
    "end\n",
    "\n",
    "# one-time equation for condensate amplitude\n",
    "function f1!(out, ts, h1, t)\n",
    "    out[1] = -im * ω₀ * φ[t] - (1/2) * (λ - γ + 2D) * φ[t]\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc1e11fb",
   "metadata": {},
   "source": [
    "## Solving the example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32403551",
   "metadata": {},
   "outputs": [],
   "source": [
    "# call the solver\n",
    "sol = kbsolve!(fv!, fd!, [GK,], (0.0, 1.0); atol=1e-6, rtol=1e-4, v0 = [φ,], f1! =f1!);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80efc51a",
   "metadata": {},
   "source": [
    "## Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c876e1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "let\n",
    "    fig = figure(figsize=(4, 3))\n",
    "\n",
    "    ax = subplot(111)\n",
    "    plot(ω₀ .* sol.t, abs2.(φ[:]) / 2 |> real, ls=\"-\", c=\"C0\", label=\"\\$\\\\varphi(t)\\$\", lw=1.5)\n",
    "    plot(ω₀ .* sol.t, [(im * GK[k, k] - 1) / 2 for k in 1:size(sol.t)[1]] |> real, ls=\"-\", c=\"C1\", label=\"\\$G^K(t, t)\\$\", lw=1.5)\n",
    "    ax.set_xlabel(\"\\$ \\\\omega_0 t\\$\")\n",
    "    ax.set_xlim(0, sol.t[end])\n",
    "    ax.set_ylim(0.0, 1.0)\n",
    "    ax.legend(frameon=false)\n",
    "    tight_layout()\n",
    "    fig\n",
    "end;"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.0",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}