diff --git a/notebooks/numerical-derivation-of-master-equation.ipynb b/notebooks/numerical-derivation-of-master-equation.ipynb new file mode 100644 index 0000000..479d585 --- /dev/null +++ b/notebooks/numerical-derivation-of-master-equation.ipynb @@ -0,0 +1,195 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*This notebook can be found on* [github](https://github.com/qojulia/QuantumOptics.jl-examples/tree/master/notebooks/numerical-derivation-of-master-equation.ipynb)\n", + "\n", + "# Numerical derivation of a master equation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Open quantum systems, i.e. systems coupled to an environment, can under certain conditions be described very well by a master equation. Its derivation is a bit technical but can very easily be reproduced numerically. Roughly following the approach in [1], we start with a system consisting of a single 2-level atom coupled to infinitely many modes of the electromagnetic field. The time dynamics of this system are governed by a Schrödinger equation with Hamiltonian\n", + "$$\n", + "H = \\Omega \\sigma^+ \\sigma^-\n", + " + \\sum_k \\omega_k a_k^\\dagger a_k\n", + " + \\sum_k \\kappa(\\omega_k) (\\sigma^+ a_k + \\sigma^- a_k^\\dagger)\n", + "$$\n", + "For simplicity we assume that $\\kappa(\\omega) = \\kappa_0$ for all relevant frequencies and that the selected frequencies are evenly spaced. Further on we change in a rotating frame and obtain\n", + "$$\n", + "H = \\sum_k \\Delta_k a_k^\\dagger a_k\n", + " + \\sum_k \\kappa_0 (\\sigma^+ a_k + \\sigma^- a_k^\\dagger).\n", + "$$\n", + "We will now demonstrate that for appropriate choices of $\\Delta_k/\\kappa_0$ the dynamics of the reduced system, where the electromagnetic modes are traced out, are given by the master equation\n", + "$$\n", + "\\dot \\rho = -i [\\Omega \\sigma^+ \\sigma^-, \\rho]\n", + " + \\gamma (\\sigma^- \\rho \\sigma^+\n", + " - \\frac{1}{2} \\sigma^+ \\sigma^- \\rho\n", + " - \\frac{1}{2} \\rho \\sigma^+ \\sigma^-)\n", + "$$\n", + "The decay rate $\\gamma$ is here connected to the coupling strength $\\kappa_0$ and the density of states $g(\\omega)$ by the relation\n", + "$$\n", + " \\gamma = 2\\pi g(\\Omega) |\\kappa_0|^2.\n", + "$$\n", + "\n", + "\n", + "[1] *Gardiner C., Zoller P. The Quantum World of Ultra-Cold Atoms and Light Book I: Foundations of Quantum Optics.*" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The first step is to include the necessary libraries" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "using QuantumOptics\n", + "using PyPlot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then we define all relevant parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Nmodes = 16 # Number of electromagnetic field modes\n", + "Nfock = 1 # Maximum number of photons in the electromagnetic modes\n", + "γ = 1 # Decay rate in the master equation\n", + "κ0 = 0.5 # Atom-field interaction strength for a single mode\n", + "g = γ/(2π*abs2(κ0)) # Density of states (Fixed implicitly by κ0 and γ)\n", + "Δ_cutoff = (Nmodes-1)/(2*g) # Frequency cutoff of numerically included modes\n", + "Δ = linspace(-Δ_cutoff, Δ_cutoff, Nmodes); # Frequencies of numerically included modes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "choose appropriate bases" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "b_fock = FockBasis(Nfock)\n", + "b_space = tensor([b_fock for i=1:Nmodes]...)\n", + "b_atom = SpinBasis(1//2)\n", + "b = b_atom ⊗ b_space;" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "and create the necessary operators" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Atom operators\n", + "σ⁺ = sigmap(b_atom)\n", + "σ⁻ = sigmam(b_atom)\n", + "\n", + "# Single mode operators\n", + "a = destroy(b_fock)\n", + "aᵗ = create(b_fock)\n", + "n = number(b_fock)\n", + "\n", + "# Hamiltonian\n", + "H = SparseOperator(b)\n", + "for k=1:Nmodes\n", + " H += κ0*embed(b, [1, k+1], [σ⁺, a]) \n", + " H += κ0*embed(b, [1, k+1], [σ⁻, aᵗ])\n", + " H += Δ[k]*embed(b, k+1, n)\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, after defining an initial state, it's straight forward to calculate the dynamics." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Initial state\n", + "ψ₀ = spinup(b_atom) ⊗ tensor([fockstate(b_fock, 0) for i=1:Nmodes]...)\n", + "\n", + "# Time evolution\n", + "tspan = [0:0.01:5;]\n", + "tout, ψₜ = timeevolution.schroedinger(tspan, ψ₀, H);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We use matplotlib to compare the result of the Schrödinger equation for the whole system to the analytic result of the reduced system described by a master equation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using PyPlot\n", + "\n", + "op = embed(b, 1, sparse(dm(spinup(b_atom))))\n", + "plot(tspan, expect(op, ψₜ))\n", + "plot(tspan, exp.(-γ*tspan), \"--k\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 0.6.0", + "language": "julia", + "name": "julia-0.6" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "0.6.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}