Derivation of numerical master equation example.

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
+}