Project_4_Instructions.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "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.6.4"
    },
    "colab": {
      "name": "Project 4 - Instructions",
      "provenance": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Tat9SUy08Kqg",
        "colab_type": "text"
      },
      "source": [
        "# Project description\n",
        "\n",
        "Submit a Jupyter notebook containing the following sections, clearly denoted with markdown headings.\n",
        "\n",
        "### Abstract\n",
        "\n",
        "A general description of the problem you studied and of your results.\n",
        "\n",
        "### Physical System and Model\n",
        "\n",
        "A description of the physical system you are studying and the mathematical model being used to study it.  Your description should be sufficiently descriptive and mathematically detailed that someone who's taken the upper-division core physics courses should be able to follow it, even if they don't have much knowledge of the system and model.\n",
        "\n",
        "### Algorithms\n",
        "\n",
        "A description of any algorithm(s) you are using to solve the problem at hand, and their implementation in the context of that problem.  This doesn't need to be as detailed as providing pseudocode but should be sufficiently detailed that it's clear how the algorithm(s) are being used to understand the mathematical model of the physical system being studied.  You should also incluide a motivation for taking an algorithmic approach as opposed to using pen-and-paper methods.\n",
        "\n",
        "### Code\n",
        "\n",
        "A python code cell that includes any main functions you use to obtain your results.  Note that there can be lots of other coding cells elsewhere in the notebook, but this cell should be where all the core code is written.\n",
        "\n",
        "### Results\n",
        "\n",
        "A description of your numerical solution to the problem along with a discussion of how you can be confident that you can trust your results.  In particular, include any code tests that you believe are relevant and a discussion of how you know that the output of your code for the system at hand makes qualitative (and perhaps quantitative) sense.\n",
        "\n",
        "### References\n",
        "\n",
        "All resources you use should be cited.  This includes but is not limited to, books, papers, and websites."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "S-kySxGL8Kqi",
        "colab_type": "text"
      },
      "source": [
        "# Choose your own project adventure.\n",
        "\n",
        "## Option 1: Precession of the perihelion of Mercury\n",
        "\n",
        "According to Newtonian gravitation, if there were only a single planet orbiting the Sun, it would orbit in an ellipse.  This is a fact we verified numerically in the first project.  The point of closest approach of such a planet is called the **perihelion** of its orbit.\n",
        "\n",
        "In our solar system, there isn't just a single planet orbiting the sun, there are 8.  One might wonder what effect the gravitational interactions between these planets have on their orbits.  If you're an astronomer, you don't have to wonder; you can't just go out and observe the planetary orbits.  If you do this precisely enough, you'll find that Mercury's perihelion precesses -- it doesn't stay stationary but rotates at some angular velocity.\n",
        "\n",
        "Observations reveal that the total precession is about 574 arcseconds per century.  The dominant contribution to this precession, more than 90% of it, is due to the gravitational interactions with other planets.  The rest of the effect is primarily due to corrections to Newtonian gravitation from general relativity.\n",
        "\n",
        "Show, using numerical simulation, that most of Mercury's perihelion precession can be accounted for by Newtonian gravitational interactions with the other planets by simulating the entire solar system for sufficiently long that the precession becomes measurable.  Do not assume that the sun is stationary -- this approximate fact should emerge from you simulation.\n",
        "\n",
        "If you're familiar with general relativity, you might try to also take on the challenge of showing numerically that including general relativistic effects allows one to account for the remainder of the precession not accounted for by Newtonian gravitation.\n",
        "\n",
        "## Option 2: Statistical mechanics of the quantum harmonic oscillator using path integral MCMC\n",
        "\n",
        "When the quantum harmonic oscillator is put into contact with a heat bath at absolute temperature $T$, there is some probability $p_n(T)$ that the system will occupy each of its energy eigenstates proportional to the Boltzmann factor associated with that state.  This can be used to determine an expression for the average energy the system will have in equilibrium:\n",
        "\n",
        "\\begin{align}\n",
        "    \\langle \\hat H\\rangle = \\sum_{n=0}^\\infty p_n E_n = -\\frac{1}{Z}\\frac{\\partial Z}{\\partial \\beta}\n",
        "\\end{align}\n",
        "\n",
        "We have already seen one way to numerically determine this ensemble average energy using MCMC in the case of the Ising model: we can generate a sequence of the system's states using Metropolis-Hastings and use this sequence of states to compute an approximation to the ensemble average value of any observable we desire such as the equilibrium average energy.\n",
        "\n",
        "It turns out that there is another way to do this using Feynman's path integral formulation of quantum mechanics and MCMC -- so-called **path integral monte carlo**.  See the appendix below for notes on the path integral that you may find useful.\n",
        "\n",
        "Compute and plot $\\langle \\hat H\\rangle$ versus $T$ (or $\\beta$) for the quantum harmonic oscillator analytically first by simply performing the sum over states written above (or some equivalent procedure you might be familiar with from stat mech), then do the same using path integral monte carlo.  Make sure your results agree reasonably well.  \n",
        "\n",
        "Advice: when I was first doing this, I found it useful to make sure I understood how to sample from a familiar continuous probability distribution on $\\mathbb R^N$ before tackling the path integral.  A good choice is a multi-variate Gaussian because you can check your sampling results against built-in Gaussian samplers in scientific Python packages e.g. numpy/scipy.\n",
        "\n",
        "You may also find the following papers useful, but you may consult any resources you find useful as long as you cite them:\n",
        "\n",
        "[The harmonic oscillator at finite temperature using path integrals](https://aapt.scitation.org/doi/abs/10.1119/1.15737), Larsen, Ravndal\n",
        "[Path integral for the quantum harmonic oscillator using elementary methods](https://aapt.scitation.org/doi/abs/10.1119/1.18896), Cohen\n",
        "\n",
        "## Option 3: Build your own\n",
        "\n",
        "You may find an interesting physics problem that lends itself to numerical solution that you'd like to solve and build your own project on that problem.  If you want to do this, Josh recommends running it by him first."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "m5c2ZWXcFnNH",
        "colab_type": "text"
      },
      "source": [
        "# Appendix: Path integral details"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0CxbatZMF0-L",
        "colab_type": "text"
      },
      "source": [
        "## Preliminary facts and conventions\n",
        "\n",
        "Quantum Stat mech\n",
        "\n",
        "\\begin{align}\n",
        "    Z(\\beta) = \\mathrm{tr}(e^{-\\beta \\hat H})\n",
        "\\end{align}\n",
        "\n",
        "Position space representation of the trace of an operator for one-dimensional system\n",
        "\n",
        "\\begin{align}\n",
        "    \\mathrm{tr}\\,\\hat O = \\int dq \\langle q|\\hat O|q\\rangle\n",
        "\\end{align}\n",
        "\n",
        "Operator exponential facts\n",
        "\n",
        "\\begin{align}\n",
        "    e^{\\hat A + \\hat B} = e^{\\hat A}e^{\\hat B}, \\qquad [\\hat A, \\hat B] = 0.\n",
        "\\end{align}\n",
        "\n",
        "More generally for $\\epsilon > 0$\n",
        "\n",
        "\\begin{align}\n",
        "    e^{\\epsilon(\\hat A + \\hat B)} = e^{\\epsilon \\hat A}e^{\\epsilon \\hat B} + O(\\epsilon^2)\n",
        "\\end{align}\n",
        "\n",
        "where the $O(\\epsilon^2)$ terms depend on the commutator $[\\hat A, \\hat B]$ in such a way that they vanish when the commutator vanishes.\n",
        "\n",
        "Quantum mechanics representations of the identity in position and momentum:\n",
        "\n",
        "\\begin{align}\n",
        "    \\hat I = \\int dq\\, |q\\rangle\\langle q|, \\qquad \\hat I = \\int \\frac{dp}{2\\pi\\hbar}\\, |p\\rangle\\langle p|\n",
        "\\end{align}\n",
        "\n",
        "Inner products of position and momentum eigenstates\n",
        "\n",
        "\\begin{align}\n",
        "    \\langle q|p\\rangle = e^{iqp/\\hbar}, \\qquad \\langle p|q\\rangle = e^{-ipq/\\hbar}\n",
        "\\end{align}\n",
        "\n",
        "Analytic functions of operators actions eigenstates:\n",
        "\n",
        "\\begin{align}\n",
        "    f(\\hat Q)|q\\rangle = f(q)|q\\rangle.\n",
        "\\end{align}\n",
        "\n",
        "Fourier transform of a Gaussian is a Gaussian:\n",
        "\n",
        "\\begin{align}\n",
        "    \\int \\frac{e^{-\\frac{k^2}{2\\sigma^2}}}{\\sqrt{2\\pi\\sigma^2}} e^{ikx} \\, dk\n",
        "    &= e^{-x^2/(2/\\sigma^2)}\n",
        "\\end{align}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RNsfOrnvGSpT",
        "colab_type": "text"
      },
      "source": [
        "## Step 1: Quantum partition function as position-space trace\n",
        "\n",
        "The partition function is\n",
        "\n",
        "$$\n",
        "    Z(\\beta) = \\mathrm{tr}(e^{-\\beta \\hat H})\n",
        "$$\n",
        "\n",
        "If $\\{|n\\rangle\\}$ is an ONB then for any operator (trace class) $\\hat O$ we have:\n",
        "\n",
        "\\begin{align}\n",
        "  \\mathrm{tr}(\\hat O) \n",
        "  &= \\sum_n\\langle n|\\hat O|n\\rangle \\\\\n",
        "  &= \\sum_n \\langle n|\\int dq |q\\rangle\\langle q| \\hat O\\int dq' |q'\\rangle\\langle q'||n\\rangle\\\\\n",
        "  &= \\int dq\\int dq'\\langle q|\\hat O|q'\\rangle \\sum_n\\langle q'|n\\rangle \\langle n|q\\rangle \\\\\n",
        "  &= \\int dq\\int dq'\\langle q|\\hat O|q'\\rangle \\delta(q-q') \\\\\n",
        "  &= \\int dq \\langle q|\\hat O|q\\rangle\n",
        "\\end{align}\n",
        "\n",
        "It follows that for the partition function we have\n",
        "\n",
        "$$\n",
        "    Z(\\beta) = \\int dq \\,\\langle q|e^{-\\beta \\hat H}|q\\rangle\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KGFIDBaLGh5T",
        "colab_type": "text"
      },
      "source": [
        "## Step 2: Partition function as a multidimensional integral by splitting\n",
        "\n",
        "We recall that if $[\\hat A, \\hat B] = 0$ then\n",
        "\n",
        "\\begin{align}\n",
        "    e^{\\hat A + \\hat B} = e^{\\hat A}e^{\\hat B}.\n",
        "\\end{align}\n",
        "\n",
        "In particular\n",
        "\n",
        "\\begin{align}\n",
        "    e^{\\hat O + \\hat O} = e^{\\hat O}e^{\\hat O}\n",
        "\\end{align}\n",
        "\n",
        "Thus\n",
        "\n",
        "\\begin{align}\n",
        "    e^{-\\beta \\hat H}\n",
        "    &= e^{-\\left(\\frac{\\beta}{N} + \\frac{\\beta}{N} + \\cdots + \\frac{\\beta}{N}\\right) \\hat H} \\\\\n",
        "    &= e^{-\\frac{\\beta}{N} \\hat H}e^{-\\frac{\\beta}{N} \\hat H}\\cdots e^{-\\frac{\\beta}{N} \\hat H}\n",
        "\\end{align}\n",
        "\n",
        "It follows that:\n",
        "\n",
        "\\begin{align}\n",
        "    Z(\\beta) \n",
        "    &= \\int dq \\,\\langle q|e^{-\\beta \\hat H}|q\\rangle \\\\\n",
        "    &= \\int dq_1 \\,\\langle q_1|e^{-\\frac{\\beta}{N} \\hat H}e^{-\\frac{\\beta}{N} \\hat H}\\cdots e^{-\\frac{\\beta}{N} \\hat H}|q_1\\rangle \\\\\n",
        "    &= \\int dq_1\\int dq_2\\cdots \\int dq_N \\,\\langle q_1|e^{-\\frac{\\beta}{N} \\hat H}|q_N\\rangle\\langle q_N|e^{-\\frac{\\beta}{N} \\hat H}|q_{N-1}\\rangle\\langle q_{N-1}|\\cdots |q_2\\rangle\\langle q_2| e^{-\\frac{\\beta}{N} \\hat H}|q_1\\rangle \\\\\n",
        "    &= \\int d^N q \\,\\langle q_1|e^{-\\frac{\\beta}{N} \\hat H}|q_N\\rangle\\langle q_N|e^{-\\frac{\\beta}{N} \\hat H}|q_{N-1}\\rangle\\langle q_{N-1}|\\cdots |q_2\\rangle\\langle q_2| e^{-\\frac{\\beta}{N} \\hat H}|q_1\\rangle\n",
        "\\end{align}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Bpen8gLbGr90",
        "colab_type": "text"
      },
      "source": [
        "## Step 3: Operator exponential splitting for small arguments\n",
        "\n",
        "Suppose that $\\epsilon>0$ is a real number, then notice that\n",
        "\n",
        "\\begin{align}\n",
        "    e^{\\epsilon(\\hat A + \\hat B)}\n",
        "    &= \\hat I + \\epsilon(\\hat A + \\hat B) + \\frac{\\epsilon^2}{2!}(\\hat A + \\hat B)^2 + \\cdots\n",
        "\\end{align}\n",
        "\n",
        "and\n",
        "\n",
        "\\begin{align}\n",
        "    e^{\\epsilon\\hat A}e^{\\epsilon\\hat B} \n",
        "    &= (\\hat I + \\epsilon \\hat A + O(\\epsilon^2))(\\hat I + \\epsilon \\hat B + O(\\epsilon^2)) \\\\\n",
        "    &= \\hat I + \\epsilon\\hat A + \\epsilon \\hat B + O(\\epsilon^2) \\\\\n",
        "    &= \\hat I + \\epsilon(\\hat A + \\hat B) + O(\\epsilon^2)\n",
        "\\end{align}\n",
        "\n",
        "It follows that\n",
        "\n",
        "\\begin{align}\n",
        "    e^{\\epsilon(\\hat A + \\hat B)} = e^{\\epsilon\\hat A}e^{\\epsilon\\hat B} + O(\\epsilon^2)\n",
        "\\end{align}\n",
        "\n",
        "Therefore, if $\\epsilon$ is small, we have approximately\n",
        "\n",
        "\\begin{align}\n",
        "    e^{\\epsilon(\\hat A + \\hat B)} \\approx e^{\\epsilon\\hat A}e^{\\epsilon\\hat B}\n",
        "\\end{align}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "S-uuNc-wGxv8",
        "colab_type": "text"
      },
      "source": [
        "## Step 4: Evaluation of each split factor in the partition function\n",
        "\n",
        "Define $\\epsilon = \\beta/N$ and suppose then\n",
        "\n",
        "\\begin{align}\n",
        "    \\hat H \n",
        "    &= \\hat T + \\hat V \\\\\n",
        "    &= \\frac{1}{2m}\\hat P^2 + V(\\hat Q)\n",
        "\\end{align}\n",
        "\n",
        "\\begin{align}\n",
        "    \\langle q_{i+1}|e^{-\\frac{\\beta}{N}\\hat H}|q_{i}\\rangle\n",
        "    &= \\langle q_{i+1}|e^{-\\epsilon\\hat H}|q_{i}\\rangle \\\\\n",
        "    &= \\langle q_{i+1}|e^{-\\epsilon(\\hat T + \\hat V)}|q_{i}\\rangle \\\\\n",
        "    &\\approx \\langle q_{i+1}|e^{-\\epsilon\\hat T}e^{-\\epsilon\\hat V}|q_{i}\\rangle \\\\\n",
        "    &= \\int \\frac{dp}{2\\pi\\hbar} \\langle q_{i+1}|e^{-\\epsilon\\hat T}|p\\rangle\\langle p|e^{-\\epsilon V(\\hat Q)}|q_{i}\\rangle \\\\\n",
        "    &= \\int \\frac{dp}{2\\pi\\hbar} \\,e^{-\\epsilon p^2/2m}e^{-\\epsilon V(q_{i})}\\langle q_{i+1}|p\\rangle\\langle p|q_{i}\\rangle \\\\\n",
        "    &= \\int \\frac{dp}{2\\pi\\hbar} \\,e^{-\\epsilon p^2/2m}e^{-\\epsilon V(q_{i})}e^{iq_{i+1}p/\\hbar}e^{-iq_ip/\\hbar} \\\\\n",
        "    &= \\frac{e^{-\\epsilon V(q_{i})}}{2\\pi\\hbar}\\int dp \\,e^{-\\epsilon p^2/2m}e^{i(q_{i+1} - q_i)p/\\hbar} \\\\\n",
        "    &= \\frac{e^{-\\epsilon V(q_{i})}}{2\\pi\\hbar} \\sqrt{\\frac{2\\pi m}{\\epsilon}}e^{-m(q_{i+1} - q_i) ^2/(2\\epsilon\\hbar ^2)} \\\\\n",
        "    &= \\frac{e^{-\\epsilon V(q_{i})}}{(2\\pi)^{1/2}} \\sqrt{\\frac{ m}{\\epsilon\\hbar^2}}\\exp\\left(-\\frac{1}{2}\\frac{(q_{i+1} - q_i)^2}{\\epsilon\\hbar ^2/m}\\right) \\\\\n",
        "    &= \\frac{1}{(2\\pi)^{1/2}A(\\epsilon)} \\exp\\left[-\\frac{1}{2}\\left(\\frac{q_{i+1} - q_i}{A(\\epsilon)}\\right)^2-\\epsilon V(q_{i})\\right]\n",
        "\\end{align}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "03R1cvGPHE4r",
        "colab_type": "text"
      },
      "source": [
        "## Step 5: Putting it all together\n",
        "\n",
        "The $N$-discretized partition function of a quantum particle of mass $m$ moving in one dimension at inverse temperature $\\beta$ subject to the potential $V(q)$ is\n",
        "\n",
        "$$\n",
        "    Z_N(\\beta) = \\int_{\\mathbb R^N}\\frac{d^Nq}{(2\\pi)^{N/2}A(\\epsilon)^N}\\,\\exp\\left(-\\sum_{i=1}^N\\left[\\frac{1}{2}\\left(\\frac{q_{i+1}-q_i}{A(\\epsilon)}\\right)^2 + \\epsilon V(q_i)\\right]\\right), \\qquad A(\\epsilon) \\equiv \\left(\\frac{\\epsilon\\hbar^2}{m}\\right)^{1/2}, \\qquad \\epsilon \\equiv \\frac{\\beta}{N}, \\qquad q_{N+1} \\equiv q_1\n",
        "$$\n",
        "\n",
        "Notice that interestingly, the partition function only depends $\\beta$ through the ratio $\\epsilon = \\beta/N$.  Letting, $u_i = q_i/A(\\epsilon)$, the measure becomes independent of $\\beta$, as does the \"kinetic term\" in the exponential:\n",
        "\n",
        "$$\n",
        "    Z_N(\\beta) = \\int_{\\mathbb R^N}\\frac{d^Nu}{(2\\pi)^{N/2}}\\,\\exp\\left(-\\sum_{i=1}^N\\left[\\frac{1}{2}\\left(u_{i+1}-u_i\\right)^2 + \\epsilon V\\Big(A(\\epsilon)u_i\\Big)\\right]\\right)\n",
        "$$\n",
        "\n",
        "The ensemble average energy is\n",
        "\n",
        "$$\n",
        "  \\langle \\hat H\\rangle_N \n",
        "  = - \\frac{\\partial}{\\partial\\beta} \\ln Z_N \n",
        "  = - \\frac{1}{Z_N} \\frac{\\partial Z_N}{\\partial\\beta} \n",
        "  = \\frac{\\int_{\\mathbb R^N}d^Nu\\,\\frac{\\partial}{\\partial\\beta}\\left(\\sum_{i=1}^N\\epsilon V(A(\\epsilon) u_i)\\right) \\exp\\left(-\\sum_{i=1}^N\\left[\\frac{1}{2}\\left(u_{i+1}-u_i\\right)^2 + \\epsilon V\\Big(A(\\epsilon)u_i\\Big)\\right]\\right)}{\\int_{\\mathbb R^N}d^Nu\\,\\exp\\left(-\\sum_{i=1}^N\\left[\\frac{1}{2}\\left(u_{i+1}-u_i\\right)^2 + \\epsilon V\\Big(A(\\epsilon)u_i\\Big)\\right]\\right)}\n",
        "$$\n",
        "\n",
        "Now notice that\n",
        "\n",
        "$$\n",
        "  \\frac{\\partial}{\\partial\\beta} = \\frac{\\partial}{\\partial\\epsilon}\\frac{\\partial\\epsilon}{\\partial\\beta} = \\frac{1}{N}\\frac{\\partial}{\\partial\\epsilon}\n",
        "$$\n",
        "\n",
        "and therefore the term in the integrand of the numerator is\n",
        "\n",
        "$$\n",
        "  \\frac{\\partial}{\\partial\\beta}\\left(\\sum_{i=1}^N\\epsilon V(A(\\epsilon) u_i)\\right) = \\frac{1}{N}\\frac{\\partial}{\\partial\\epsilon}\\left(\\sum_{i=1}^N\\epsilon V(A(\\epsilon) u_i)\\right) = \\frac{1}{N}\\sum_{i=1}^N \\left[V(A(\\epsilon)u_i) + \\epsilon V'(A(\\epsilon) u_i)A'(\\epsilon)u_i\\right]\n",
        "  = \\frac{1}{N}\\sum_{i=1}^N \\left[V(A(\\epsilon) u_i) + \\frac{1}{2}(A(\\epsilon)u_i)V'(A(\\epsilon) u_i)\\right]\n",
        "$$\n",
        "\n",
        "Plugging this back into the expression for the ensemble average energy, and reverting the integrals to the old variables $q_i$ and using the notation $\\vec q = (q_1, q_2, \\dots, q_N)$ gives\n",
        "\n",
        "$$\n",
        "  \\langle \\hat H\\rangle_N = \\frac{\\int_{\\mathbb R^N}d^Nq\\,\\Phi(\\vec q) \\,e^{-U(\\vec q)}}{\\int_{\\mathbb R^N}d^Nq\\,e^{-U(\\vec q)}}, \\qquad \\Phi(\\vec q) = \\frac{1}{N}\\sum_{i=1}^N\\left[V(q_i) + \\frac{1}{2}q_i V'(q_i)\\right], \\qquad U(\\vec q) = \\sum_{i=1}^N\\left[\\frac{m}{2(\\beta/N)\\hbar^2}\\left(q_{i+1}-q_i\\right)^2 + (\\beta/N) V(q_i)\\right]\n",
        "$$\n",
        "\n",
        "This is now presented in a way that MCMC can be applied for a given potential since we have written the desired quantity as the expectation value of a certain function $\\Phi$ on the space of vectors $\\vec q$.  To apply the algorithm, one needs to sample points from this space according to the un-normalized probability distribution $e^{-U(\\vec q)}$ using the Metropolis-Hastings algorithm.  Note that the paramaters/functions which need to be controlled are\n",
        "\n",
        "$$\n",
        "  N, \\beta, V, m, \\hbar\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SHiv4r4GDDoB",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}