Chapter 13 - Hydrogen Solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d6a7d5af",
   "metadata": {},
   "source": [
    "#Hydrogen functions\n",
    "Start with some imports fro Symbolic Python library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2b770da4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy.physics.hydrogen import R_nl\n",
    "from sympy.functions.special.spherical_harmonics import Ynm\n",
    "from sympy import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f66c789",
   "metadata": {},
   "source": [
    "Define some variables, radial, polar, azimuthal, time, and two frequencies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "81b009dc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(r, theta, phi, t, w1, w2)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var(\"r theta phi t w1 w2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6711d13",
   "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",
   "execution_count": 3,
   "id": "6c8d3b23",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2*exp(-r)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R_nl(1, 0, r, 1)  # the n = 1, l = 0 radial function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c2d5dfa9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sqrt(3)*cos(theta)/(2*sqrt(pi))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ynm(1,0,theta,phi).expand(func=True)  # the l = 0, m = 0 spherical harmonic"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11845571",
   "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",
   "execution_count": 5,
   "id": "1e84ea7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is the |100> state:\n",
    "psi100 = R_nl(1, 0, r, 1)*Ynm(0,0,theta,phi).expand(func=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "61f59024",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "exp(-r)/sqrt(pi)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psi100  # check to see how it looks as an expression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d3eeee1",
   "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",
   "execution_count": 7,
   "id": "661634bf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(r**2*sin(theta) * (psi100)**2 ,(r,0,oo),(theta,0,pi),(phi,0,2*pi))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01ab7baa",
   "metadata": {},
   "source": [
    "## Now do the $|210\\rangle$ state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ff69d290",
   "metadata": {},
   "outputs": [],
   "source": [
    "psi210 = R_nl(2, 1, r, 1)*Ynm(1,0,theta,phi).expand(func=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d6706968",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sqrt(2)*r*exp(-r/2)*cos(theta)/(8*sqrt(pi))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psi210  # check how it looks"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1abc84d0",
   "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",
   "id": "3f9595ca",
   "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",
   "execution_count": 10,
   "id": "8d3aa975",
   "metadata": {},
   "outputs": [],
   "source": [
    "expect = integrate(r**2*sin(theta)* (r*cos(theta)) * (psi100*psi100),(r,0,oo),(theta,0,pi),(phi,0,2*pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "13c7fd2c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expect"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aabd0e6f",
   "metadata": {},
   "source": [
    "No surprise, the average z position of the electron in the hydrogen atom is 0."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "853dd5fa",
   "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",
   "execution_count": 12,
   "id": "039bb325",
   "metadata": {},
   "outputs": [],
   "source": [
    "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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "477da873",
   "metadata": {},
   "outputs": [],
   "source": [
    "expect2 = integrate(r**2*sin(theta)* (r*cos(theta)) * psi*psi_conj,(r,0,oo),(theta,0,pi),(phi,0,2*pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "2bd1cf54",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<sympy.plotting.plot.Plot at 0x108eee190>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(re(expect2.subs({w2:2, w1:1})),(t,0,10))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3a4ff43",
   "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",
   "id": "bb0c4935",
   "metadata": {},
   "source": [
    "## Your assignment:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c046d6f",
   "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",
   "id": "77def1f0",
   "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",
   "execution_count": 30,
   "id": "4f61b93a",
   "metadata": {},
   "outputs": [],
   "source": [
    "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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "37d4bc9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "expect3 = integrate(r**2*sin(theta)* (r*cos(theta)) * psi2*psi2_conj,(r,0,oo),(theta,0,pi),(phi,0,2*pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "42aafb51",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "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))"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expect3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "a7047edd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<sympy.plotting.plot.Plot at 0x1059f0490>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(re(expect3.subs({w2:2, w1:1})),(t,0,10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b489ee4a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}