{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Bosonic Mixture with Excitation Transfer" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using Pkg; Pkg.activate()\n", "using KadanoffBaym, FFTW, Interpolations\n", "using LinearAlgebra\n", "\n", "using PyPlot\n", "PyPlot.plt.style.use(\"./paper.mplstyle\")\n", "using LaTeXStrings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Hamiltonian\n", "\n", "\\begin{align}\n", " \\hat{H} &= \\omega_0 \\sum_{i=1}^{L} \\left(\\hat{b}^\\dagger_{i} \\hat{b}^{\\phantom{\\dagger}}_{i} - \\hat{a}^\\dagger_{i} \\hat{a}^{\\phantom{\\dagger}}_{i} \\right) + J \\sum_{\\langle i, j \\rangle} \\hat{a}^\\dagger_{i}\\hat{b}^{\\phantom{\\dagger}}_{i}\\hat{b}^{\\dagger}_{j}\\hat{a}^{\\phantom{\\dagger}}_{j}\n", "\\end{align}\n", "\n", "### Green functions\n", "\n", "\\begin{align}\n", " \\mathcal{A}_{i,\\,j}(t, t') &= -i \\delta_{ij}\\left\\langle{\\mathcal{T}_\\mathcal{C}\\hat{a}_i^{}(t)\\hat{a}_i^\\dagger(t')}\\right\\rangle \\\\\n", " \\mathcal{B}_{i,\\,j}(t, t') &= -i \\delta_{ij}\\left\\langle{\\mathcal{T}_\\mathcal{C}\\hat{b}_i^{}(t)\\hat{b}_i^\\dagger(t')}\\right\\rangle\n", "\\end{align} \n", "\n", "### Self-energies\n", "\n", "\\begin{align}\n", " \\left[{\\boldsymbol{\\Sigma}^\\lessgtr_{\\mathcal{A}}(t, t')}\\right]_{i,\\, i} &= - J^2\\sum_{j \\in {\\mathcal{N}}_i}\\mathcal{B}^\\lessgtr_{i,\\, i}(t, t') \\mathcal{A}^\\lessgtr_{j,\\, j}(t, t')\\mathcal{B}^\\gtrless_{j,\\, j}(t', t), \\\\\n", " \\left[{\\boldsymbol{\\Sigma}^\\lessgtr_{\\mathcal{B}}(t, t')}\\right]_{i,\\, i} &= - J^2\\sum_{j \\in {\\mathcal{N}}_i}\\mathcal{A}_{i,\\, i}^\\lessgtr(t, t') \\mathcal{A}^\\gtrless_{j,\\, j}(t', t)\\mathcal{B}^\\lessgtr_{j,\\, j}(t, t')\n", "\\end{align} " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Container for problem data\n", "struct ProblemData\n", " GL::GreenFunction{ComplexF64, 4, Array{ComplexF64, 4}, SkewHermitian}\n", " GG::GreenFunction{ComplexF64, 4, Array{ComplexF64, 4}, SkewHermitian}\n", " ΣL::GreenFunction{ComplexF64, 4, Array{ComplexF64, 4}, SkewHermitian}\n", " ΣG::GreenFunction{ComplexF64, 4, Array{ComplexF64, 4}, SkewHermitian}\n", " L::Int64\n", " H::Matrix{ComplexF64}\n", " J::Float64\n", " \n", " # Initialize problem\n", " function ProblemData(GL0::Matrix{ComplexF64}, L::Int64, H::Matrix{ComplexF64}, J::Float64)\n", " @assert H == H' \"Non-hermitian Hamiltonian\"\n", "\n", " data = new(\n", " GreenFunction(reshape(GL0, size(GL0)..., 1, 1), SkewHermitian),\n", " GreenFunction(reshape(GL0 - 1.0im * I, size(GL0)..., 1, 1), SkewHermitian),\n", " GreenFunction(zeros(ComplexF64, size(GL0)..., 1, 1), SkewHermitian),\n", " GreenFunction(zeros(ComplexF64, size(GL0)..., 1, 1), SkewHermitian),\n", " L,\n", " H,\n", " J\n", " )\n", "\n", " # Initialize self-energies\n", " self_energies!(data, 1.0zeros(1), 1.0zeros(1), 1.0zeros(1), 1, 1)\n", "\n", " return data\n", " end\n", "end\n", "\n", "# Self-energy\n", "function self_energies!(data, _, _, _, t1, t2)\n", " (; GL, GG, ΣL, ΣG, L, H, J) = data\n", "\n", " # adjust array size\n", " if (n = size(GL, 3)) > size(ΣL, 3)\n", " resize!(ΣL, n)\n", " resize!(ΣG, n)\n", " end\n", " \n", " # index mapping\n", " NNs = [1, 0] # nearest neighbours\n", " N_mu = mu -> NNs[mu % L + 1] + mu - mu % L\n", " idxs = [[(mu + L) % 2L, N_mu(mu), (N_mu(mu) + L) % 2L] .+ 1 for mu in 0:2L-1]\n", " idxs = hcat(idxs...)\n", " \n", " # self-energies using the index mapping\n", " ΣL[t1, t2] = -J^2 * diag(GL[t1, t2])[idxs[1, :]] .* diag(GL[t1, t2])[idxs[2, :]] .* diag(GG[t2, t1])[idxs[3, :]] |> diagm\n", " ΣG[t1, t2] = -J^2 * diag(GG[t1, t2])[idxs[1, :]] .* diag(GG[t1, t2])[idxs[2, :]] .* diag(GL[t2, t1])[idxs[3, :]] |> diagm\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# integration boundaries\n", "tspan = (0.0, 5.0)\n", "\n", "# problem data\n", "data = begin\n", "\n", " # parameters\n", " L = 2\n", " h = 5.0\n", " H = ComplexF64[-h 0 0 0; 0 -h 0 0; 0 0 h 0; 0 0 0 h];\n", " J = 1.\n", "\n", " # initial condition\n", " delta_n = 0.1\n", " GL0 = -1.0im .* 2e-1 .* [2 0 0 0; 0 1 0 0; 0 0 2 0; 0 0 0 0.5]\n", "\n", " ProblemData(GL0, L, H, J)\n", "end;" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# right-hand side for the \"vertical\" evolution\n", "function fv!(out, data, ts, h1, h2, t1, t2)\n", " (; GL, GG, ΣL, ΣG, L, H) = data\n", " \n", " # real-time collision integral\n", " ∫dt1(A,B,C) = sum(h1[s] * ((A[t1, s] - B[t1, s]) * C[s, t2]) for s in 1:t1)\n", " ∫dt2(A,B,C) = sum(h2[s] * (A[t1, s] * (B[s, t2] - C[s, t2])) for s in 1:t2)\n", "\n", " out[1] = -1.0im * (H * GL[t1, t2] + ∫dt1(ΣG, ΣL, GL) + ∫dt2(ΣL, GL, GG))\n", " out[2] = -1.0im * (H * GG[t1, t2] + ∫dt1(ΣG, ΣL, GG) + ∫dt2(ΣG, GL, GG))\n", "end\n", "\n", "# right-hand side for the \"diagonal\" evolution\n", "function fd!(out, data, ts, h1, h2, t1, t2)\n", " fv!(out, data, ts, h1, h2, t1, t2)\n", " out .-= adjoint.(out)\n", "end\n", "\n", "# call the solver\n", "sol = kbsolve!(\n", " (out, x...) -> fv!(out, data, x...), \n", " (out, x...) -> fd!(out, data, x...), \n", " [data.GL, data.GG], \n", " tspan; \n", " callback = (x...) -> self_energies!(data, x...), \n", " atol=1e-10, \n", " rtol=1e-8,\n", " stop = x -> (println(\"t: $(x[end])\"); flush(stdout); false)\n", ");" ] }, { "cell_type": "markdown", "metadata": { "heading_collapsed": true }, "source": [ "## QuTiP benchmark" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hidden": true }, "outputs": [], "source": [ "using PyCall\n", "qt = pyimport(\"qutip\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hidden": true, "scrolled": true }, "outputs": [], "source": [ "# time parameters\n", "times = range(first(sol.t), stop=last(sol.t), length=128+1)\n", "n = length(times) - 1" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "### Hamiltonian" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hidden": true }, "outputs": [], "source": [ "n_max = 2; # Fock-space truncation\n", "\n", "# initial state\n", "psi0 = qt.tensor([(sqrt(1 - 1.0im * data.GL[1, 1][k, k]) * qt.basis(n_max + 1, 0) \n", " + sqrt(1.0im * data.GL[1, 1][k, k]) * qt.basis(n_max + 1, 1)).unit() for k in 1:2*L]);\n", "\n", "# operators\n", "ids = [qt.qeye(n_max + 1), qt.qeye(n_max + 1), qt.qeye(n_max + 1), qt.qeye(n_max + 1)]\n", "ops = [deepcopy(ids) for _ in 1:4]\n", "\n", "for (i, op) in enumerate(ops)\n", " op[i] = qt.destroy(n_max + 1)\n", "end\n", "\n", "ops = qt.tensor.(ops)\n", "b1, b2, a1, a2 = ops;" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hidden": true }, "outputs": [], "source": [ "# make diagonal density matrix (i.e. dropping coherences)\n", "rho0 = psi0 * psi0.dag();\n", "for j in 1:rho0.shape[2] , i in 1:rho0.shape[1]\n", " i != j ? rho0.data[i, j] = 0.0 : continue\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hidden": true }, "outputs": [], "source": [ "# Hamiltonian (hard-coded for two sites)\n", "H = -h * b1.dag() * b1\n", "H += -h * b2.dag() * b2\n", "H += +h * a1.dag() * a1\n", "H += +h * a2.dag() * a2\n", "H += J * (a1.dag() * b1 * b2.dag() * a2 + b1.dag() * a1 * a2.dag() * b2)\n", "\n", "# observables\n", "obs = [b1.dag() * b1, b2.dag() * b2, a1.dag() * a1, a2.dag() * a2];" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "### Simulation" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hidden": true }, "outputs": [], "source": [ "# quickly solve once for observables\n", "me = qt.mesolve(H, rho0, times, [], obs)\n", "\n", "# solve for the time-dependent density matrix\n", "t_sols = qt.mesolve(H, rho0, times); # t_sols.states returns density matrices" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "#### Two times" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hidden": true }, "outputs": [], "source": [ "# compute two-time average -i\n", "function two_time_average(X, Y, rho_t1)\n", " X_t2_Y_t1 = zeros(ComplexF64, n + 1, n + 1)\n", " for k in 1:(n + 1)\n", " Y_rho_t1 = qt.mesolve(H, Y * rho_t1.states[k], times).states\n", " for l in 1:(n + 1)\n", " X_t2_Y_t1[k, l] = -1.0im * (X * Y_rho_t1[l]).tr() \n", " end\n", " end\n", " \n", " # transform to regular two-time square used by the KB solver\n", " rotated_X_t2_Y_t1 = zeros(ComplexF64, n + 1, 2*(n + 1) - 1)\n", " for (k, x) in enumerate([X_t2_Y_t1[k, :] for k in 1:(n + 1)])\n", " for (l, y) in enumerate(x)\n", " ind = k + l - 1\n", " rotated_X_t2_Y_t1[k, ind] = y \n", " end\n", " end \n", " \n", " return rotated_X_t2_Y_t1[:, 1:n+1]\n", "end;" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hidden": true }, "outputs": [], "source": [ "create_destroy = [zeros(ComplexF64, n + 1, n + 1) for _ in 1:1]#2L]\n", "destroy_create = [zeros(ComplexF64, n + 1, n + 1) for _ in 1:1]#2L];" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hidden": true }, "outputs": [], "source": [ "destroyers = [b1, b2, a1, a2]\n", "for k in 1:1#2L\n", " create_destroy[k] = two_time_average(destroyers[k].dag(), destroyers[k], t_sols)\n", " destroy_create[k] = two_time_average(destroyers[k], destroyers[k].dag(), t_sols)\n", " \n", " # flip real part when creation operator is to the right (equivalent to flipping the sign of tau)\n", " destroy_create[k] = -conj(destroy_create[k])\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hidden": true }, "outputs": [], "source": [ "b1_dag_b1 = create_destroy[1]\n", "b1_b1_dag = destroy_create[1];" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hidden": true }, "outputs": [], "source": [ "# compute missing two-time triangle from symmetry\n", "full_square = A -> (A - adjoint(A) - (A |> diag |> diagm));" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hidden": true }, "outputs": [], "source": [ "figure(figsize=(6, 2))\n", "subplot(121)\n", "imshow(real(full_square(b1_dag_b1)), cmap=\"plasma\")\n", "\n", "subplot(122)\n", "imshow(real(full_square(b1_b1_dag)), cmap=\"plasma\")\n", "\n", "tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# final time\n", "T = sol.t[end]\n", "\n", "# scale parameters\n", "t_scale = (J == 0. ? 1. : J)\n", "ω_scale = (J == 0. ? 1. : 1/J);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Equal time" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xpad = 8\n", "ypad = 5\n", "\n", "figure(figsize=(7, 3))\n", "\n", "ax = subplot(121)\n", "idx = 1\n", "ax.plot(sol.t, data.GL.data[idx, idx ,:,:] |> ((-) ∘ imag ∘ diag), ls=\"-\", label=\"\\$i=\"*string((idx - 1) % 2 + 1)*\"\\$\", c=\"C0\")\n", "ax.plot(times, me.expect[idx], \"--\", c=\"C0\", lw=3, alpha=0.5)\n", "idx = 2\n", "ax.plot(sol.t, data.GL.data[idx, idx ,:,:] |> ((-) ∘ imag ∘ diag), ls=\":\", label=\"\\$i=\"*string((idx - 1) % 2 + 1)*\"\\$\", c=\"C1\")\n", "ax.plot(times, me.expect[idx], \"-.\", c=\"C1\", lw=3, alpha=0.5)\n", "ax.set_xlabel(\"\\$Jt\\$\") \n", "ylabel(\"\\$-\\\\mathrm{Im}\\\\; \\\\mathcal{A}^<_{i,\\\\, i}(t, t)\\$\", labelpad=10)\n", "ax.set_xticks([0, 2.5, 5])\n", "ax.set_xlim(0, sol.t[end]) \n", "ax.set_ylim(0, 0.5)\n", "ax.legend(loc=\"best\", handlelength=1.4, frameon=false, borderpad=0, labelspacing=0.25)\n", "\n", "ax = subplot(122)\n", "idx = 3\n", "ax.plot(sol.t, data.GL.data[idx, idx ,:,:] |> ((-) ∘ imag ∘ diag), ls=\"-\", label=\"\\$i=\"*string((idx - 1) % 2 + 1)*\"\\$\", c=\"C2\")\n", "ax.plot(times, me.expect[idx], \"--\", c=\"C2\", lw=3, alpha=0.5)\n", "idx = 4\n", "ax.plot(sol.t, data.GL.data[idx, idx ,:,:] |> ((-) ∘ imag ∘ diag), ls=\":\", label=\"\\$i=\"*string((idx - 1) % 2 + 1)*\"\\$\", c=\"C3\")\n", "ax.plot(times, me.expect[idx], \"-.\", c=\"C3\", lw=3, alpha=0.5)\n", "ax.set_xlabel(\"\\$Jt\\$\") \n", "ylabel(\"\\$-\\\\mathrm{Im}\\\\; \\\\mathcal{B}^<_{i,\\\\, i}(t, t)\\$\", labelpad=16)\n", "ax.yaxis.set_label_position(\"right\")\n", "ax.set_xticks([0, 2.5, 5])\n", "ax.set_yticklabels([])\n", "ax.set_xlim(0, sol.t[end]) \n", "ax.set_ylim(0, 0.5)\n", "ax.legend(loc=\"best\", handlelength=1.4, frameon=false, borderpad=0, labelspacing=0.25)\n", "tight_layout(pad=0.25, w_pad=1, h_pad=0)\n", "\n", "# savefig(\"interacting_bosons_example_1.pdf\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Relative time" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# quantum number to look at\n", "idx = 1;" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "ρ_11_kb = interpolate((sol.t, sol.t), view(data.GL.data .- data.GG.data, idx, idx, :, : ), Gridded(Linear()));\n", "ρ_11_qt = interpolate((times, times), view(full_square(create_destroy[idx] .- destroy_create[idx]), :, : ), Gridded(Linear()));" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "new_times = range(first(sol.t), stop=last(sol.t), length=2048);\n", "\n", "ρ_11_kb_wigner, (taus, ts) = wigner_transform([ρ_11_kb(t1, t2) for t1 in new_times, t2 in new_times]; ts=new_times, fourier=false);\n", "ρ_11_qt_wigner, (taus_qt, ts_qt) = wigner_transform([ρ_11_qt(t1, t2) for t1 in new_times, t2 in new_times]; ts=new_times, fourier=false);" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xpad = 8\n", "ypad = 5\n", "\n", "center = floor(length(new_times) / 2) |> Int\n", "\n", "figure(figsize=(7, 3))\n", "\n", "ax = subplot(121)\n", "plot(t_scale * taus, ρ_11_kb_wigner[:, center] |> imag, ls=\"-\", c=\"C0\", lw=1.5)\n", "plot(t_scale * taus_qt, ρ_11_qt_wigner[:, center] |> imag, ls=\"--\", c=\"C0\", lw=2.5, alpha=0.5)\n", "\n", "ax.xaxis.set_tick_params(pad=xpad)\n", "ax.yaxis.set_tick_params(pad=ypad)\n", "ax.set_xlabel(L\"J \\tau\")\n", "ax.set_ylabel(L\"\\mathrm{Re}\\,A_{\\mathcal{A}_{1,\\, 1}}(T, \\tau)_W\")\n", "ax.set_xlim(-t_scale * T, t_scale * T)\n", "ax.set_ylim(-1.0, 1.0)\n", "\n", "ax = subplot(122)\n", "plot(t_scale * taus, -ρ_11_kb_wigner[:, center] |> real, ls=\"-\", c=\"C0\", lw=1.5)\n", "plot(t_scale * taus_qt, -ρ_11_qt_wigner[:, center] |> real, ls=\"--\", c=\"C0\", lw=2.5, alpha=0.5)\n", "ax.set_xlabel(L\"J \\tau\")\n", "ax.set_xlim(-t_scale * T, t_scale * T)\n", "ax.set_ylim(-1.0, 1.0)\n", "ax.set_yticklabels([])\n", "ax.xaxis.set_tick_params(pad=xpad)\n", "ax.yaxis.set_tick_params(pad=ypad)\n", "ax.set_ylabel(L\"\\mathrm{Im}\\,A_{\\mathcal{A}_{1,\\, 1}}(T, \\tau)_W\", labelpad=16)\n", "ax.yaxis.set_label_position(\"right\")\n", "\n", "tight_layout(pad=0.1, w_pad=1, h_pad=0)\n", "\n", "# savefig(\"interacting_bosons_example_2.pdf\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.7.0", "language": "julia", "name": "julia-1.7" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.7.0" } }, "nbformat": 4, "nbformat_minor": 4 }