Chapter 13 - Hydrogen Solutions.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:7da6401e22419bd3f79ff6eed92f780c9a904fa8837117c9c4e9b9ee4470be79"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Hydrogen functions\n",
      "Start with some imports fro Symbolic Python library:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sympy.physics.hydrogen import R_nl\n",
      "from sympy.functions.special.spherical_harmonics import Ynm\n",
      "from sympy import *"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Define some variables, radial, polar, azimuthal, time, and two frequencies:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "var(\"r theta phi t w1 w2\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "(r, theta, phi, t, w1, w2)"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Look at a few of the radial equations and the spherical harmonics\n",
      "Notice that instead of Ylm the name is Ynm... the arguments to the function are still the quantum numbers `l` and `m`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "R_nl(1, 0, r, 1)  # the n = 1, l = 0 radial function"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "2*exp(-r)"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Ynm(1,0,theta,phi).expand(func=True)  # the l = 0, m = 0 spherical harmonic"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "sqrt(3)*cos(theta)/(2*sqrt(pi))"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Write the equation for the $|nlm\\rangle = |100\\rangle$ state. Use the sympy method `.expand(func=True)` to convert to the actual expression. To create this state, we combine the Radial function and the Ylm function. Make sure to set n, l, and m to the correct values. The fourth argument to `R_nl` is `Z` which we set to 1 since we are talking about a 1-proton nucleus.\n",
      "\n",
      "The combination of R_nl and Ynm should look like the following (replace N, L, and M with the appropriate values):\n",
      "\n",
      "`R_nl(N, L, r, 1)*Ynm(L, M, theta, phi).expand(func=True)`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# this is the |100> state:\n",
      "psi100 = R_nl(1, 0, r, 1)*Ynm(0,0,theta,phi).expand(func=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "psi100  # check to see how it looks as an expression"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "exp(-r)/sqrt(pi)"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Integrating over all space\n",
      "Remember spherical coordinate integrals of function $f(r,\\theta,\\phi)$ over all space look like: $$\\int_0^\\infty\\int_0^\\pi\\int_0^{2\\pi}r^2\\sin(\\theta)drd\\theta d\\phi \\,\\,f(r,\\theta,\\phi)$$ so you alwasy need to add a factor of `r**2*sin(theta)` and then integrate `r` from 0 to infinity, `theta` from $0-\\pi$ and `phi` from $0-2\\pi$. As a check, you should integrate the square of the `psi100` wavefunction over all space to see that it equals 1 (i.e. it is normalized)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "integrate(r**2*sin(theta) * (psi100)**2 ,(r,0,oo),(theta,0,pi),(phi,0,2*pi))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "1"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Now do the $|210\\rangle$ state:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "psi210 = R_nl(2, 1, r, 1)*Ynm(1,0,theta,phi).expand(func=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "psi210  # check how it looks"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "sqrt(2)*r*exp(-r/2)*cos(theta)/(8*sqrt(pi))"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note, if you compare these to listed solutions (for example at http://hyperphysics.phy-astr.gsu.edu/hbase/quantum/hydwf.html#c3) you see that there are not any factors of $a_0$. This is because the `R_nl` function is defined in units of $a_0$. $a_0$ is the Bohr Radius: http://en.wikipedia.org/wiki/Bohr_radius"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Now calculate $\\langle z \\rangle$:\n",
      "To calculate $\\langle z \\rangle$ we need to convert to spherical coordinates: $z = r\\cos\\theta$. The terms in the following integral are the $r^2\\sin\\theta$ then $z$ (in spherical coords) then the wave function squared."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "expect = integrate(r**2*sin(theta)* (r*cos(theta)) * (psi100*psi100),(r,0,oo),(theta,0,pi),(phi,0,2*pi))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "expect"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "0"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "No surprise, the average z position of the electron in the hydrogen atom is 0."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Now for problem 13.21\n",
      "find $\\langle z \\rangle(t)$. Use the same integral, but add a time-dependent piece to each term in the wavefunction, add them together and multiply by the complex conjugate."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "psi = 1/sqrt(2)*(psi100*exp(1j*w1*t) + psi210*exp(1j*w2*t))\n",
      "psi_conj = 1/sqrt(2)*(psi100*exp(-1j*w1*t) + psi210*exp(-1j*w2*t))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "expect2 = integrate(r**2*sin(theta)* (r*cos(theta)) * psi*psi_conj,(r,0,oo),(theta,0,pi),(phi,0,2*pi))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(re(expect2.subs({w2:2, w1:1})),(t,0,10))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 24,
       "text": [
        "<sympy.plotting.plot.Plot at 0x108eee190>"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We need to interpret this result. First you should show that this expression is simply a constant amplitude factor times $\\cos((w2-w1)t)$, in other words $\\langle z \\rangle$ oscillates at frequency `w2-w1`."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Your assignment:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Explore other combinations of states and draw conclusions about the z behavior from the results. You may not be able to get these expressions to simplify, but the important thing is to look for the time dependence and simplify that part."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- Does $\\langle z \\rangle$ oscillate for any combination of two Hydrogen states $|nlm\\rangle$?\n",
      "- Are there restrictions on what n values give oscillating $\\langle z \\rangle$ expressions? (hint, to keep it simple, always let one state be the n=1 state)\n",
      "- How does $\\langle z \\rangle$ change with different l and m values are used in the state?\n",
      "\n",
      "Hints for interpreting your results:\n",
      "- What are the relavant frequencies in your expression for $\\langle z \\rangle$ and why?\n",
      "- Simplify one of your $\\langle z \\rangle$ expressions and write the time dependence in terms of the frequencies w2 and w1."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "psi320 = R_nl(3, 2, r, 1)*Ynm(2,0,theta,phi).expand(func=True)\n",
      "psi2 = 1/sqrt(2)*(psi100*exp(1j*w1*t) + psi310*exp(1j*w2*t))\n",
      "psi2_conj = 1/sqrt(2)*(psi100*exp(-1j*w1*t) + psi310*exp(-1j*w2*t))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "expect3 = integrate(r**2*sin(theta)* (r*cos(theta)) * psi2*psi2_conj,(r,0,oo),(theta,0,pi),(phi,0,2*pi))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "expect3"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 34,
       "text": [
        "2*pi*(27*sqrt(2)*exp(-1.0*I*t*w1)*exp(1.0*I*t*w2)/(512*pi) + 27*sqrt(2)*exp(1.0*I*t*w1)*exp(-1.0*I*t*w2)/(512*pi))"
       ]
      }
     ],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(re(expect3.subs({w2:2, w1:1})),(t,0,10))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 29,
       "text": [
        "<sympy.plotting.plot.Plot at 0x1059f0490>"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}