Polaron_effective_mass_theories.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The work below relates to the effective mass theory sections of the 2017 publication _Slow cooling of hot polarons in halide perovskite solar cells_. The notebook was written by Jarvist M. Frost."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0e-10"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Physics Constants for electronic structure + stat physics\n",
    "# defined as const to help with type stability\n",
    "# 2016-07 - Updated accuracy from https://github.com/DataWookie/PhysicalConstants.jl\n",
    "const c = 2.99792458e8;                            # m / s  # speed of light in this universe\n",
    "const h= 6.62606896e-34;                           # kg m2 / s \n",
    "const hbar = const ħ = 1.05457162825e-34;          # kg m2 / s \n",
    "\n",
    "const eV = const q = const ElectronVolt = 1.602176487e-19;                         # kg m2 / s2 \n",
    "const MassElectron = const melectron = 9.10938188e-31;                          # kg\n",
    "const ε0 = const ε_0 = const VacuumPermittivity = 8.854187817e-12;                   # A2 s4 / kg m3 \n",
    "\n",
    "k_b=GSLCONSTMKSABOLTZMANN = 1.3806504e-23 # kg m2 / K s2\n",
    "\n",
    "# Units\n",
    "Å=1E-10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.12"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Material specifics - for MAPbI3\n",
    "d=6.2Å # lattice spacing, MAPbI\n",
    "\n",
    "# Dielectric constants\n",
    "ε_Inf=4.5\n",
    "ε_S=24.1\n",
    "\n",
    "f=2.25E12 # Characteristic phonon mode frequency\n",
    "ω=2*π*f\n",
    "\n",
    "# Bare-band effective masses.\n",
    "# From QSGW calc, MAPI\n",
    "m_eff_e=0.12\n",
    "m_eff_h=0.15\n",
    "m_eff = m_eff_e # default to using electron effective mass. No particular reason why - but you've got to choose one!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " Cross check vs. literature values.\n",
      " NaCl Frohlich paper α=5.262322769751378\n",
      " CdTe  α=0.3505417171832051\n",
      " GaAs  α=0.07081920905063181\n",
      "\n",
      " MAPI  4.5, 24.1, 2.25THz (~75 cm^-1) ; α=2.3939407113094093"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2.3939407113094093"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# via Feynman 1955\n",
    "#   http://dx.doi.org/10.1103/PhysRev.97.660\n",
    "function feynmanalpha(ε_Inf,ε_S,freq,m_eff)\n",
    "    ω=freq*2*pi #frequency to angular velocity\n",
    "    # Note to self - you've introduced the 4*pi factor into the dielectric constant; \n",
    "    # This gives numeric agreement with literature values, but I'm uncertain of the justification.\n",
    "    # Such a factor also seems to be necessary for calculation of Exciton binding energies (see below).\n",
    "    # Maybe scientists in the 50s assumed that the Epsilon subsumed the 4*pi factor?\n",
    "    α=0.5* (1/ε_Inf - 1/ε_S)/(4*pi*ε_0) * (q^2/(hbar*ω)) * (2*MassElectron*m_eff*ω/hbar)^0.5\n",
    "    return (α)\n",
    "end\n",
    "\n",
    "print(\"\\n Cross check vs. literature values.\")\n",
    "print(\"\\n NaCl Frohlich paper α=\",feynmanalpha(2.3, 5.6, (4.9E13/(2*pi)), 1.0)) # should be ~about 5 (Feynman1955)\n",
    "print(\"\\n CdTe  α=\",feynmanalpha(7.1,   10.4,  5.08E12, 0.095)) #Stone 0.39 / Devreese 0.29\n",
    "print(\"\\n GaAs  α=\",feynmanalpha(10.89, 12.9,  8.46E12, 0.063)) # Devreese 0.068\n",
    "\n",
    "print(\"\\n\")\n",
    "print(\"\\n MAPI  4.5, 24.1, 2.25THz (~75 cm^-1) ; α=\",feynmanalpha(4.5,   24.1,  2.25E12,    0.12))\n",
    "\n",
    "#This value used for further calculations...\n",
    "α=feynmanalpha(ε_Inf, ε_S ,  f,    m_eff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bare mass (no polaronic...) 0.12\n",
      "Effective mass Large polaron limit, *(1+α/6) , 0.1678788142261882\n",
      "Effective mass Small polaron limit, *0.02*α^6 , 0.45174315859021225\n"
     ]
    }
   ],
   "source": [
    "#α=1.1932891264417276 \n",
    "\n",
    "println(\"Bare mass (no polaronic...) \",m_eff)\n",
    "# By Feynman approximation for polaron mass. See Feynman 1955 for original literature.\n",
    "#  Also Eqn (9.6.21); p. 955\n",
    "#  Jones and March, Theoretical solid state physics, Vol 2\n",
    "println(\"Effective mass Large polaron limit, *(1+α/6) , \",m_eff*(1+α/6))\n",
    "\n",
    "println(\"Effective mass Small polaron limit, *0.02*α^6 , \",m_eff*0.02*α^6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ε_Red: 0.18072844628861226\n",
      "Bare effective mass:0.12 , ε_Red:0.18072844628861226\n",
      "Rescaled (+α/6) m_eff: 0.1678788142261882)Large polaron radius\n",
      "Jones and March, rp = 2.7396722621139184e-9 m\n",
      "\t = 27.396722621139183 Å = 4.418826229215997 Lattice units\n",
      "Volume polaron: 8.613594990536387e-26 m (SI), 86135.94990536386 Å^3  361.4176676905939 unit cells^3\n",
      "Large Polaron binding energy: 4.781232783166302e-20 J 0.2984211054126088 eV\n"
     ]
    }
   ],
   "source": [
    "# Large polaron radius via p.942 in Theoretical Solid State Physics (Jones and March), Vol 2\n",
    "# ISBN: 0486650162\n",
    "ε_Red=(1/ε_Inf-1/ε_S)\n",
    "print(\"ε_Red: \",ε_Red,\"\\n\")\n",
    "\n",
    "# Not sure if you use dressed polaron effective mass here; or direct band structure one\n",
    "println(\"Bare effective mass:$(m_eff) , ε_Red:$(ε_Red)\")\n",
    "m_eff=m_eff*(1+α/6)\n",
    "print(\"Rescaled (+α/6) m_eff: $(m_eff))\")\n",
    "\n",
    "# Eqn (9.3.7); p.942; \n",
    "#  Jones and March, Theoretical solid state physics, Vol 2\n",
    "rp=2*pi^2*hbar^2*ε_0*ε_Red^-1/(q*q*MassElectron*m_eff)\n",
    "\n",
    "println(\"Large polaron radius\")\n",
    "println(\"Jones and March, rp = \",rp,\" m\")\n",
    "println(\"\\t = \",rp/Å,\" Å = \",rp/d,\" Lattice units\")\n",
    "\n",
    "V=(4/3)*π * rp^3\n",
    "println(\"Volume polaron: $V m (SI), \",V/Å^3, \" Å^3  \",V/d^3,\" unit cells^3\")\n",
    "\n",
    "# Eqn (9.3.8); p.942; \n",
    "#  Jones and March, Theoretical solid state physics, Vol 2\n",
    "W=(0.25*q^2/rp)*ε_Red/ε_0\n",
    "println(\"Large Polaron binding energy: $W J \",W/q,\" eV\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rfb (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This from T.D.Schultz 1959, Physical Review, Vol 116, Number 3, pp526--543\n",
    "# Schultz (2.5a)\n",
    "#rfa(α)=(3/(0.44*α))^(1/2)\n",
    "# After re-derivation from Feynman1955 ( see 8-8-2017; Yellow 2017.B Notebook pp.33-34 ), \n",
    "# realise this is actually 4/9=0.4444...\n",
    "rfa(α)=(3/((4/9)*α))^(1/2)\n",
    "# Also spot this, in the variational form from Feynman in Schultz, you need a term of \n",
    "# meff in the scaling factor, as the polaron state is expressed in terms of the effective mass\n",
    "\n",
    "# Schultz (2.5b)\n",
    "rfb(α)=3*(π/2)^(1/2)*α\n",
    "#rfb(α)=3*(π/(2*α))^(1/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Let's try and reproduce Schultz little table below (2.5b).\n",
      "  α=3, 1.42   ||  α=5, 1.00  ||  α=7, 0.748  ||  α=9, 0.557  ||  α=11, 0.443  \n",
      "α=3 Schultz(2.5a) α→0 rf=1.50 (2mω)^1/2 Schultz (2.5b) α→∞ rf=11.28 (2mω)^1/2\n",
      "α=5 Schultz(2.5a) α→0 rf=1.16 (2mω)^1/2 Schultz (2.5b) α→∞ rf=18.80 (2mω)^1/2\n",
      "α=7 Schultz(2.5a) α→0 rf=0.98 (2mω)^1/2 Schultz (2.5b) α→∞ rf=26.32 (2mω)^1/2\n",
      "α=9 Schultz(2.5a) α→0 rf=0.87 (2mω)^1/2 Schultz (2.5b) α→∞ rf=33.84 (2mω)^1/2\n",
      "α=11 Schultz(2.5a) α→0 rf=0.78 (2mω)^1/2 Schultz (2.5b) α→∞ rf=41.36 (2mω)^1/2\n",
      "Unsure of origin of disagreement. (small α is by far the closest) Schultz table may be numerical integration (of variational problem) results, rather than using the asymtotic Feynman solutions.\n"
     ]
    }
   ],
   "source": [
    "println(\"Let's try and reproduce Schultz little table below (2.5b).\")\n",
    "println(\"  α=3, 1.42   ||  α=5, 1.00  ||  α=7, 0.748  ||  α=9, 0.557  ||  α=11, 0.443  \")\n",
    "for α=3:2:11\n",
    "    print(\"α=$α \")\n",
    "    @printf(\"Schultz(2.5a) α→0 rf=%.2f (2mω)^1/2 \",rfa(α))\n",
    "    @printf(\"Schultz (2.5b) α→∞ rf=%.2f (2mω)^1/2\\n\",rfb(α))\n",
    "end\n",
    "println(\"Unsure of origin of disagreement. (small α is by far the closest) Schultz table may be numerical integration (of variational problem) results, rather than using the asymtotic Feynman solutions.\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "MAPI(meff=0.12): α=2.3939407113094093\n",
      "v=3.5319868247354242 w=3 by the Feynman1955 small-α expansion\n",
      "Schultz (2.5a) α→0 rf=1.679172028574159 (internal units) = 2.952068697325166e-9 (SI, m)\n",
      "Schultz (2.5b) α→∞ rf=9.001079212137622 (internal units) = 1.5824348983980366e-8 (SI, m)\n",
      "Scale of internal units: 1.7580501860978867e-9 (SI,m)\n",
      "\n",
      "MAPI(meff=0.15): α=2.6765070822960193\n",
      "v=3.5947793516213373 w=3 by the Feynman1955 small-α expansion\n",
      "Schultz (2.5a) α→0 rf=1.588062856096852 (internal units) = 3.1214326082540507e-9 (SI, m)\n",
      "Schultz (2.5b) α→∞ rf=10.063512494599985 (internal units) = 1.9780436229975453e-8 (SI, m)\n",
      "Scale of internal units: 1.965559861985515e-9 (SI,m)\n"
     ]
    }
   ],
   "source": [
    "for meff in [m_eff_e,m_eff_h]\n",
    "    α=feynmanalpha(ε_Inf, ε_S ,  f,    meff)\n",
    "    println(\"\\nMAPI(meff=$meff): α=$α\")\n",
    "\n",
    "    scale=sqrt(2*meff*melectron*ω )\n",
    "    \n",
    "    w=3 ; v=w+(2/9)*α\n",
    "    println(\"v=$v w=$w by the Feynman1955 small-α expansion\")\n",
    "    println(\"Schultz (2.5a) α→0 rf=\",rfa(α),\" (internal units) = \",rfa(α)*scale,\" (SI, m)\")\n",
    "    println(\"Schultz (2.5b) α→∞ rf=\",rfb(α),\" (internal units) = \",rfb(α)*scale,\" (SI, m)\")\n",
    "    println(\"Scale of internal units: \",scale,\" (SI,m)\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.246950256514266e18"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Polaron volume / touching\n",
    "\n",
    "# Polaron radius\n",
    "rp=26.8e-10 # SI units, metres\n",
    "# This is from Finite Temperature solution at 300 K for electrons (26.8)\n",
    "# Volume of a cube with 2*polaron radius to a side\n",
    "Vcube=(2*rp)^3\n",
    "# Double count as each excitation generates a (hole) and (electron) polaron\n",
    "Vcube=2*Vcube\n",
    "\n",
    "Vcube^-1 # = density per metre cubed, close packed\n",
    "Vcube^-1 / (100^3) # = density per centimetre cubed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Excitons with ε_Inf:4.5\n",
      "Effective masses: me: 1.0931258256e-31 mh: 1.366407282e-31 mu: 6.072921253333333e-32\n",
      "Exciton radius a0: 3.5719462619957756e-9 m 35.71946261995775 Å 5.761203648380283 Lattice Spacings\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# EXCITONS\n",
    "# Following p.973 Jones & March, Theoretical Solid State Phys Vol 2.\n",
    "# ISBN: 0486650162\n",
    "# Hygronic like Wannier exciton; Eqn (9.8.10); p.973\n",
    "\n",
    "ε=ε_Inf * ε_0 * 4*pi \n",
    "# term of 4pi appears because in the derivation of Dresselhaus 1956, \n",
    "# the Epsilon is used as a general prefactor in the denominator of a Coulomb term in the Hamiltonian\n",
    "\n",
    "me=m_eff_e * MassElectron\n",
    "mh=m_eff_h * MassElectron\n",
    "\n",
    "println(\"Excitons with ε_Inf:$(ε_Inf)\")\n",
    "\n",
    "# reduced mass - consider using polaron masses here\n",
    "mu=1/(1/me + 1/mh)\n",
    "\n",
    "\n",
    "println(\"Effective masses: me: $me mh: $mh mu: $mu\")\n",
    "\n",
    "En(k,n) = -mu*q^4/(2*ħ^2*ε^2*n^2) + ħ^2*k^2/(2*(me+mh))\n",
    "a0 = (ε*ħ^2)/(mu*q^2)\n",
    "d=6.2Å # lattice spacing\n",
    "println(\"Exciton radius a0: \",a0,\" m \",a0/Å,\" Å \",a0/d,\" Lattice Spacings\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exciton radius a0: 3.5719462619957756e-9 m 35.71946261995775 Å 5.761203648380283 Lattice Spacings\n",
      "Volume exciton: 1.9089887055693506e-25 m (SI), 190898.87055693503 Å^3  800.9922063581914 unit cells^3\n"
     ]
    }
   ],
   "source": [
    "# From above, exciton radius (SI) is a0\n",
    "V=(4/3)*π * a0^3\n",
    "println(\"Exciton radius a0: \",a0,\" m \",a0/Å,\" Å \",a0/d,\" Lattice Spacings\")\n",
    "println(\"Volume exciton: $V m (SI), \",V/Å^3, \" Å^3  \",V/d^3,\" unit cells^3\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Cross-checked exciton binding energy (figuring out the old-skool units) against this classic paper,\n",
    "# From http:/dx.doi.org/10.1103/PhysRev.116.573 - Thomas and Hopfield, PhysRev 1959 CdS\n",
    "# Eb=28 meV; Dielectric = 9.3; Mu = 0.18\n",
    "# mu=0.18 * melectron\n",
    "# ε_Inf=9.3\n",
    "# ε=ε_Inf * ε_0 * 4*pi \n",
    "# ---> n: 1 k: 0 En(k,n): -28.34762730330596 meV\n",
    "#   == Happy Jarvist 8^)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Hydrogenic exciton states @ the gamma point\n",
      "n: 1 k: 0.0 En(k,n): -44.79240010358902 meV\n",
      "n: 2 k: 0.0 En(k,n): -11.198100025897254 meV\n",
      "n: 3 k: 0.0 En(k,n): -4.976933344843224 meV\n",
      "n: 4 k: 0.0 En(k,n): -2.7995250064743136 meV\n",
      "n: 5 k: 0.0 En(k,n): -1.7916960041435608 meV\n",
      "n: 6 k: 0.0 En(k,n): -1.244233336210806 meV\n",
      "n: 7 k: 0.0 En(k,n): -0.9141306143589597 meV\n",
      "n: 8 k: 0.0 En(k,n): -0.6998812516185784 meV\n",
      "n: 9 k: 0.0 En(k,n): -0.5529925938714694 meV\n",
      "n: 10 k: 0.0 En(k,n): -0.4479240010358902 meV\n"
     ]
    }
   ],
   "source": [
    "println(\" Hydrogenic exciton states @ the gamma point\")\n",
    "k=0.0\n",
    "\n",
    "for n in 1:10\n",
    "   println(\"n: \",n,\" k: \",k,\" En(k,n): \",1000*En(k,n)/eV,\" meV\")\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Hydrogenic exciton states.\n",
      " n is excitation level; k is crystal momentum.\n"
     ]
    }
   ],
   "source": [
    "println(\" Hydrogenic exciton states.\\n n is excitation level; k is crystal momentum.\")\n",
    "\n",
    "# Very messy + not currently used anywhere\n",
    "#kmax=π/d\n",
    "#for k in -kmax:kmax/50:kmax\n",
    "#    for n in 1\n",
    "#        println(\"n: \",n,\" k: \",k,\" En(k,n): \",1000*En(k,n)/eV,\" meV\")\n",
    "#    end\n",
    "#end\n",
    "\n",
    "#krange=-kmax:kmax/50:kmax\n",
    "#n=1\n",
    "#E=[1000*En(k,n)/eV for k in krange]\n",
    "\n",
    "#using Plots\n",
    "#gr()\n",
    "#plot(krange,E,label=\"Exciton dispersion\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# OK, time to go fetch the phonons - these data are frequencies (not reduced) in THz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Array{Float64,1}:\n",
       " 0.110239\n",
       " 0.255401\n",
       " 0.725114\n",
       " 2.5411  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# via \n",
    "# jarvist@titanium:~/phonopy-work/2017-03-Scott-Distortions/0010-CsPbI3\n",
    "# > grep frequency modulation.yaml  | uniq | sed s/frequency://\n",
    "CsPbI3=[      \n",
    "    #-0.7638364901 # this mode is soft; but including it will mess the summation!\n",
    "      0.110239\n",
    "      0.2554009591\n",
    "      0.7251139061\n",
    "      2.5411013872]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15-element Array{Float64,1}:\n",
       " 0.560101\n",
       " 0.79995 \n",
       " 0.862435\n",
       " 0.978043\n",
       " 0.983314\n",
       " 1.03744 \n",
       " 1.7466  \n",
       " 1.97199 \n",
       " 2.0681  \n",
       " 2.2484  \n",
       " 2.42856 \n",
       " 2.74867 \n",
       " 3.51382 \n",
       " 3.83542 \n",
       " 3.98613 "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# OK, real data here, frequencies - For MAPI\n",
    "MAPI=[ #    -0.0000005420\n",
    "    #  0.0000001842\n",
    "    #  0.0000004295 # acoustic modes\n",
    "      0.5601007762\n",
    "      0.7999497115\n",
    "      0.8624353643\n",
    "      0.9780432985\n",
    "      0.9833140605\n",
    "      1.0374436221\n",
    "      1.7466045486\n",
    "      1.9719855248\n",
    "      2.0680970668\n",
    "      2.2483954169\n",
    "      2.4285633345\n",
    "      2.7486672126\n",
    "      3.5138222502\n",
    "      3.8354218617\n",
    "    3.9861276966]\n",
    "\n",
    "PhononFreqs=MAPI # for work further down below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MAPI Phonon energies, in meV: [2.31639, 3.30833, 3.56675, 4.04486, 4.06666, 4.29052, 7.22338, 8.15548, 8.55296, 9.29862, 10.0437, 11.3676, 14.532, 15.862, 16.4853]\n",
      "Sum of (cage mode, 0-15) MAPI phonon energies, in meV: 123.11456400472724\n"
     ]
    }
   ],
   "source": [
    "PhononEnergy = PhononFreqs.*1E12* 2*π * ħ # simply ħω ; including unit conversions\n",
    "println(\"MAPI Phonon energies, in meV: \",1000*PhononEnergy/q)\n",
    "println(\"Sum of (cage mode, 0-15) MAPI phonon energies, in meV: \",sum(1000*PhononEnergy/q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11.08608935565894"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T=300\n",
    "Z=sum(exp.(-PhononEnergy/(k_b*T)))\n",
    "# Now, we've just calculated the partition function (Z) for a single unit cell, at a certain temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15-element Array{Float64,1}:\n",
       " 0.0824723\n",
       " 0.0793678\n",
       " 0.0785784\n",
       " 0.0771385\n",
       " 0.0770735\n",
       " 0.0764089\n",
       " 0.0682141\n",
       " 0.0657985\n",
       " 0.0647945\n",
       " 0.0629523\n",
       " 0.0611638\n",
       " 0.0581106\n",
       " 0.0514156\n",
       " 0.0488373\n",
       " 0.0476739"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Boltzmann / Thermodynamic occupation, (assuming fermi statistics?)\n",
    "ThermalOccupation=exp.(-PhononEnergy/(k_b*T))/Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(ThermalOccupation) # Sanity check - partition function should normalise this to 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(Average Thermal) Energy in phonons (Boltzmann stats) of unitcell at T= 300, [0.191038, 0.262574, 0.280269, 0.312014, 0.313432, 0.327834, 0.492736, 0.536618, 0.554185, 0.58537, 0.614313, 0.660576, 0.747172, 0.774659, 0.785919] meV\n",
      "Sum: 7.438708995926876 meV\n"
     ]
    }
   ],
   "source": [
    "ThermalOccupation.*PhononEnergy\n",
    "println(\"(Average Thermal) Energy in phonons (Boltzmann stats) of unitcell at T= $T, \", 1000/q*ThermalOccupation.*PhononEnergy, \" meV\")\n",
    "println(\"Sum: \",sum(1000/q*ThermalOccupation.*PhononEnergy),\" meV\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BE thermal occupation of phonon states...\n",
      "Test for k_B T=25 meV; i.e. thermal energy phonon occupation at STP.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.6134392648438429"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "println(\"BE thermal occupation of phonon states...\")\n",
    "BE(T,En)=1./(exp.(En/(k_b*T)) - 1.0) # Bose-Einstein Occupation\n",
    "\n",
    "println(\"Test for k_B T=25 meV; i.e. thermal energy phonon occupation at STP.\")\n",
    "BE(300,25E-3*q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15-element Array{Float64,1}:\n",
       " 10.6679 \n",
       "  7.32489\n",
       "  6.75957\n",
       "  5.90436\n",
       "  5.87017\n",
       "  5.53921\n",
       "  3.10219\n",
       "  2.69614\n",
       "  2.5501 \n",
       "  2.31011\n",
       "  2.10624\n",
       "  1.81072\n",
       "  1.32557\n",
       "  1.18062\n",
       "  1.12097"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T=300\n",
    "\n",
    "BEThermalOccupation=BE(T,PhononEnergy)\n",
    "\n",
    "# This is quite interesting to know in general - average phonon occupation statistics at 300 K for MAPI\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(Average Thermal) Energy in phonons (BE stats) of unitcell at T= 300, [24.7111, 24.2331, 24.1097, 23.8823, 23.872, 23.7661, 22.4083, 21.9883, 21.8109, 21.4808, 21.1545, 20.5834, 19.2632, 18.727, 18.4795] meV\n",
      "Sum: 330.4704107906226 meV\n"
     ]
    }
   ],
   "source": [
    "println(\"(Average Thermal) Energy in phonons (BE stats) of unitcell at T= $T, \", 1000/q*BEThermalOccupation.*PhononEnergy, \" meV\")\n",
    "println(\"Sum: \",sum(1000/q*BEThermalOccupation.*PhononEnergy),\" meV\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Temperature 0  BE Phonon E: 0.0 meV SHC: NaN meV/K\n",
      "Temperature 20  BE Phonon E: 3.4786027554348755 meV SHC: 0.17393013777174376 meV/K\n",
      "Temperature 40  BE Phonon E: 16.62831330625127 meV SHC: 0.4157078326562817 meV/K\n",
      "Temperature 60  BE Phonon E: 35.38498771320442 meV SHC: 0.5897497952200736 meV/K\n",
      "Temperature 80  BE Phonon E: 56.96727972004571 meV SHC: 0.7120909965005714 meV/K\n",
      "Temperature 100  BE Phonon E: 80.03317450812287 meV SHC: 0.8003317450812287 meV/K\n",
      "Temperature 120  BE Phonon E: 103.94352106046085 meV SHC: 0.8661960088371737 meV/K\n",
      "Temperature 140  BE Phonon E: 128.37194699946122 meV SHC: 0.9169424785675802 meV/K\n",
      "Temperature 160  BE Phonon E: 153.13839696345636 meV SHC: 0.9571149810216022 meV/K\n",
      "Temperature 180  BE Phonon E: 178.13658203554803 meV SHC: 0.9896476779752669 meV/K\n",
      "Temperature 200  BE Phonon E: 203.30012127489843 meV SHC: 1.0165006063744921 meV/K\n",
      "Temperature 220  BE Phonon E: 228.58557888287612 meV SHC: 1.0390253585585278 meV/K\n",
      "Temperature 240  BE Phonon E: 253.9634080419086 meV SHC: 1.0581808668412858 meV/K\n",
      "Temperature 260  BE Phonon E: 279.41284289558683 meV SHC: 1.0746647803676417 meV/K\n",
      "Temperature 280  BE Phonon E: 304.9188783984365 meV SHC: 1.0889959942801304 meV/K\n",
      "Temperature 300  BE Phonon E: 330.4704107906226 meV SHC: 1.101568035968742 meV/K\n",
      "Temperature 320  BE Phonon E: 356.05905184209195 meV SHC: 1.1126845370065372 meV/K\n",
      "Temperature 340  BE Phonon E: 381.678349302363 meV SHC: 1.1225833803010676 meV/K\n",
      "Temperature 360  BE Phonon E: 407.32326049473465 meV SHC: 1.131453501374263 meV/K\n",
      "Temperature 380  BE Phonon E: 432.98978834491635 meV SHC: 1.1394468114339904 meV/K\n",
      "Temperature 400  BE Phonon E: 458.6747243751311 meV SHC: 1.1466868109378279 meV/K\n",
      "Temperature 420  BE Phonon E: 484.3754637846359 meV SHC: 1.1532749137729426 meV/K\n",
      "Temperature 440  BE Phonon E: 510.08987012952383 meV SHC: 1.1592951593852814 meV/K\n",
      "Temperature 460  BE Phonon E: 535.8161747741378 meV SHC: 1.1648177712481258 meV/K\n",
      "Temperature 480  BE Phonon E: 561.5529011360147 meV SHC: 1.1699018773666974 meV/K\n",
      "Temperature 500  BE Phonon E: 587.2988068845252 meV SHC: 1.1745976137690504 meV/K\n",
      "Temperature 520  BE Phonon E: 613.0528393248921 meV SHC: 1.178947767932485 meV/K\n",
      "Temperature 540  BE Phonon E: 638.8141005917277 meV SHC: 1.182989075169866 meV/K\n",
      "Temperature 560  BE Phonon E: 664.5818202279285 meV SHC: 1.186753250407015 meV/K\n",
      "Temperature 580  BE Phonon E: 690.3553333852527 meV SHC: 1.1902678161814702 meV/K\n",
      "Temperature 600  BE Phonon E: 716.1340633478615 meV SHC: 1.1935567722464357 meV/K\n",
      "Temperature 620  BE Phonon E: 741.9175074116728 meV SHC: 1.196641140986569 meV/K\n",
      "Temperature 640  BE Phonon E: 767.705225391784 meV SHC: 1.1995394146746625 meV/K\n",
      "Temperature 660  BE Phonon E: 793.4968302049922 meV SHC: 1.2022679245530186 meV/K\n",
      "Temperature 680  BE Phonon E: 819.2919801034132 meV SHC: 1.2048411472109017 meV/K\n",
      "Temperature 700  BE Phonon E: 845.0903722313534 meV SHC: 1.2072719603305049 meV/K\n",
      "Temperature 720  BE Phonon E: 870.8917372498362 meV SHC: 1.2095718572914391 meV/K\n",
      "Temperature 740  BE Phonon E: 896.6958348280891 meV SHC: 1.2117511281460664 meV/K\n",
      "Temperature 760  BE Phonon E: 922.5024498432363 meV SHC: 1.2138190129516266 meV/K\n",
      "Temperature 780  BE Phonon E: 948.311389161808 meV SHC: 1.2157838322587282 meV/K\n",
      "Temperature 800  BE Phonon E: 974.1224789017882 meV SHC: 1.2176530986272354 meV/K\n",
      "Temperature 820  BE Phonon E: 999.9355620935821 meV SHC: 1.2194336123092464 meV/K\n",
      "Temperature 840  BE Phonon E: 1025.7504966737338 meV SHC: 1.221131543659207 meV/K\n",
      "Temperature 860  BE Phonon E: 1051.5671537574806 meV SHC: 1.2227525043691634 meV/K\n",
      "Temperature 880  BE Phonon E: 1077.3854161459703 meV SHC: 1.2243016092567844 meV/K\n",
      "Temperature 900  BE Phonon E: 1103.2051770318012 meV SHC: 1.2257835300353346 meV/K\n",
      "Temperature 920  BE Phonon E: 1129.0263388728201 meV SHC: 1.2272025422530655 meV/K\n",
      "Temperature 940  BE Phonon E: 1154.8488124092123 meV SHC: 1.228562566392779 meV/K\n",
      "Temperature 960  BE Phonon E: 1180.6725158030724 meV SHC: 1.2298672039615337 meV/K\n",
      "Temperature 980  BE Phonon E: 1206.4973738830097 meV SHC: 1.2311197692683773 meV/K\n",
      "Temperature 1000  BE Phonon E: 1232.323317479136 meV SHC: 1.232323317479136 meV/K\n",
      "Temperature 1020  BE Phonon E: 1258.1502828360797 meV SHC: 1.233480669447137 meV/K\n",
      "Temperature 1040  BE Phonon E: 1283.9782110935412 meV SHC: 1.2345944337437895 meV/K\n",
      "Temperature 1060  BE Phonon E: 1309.8070478254972 meV SHC: 1.235667026250469 meV/K\n",
      "Temperature 1080  BE Phonon E: 1335.6367426304914 meV SHC: 1.2367006876208253 meV/K\n",
      "Temperature 1100  BE Phonon E: 1361.4672487665027 meV SHC: 1.237697498878639 meV/K\n",
      "Temperature 1120  BE Phonon E: 1387.2985228248451 meV SHC: 1.238659395379326 meV/K\n",
      "Temperature 1140  BE Phonon E: 1413.1305244383236 meV SHC: 1.2395881793318628 meV/K\n",
      "Temperature 1160  BE Phonon E: 1438.9632160195135 meV SHC: 1.2404855310513048 meV/K\n",
      "Temperature 1180  BE Phonon E: 1464.7965625255968 meV SHC: 1.2413530190894888 meV/K\n",
      "Temperature 1200  BE Phonon E: 1490.6305312466739 meV SHC: 1.2421921093722281 meV/K\n",
      "Temperature 1220  BE Phonon E: 1516.4650916148537 meV SHC: 1.243004173454798 meV/K\n",
      "Temperature 1240  BE Phonon E: 1542.300215031786 meV SHC: 1.2437904959933757 meV/K\n",
      "Temperature 1260  BE Phonon E: 1568.1358747125817 meV SHC: 1.244552281517922 meV/K\n",
      "Temperature 1280  BE Phonon E: 1593.9720455443587 meV SHC: 1.2452906605815302 meV/K\n",
      "Temperature 1300  BE Phonon E: 1619.8087039578027 meV SHC: 1.2460066953521558 meV/K\n",
      "Temperature 1320  BE Phonon E: 1645.6458278103933 meV SHC: 1.2467013847048434 meV/K\n",
      "Temperature 1340  BE Phonon E: 1671.4833962800712 meV SHC: 1.2473756688657247 meV/K\n",
      "Temperature 1360  BE Phonon E: 1697.3213897682829 meV SHC: 1.2480304336531491 meV/K\n",
      "Temperature 1380  BE Phonon E: 1723.1597898114167 meV SHC: 1.248666514356099 meV/K\n",
      "Temperature 1400  BE Phonon E: 1748.9985789998445 meV SHC: 1.2492846992856033 meV/K\n",
      "Temperature 1420  BE Phonon E: 1774.8377409037978 meV SHC: 1.2498857330308435 meV/K\n",
      "Temperature 1440  BE Phonon E: 1800.6772600053782 meV SHC: 1.2504703194481794 meV/K\n",
      "Temperature 1460  BE Phonon E: 1826.5171216361948 meV SHC: 1.2510391244083525 meV/K\n",
      "Temperature 1480  BE Phonon E: 1852.357311920007 meV SHC: 1.251592778324329 meV/K\n",
      "Temperature 1500  BE Phonon E: 1878.1978177199887 meV SHC: 1.2521318784799924 meV/K\n",
      "Temperature 1520  BE Phonon E: 1904.038626590127 meV SHC: 1.2526569911777152 meV/K\n",
      "Temperature 1540  BE Phonon E: 1929.8797267304044 meV SHC: 1.2531686537210418 meV/K\n",
      "Temperature 1560  BE Phonon E: 1955.7211069454388 meV SHC: 1.2536673762470762 meV/K\n",
      "Temperature 1580  BE Phonon E: 1981.562756606253 meV SHC: 1.2541536434216791 meV/K\n",
      "Temperature 1600  BE Phonon E: 2007.404665614898 meV SHC: 1.2546279160093112 meV/K\n",
      "Temperature 1620  BE Phonon E: 2033.2468243717058 meV SHC: 1.2550906323282134 meV/K\n",
      "Temperature 1640  BE Phonon E: 2059.0892237449475 meV SHC: 1.2555422096005777 meV/K\n",
      "Temperature 1660  BE Phonon E: 2084.9318550426565 meV SHC: 1.2559830452064196 meV/K\n",
      "Temperature 1680  BE Phonon E: 2110.7747099864646 meV SHC: 1.256413517849086 meV/K\n",
      "Temperature 1700  BE Phonon E: 2136.617780687319 meV SHC: 1.2568339886395994 meV/K\n",
      "Temperature 1720  BE Phonon E: 2162.4610596228445 meV SHC: 1.257244802106305 meV/K\n",
      "Temperature 1740  BE Phonon E: 2188.3045396163016 meV SHC: 1.2576462871358054 meV/K\n",
      "Temperature 1760  BE Phonon E: 2214.148213816946 meV SHC: 1.2580387578505374 meV/K\n",
      "Temperature 1780  BE Phonon E: 2239.9920756817473 meV SHC: 1.258422514427948 meV/K\n",
      "Temperature 1800  BE Phonon E: 2265.8361189582897 meV SHC: 1.2587978438657166 meV/K\n",
      "Temperature 1820  BE Phonon E: 2291.6803376688326 meV SHC: 1.2591650206971607 meV/K\n",
      "Temperature 1840  BE Phonon E: 2317.524726095379 meV SHC: 1.259524307660532 meV/K\n",
      "Temperature 1860  BE Phonon E: 2343.369278765719 meV SHC: 1.2598759563256554 meV/K\n",
      "Temperature 1880  BE Phonon E: 2369.21399044039 meV SHC: 1.2602202076810585 meV/K\n",
      "Temperature 1900  BE Phonon E: 2395.0588561003988 meV SHC: 1.2605572926844204 meV/K\n",
      "Temperature 1920  BE Phonon E: 2420.903870935768 meV SHC: 1.2608874327790458 meV/K\n",
      "Temperature 1940  BE Phonon E: 2446.7490303347718 meV SHC: 1.2612108403787483 meV/K\n",
      "Temperature 1960  BE Phonon E: 2472.5943298738175 meV SHC: 1.2615277193233763 meV/K\n",
      "Temperature 1980  BE Phonon E: 2498.4397653079463 meV SHC: 1.2618382653070437 meV/K\n",
      "Temperature 2000  BE Phonon E: 2524.2853325619226 meV SHC: 1.2621426662809614 meV/K\n",
      "Temperature 2020  BE Phonon E: 2550.1310277218254 meV SHC: 1.2624411028325868 meV/K\n",
      "Temperature 2040  BE Phonon E: 2575.976847027144 meV SHC: 1.2627337485427175 meV/K\n",
      "Temperature 2060  BE Phonon E: 2601.822786863333 meV SHC: 1.2630207703220062 meV/K\n",
      "Temperature 2080  BE Phonon E: 2627.668843754827 meV SHC: 1.263302328728282 meV/K\n",
      "Temperature 2100  BE Phonon E: 2653.5150143583714 meV SHC: 1.2635785782658913 meV/K\n",
      "Temperature 2120  BE Phonon E: 2679.361295456846 meV SHC: 1.2638496676683237 meV/K\n",
      "Temperature 2140  BE Phonon E: 2705.2076839533256 meV SHC: 1.2641157401651055 meV/K\n",
      "Temperature 2160  BE Phonon E: 2731.054176865541 meV SHC: 1.264376933734047 meV/K\n",
      "Temperature 2180  BE Phonon E: 2756.9007713206197 meV SHC: 1.264633381339734 meV/K\n",
      "Temperature 2200  BE Phonon E: 2782.747464550124 meV SHC: 1.264885211159147 meV/K\n",
      "Temperature 2220  BE Phonon E: 2808.594253885324 meV SHC: 1.265132546795191 meV/K\n",
      "Temperature 2240  BE Phonon E: 2834.4411367527696 meV SHC: 1.265375507478915 meV/K\n",
      "Temperature 2260  BE Phonon E: 2860.288110670076 meV SHC: 1.2656142082610957 meV/K\n",
      "Temperature 2280  BE Phonon E: 2886.135173241912 meV SHC: 1.265848760193821 meV/K\n",
      "Temperature 2300  BE Phonon E: 2911.982322156261 meV SHC: 1.2660792705027222 meV/K\n",
      "Temperature 2320  BE Phonon E: 2937.829555180764 meV SHC: 1.2663058427503293 meV/K\n",
      "Temperature 2340  BE Phonon E: 2963.6768701593846 meV SHC: 1.26652857699119 meV/K\n",
      "Temperature 2360  BE Phonon E: 2989.5242650091222 meV SHC: 1.2667475699191195 meV/K\n",
      "Temperature 2380  BE Phonon E: 3015.3717377169883 meV SHC: 1.266962915007138 meV/K\n",
      "Temperature 2400  BE Phonon E: 3041.219286337062 meV SHC: 1.2671747026404425 meV/K\n",
      "Temperature 2420  BE Phonon E: 3067.0669089877315 meV SHC: 1.2673830202428642 meV/K\n",
      "Temperature 2440  BE Phonon E: 3092.9146038490553 meV SHC: 1.2675879523971538 meV/K\n",
      "Temperature 2460  BE Phonon E: 3118.7623691602466 meV SHC: 1.26778958095945 meV/K\n",
      "Temperature 2480  BE Phonon E: 3144.6102032173244 meV SHC: 1.267987985168276 meV/K\n",
      "Temperature 2500  BE Phonon E: 3170.458104370791 meV SHC: 1.2681832417483163 meV/K\n",
      "Temperature 2520  BE Phonon E: 3196.30607102352 meV SHC: 1.2683754250093333 meV/K\n",
      "Temperature 2540  BE Phonon E: 3222.154101628666 meV SHC: 1.2685646069404197 meV/K\n",
      "Temperature 2560  BE Phonon E: 3248.002194687692 meV SHC: 1.2687508572998796 meV/K\n",
      "Temperature 2580  BE Phonon E: 3273.8503487485336 meV SHC: 1.268934243700982 meV/K\n",
      "Temperature 2600  BE Phonon E: 3299.698562403791 meV SHC: 1.2691148316937657 meV/K\n",
      "Temperature 2620  BE Phonon E: 3325.5468342890085 meV SHC: 1.269292684843133 meV/K\n",
      "Temperature 2640  BE Phonon E: 3351.395163081087 meV SHC: 1.269467864803442 meV/K\n",
      "Temperature 2660  BE Phonon E: 3377.243547496696 meV SHC: 1.2696404313897354 meV/K\n",
      "Temperature 2680  BE Phonon E: 3403.0919862907986 meV SHC: 1.2698104426458203 meV/K\n",
      "Temperature 2700  BE Phonon E: 3428.940478255257 meV SHC: 1.2699779549093544 meV/K\n",
      "Temperature 2720  BE Phonon E: 3454.7890222174383 meV SHC: 1.2701430228740582 meV/K\n",
      "Temperature 2740  BE Phonon E: 3480.6376170389435 meV SHC: 1.2703056996492494 meV/K\n",
      "Temperature 2760  BE Phonon E: 3506.4862616143714 meV SHC: 1.2704660368168013 meV/K\n",
      "Temperature 2780  BE Phonon E: 3532.3349548701085 meV SHC: 1.2706240844856505 meV/K\n",
      "Temperature 2800  BE Phonon E: 3558.183695763212 meV SHC: 1.2707798913440043 meV/K\n",
      "Temperature 2820  BE Phonon E: 3584.0324832803385 meV SHC: 1.27093350470934 meV/K\n",
      "Temperature 2840  BE Phonon E: 3609.881316436651 meV SHC: 1.2710849705762857 meV/K\n",
      "Temperature 2860  BE Phonon E: 3635.730194274871 meV SHC: 1.2712343336625422 meV/K\n",
      "Temperature 2880  BE Phonon E: 3661.5791158643037 meV SHC: 1.2713816374528832 meV/K\n",
      "Temperature 2900  BE Phonon E: 3687.4280802999265 meV SHC: 1.271526924241354 meV/K\n",
      "Temperature 2920  BE Phonon E: 3713.27708670148 meV SHC: 1.2716702351717397 meV/K\n",
      "Temperature 2940  BE Phonon E: 3739.1261342127077 meV SHC: 1.2718116102764312 meV/K\n",
      "Temperature 2960  BE Phonon E: 3764.9752220004234 meV SHC: 1.2719510885136565 meV/K\n",
      "Temperature 2980  BE Phonon E: 3790.8243492538486 meV SHC: 1.272088707803305 meV/K\n",
      "Temperature 3000  BE Phonon E: 3816.6735151838034 meV SHC: 1.2722245050612677 meV/K\n"
     ]
    }
   ],
   "source": [
    "Ts=[]\n",
    "SHCs=[]\n",
    "\n",
    "for T in 0:20:3000\n",
    "    print(\"Temperature $T \")\n",
    "    BEThermalOccupation=BE(T,PhononEnergy)\n",
    "    SHC=sum(1000/q*BEThermalOccupation.*PhononEnergy)/T\n",
    "    println(\" BE Phonon E: \",sum(1000/q*BEThermalOccupation.*PhononEnergy),\" meV SHC: \",sum(1000/q*BEThermalOccupation.*PhononEnergy)/T,\" meV/K\")\n",
    "    \n",
    "    append!(Ts,T)\n",
    "    append!(SHCs,SHC)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: using Plots.h in module Main conflicts with an existing identifier.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 600 400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip00\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"600\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "0,400 600,400 600,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip01\">\n",
       "    <rect x=\"120\" y=\"0\" width=\"421\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "55.1181,372.441 592.126,372.441 592.126,7.87402 55.1181,7.87402 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip02\">\n",
       "    <rect x=\"55\" y=\"7\" width=\"538\" height=\"366\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  55.1181,366.972 55.1181,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  144.619,366.972 144.619,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  234.121,366.972 234.121,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  323.622,366.972 323.622,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  413.123,366.972 413.123,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  502.625,366.972 502.625,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  592.126,366.972 592.126,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  63.1732,264.206 584.071,264.206 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  63.1732,98.236 584.071,98.236 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,372.441 592.126,372.441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,372.441 55.1181,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  144.619,372.441 144.619,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  234.121,372.441 234.121,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  323.622,372.441 323.622,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  413.123,372.441 413.123,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  502.625,372.441 502.625,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  592.126,372.441 592.126,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,372.441 55.1181,7.87402 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,264.206 63.1732,264.206 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,98.236 63.1732,98.236 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 55.1181, 386.241)\" x=\"55.1181\" y=\"386.241\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 144.619, 386.241)\" x=\"144.619\" y=\"386.241\">500</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 234.121, 386.241)\" x=\"234.121\" y=\"386.241\">1000</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 323.622, 386.241)\" x=\"323.622\" y=\"386.241\">1500</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 413.123, 386.241)\" x=\"413.123\" y=\"386.241\">2000</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 502.625, 386.241)\" x=\"502.625\" y=\"386.241\">2500</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 592.126, 386.241)\" x=\"592.126\" y=\"386.241\">3000</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 49.1181, 268.706)\" x=\"49.1181\" y=\"268.706\">0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 49.1181, 102.736)\" x=\"49.1181\" y=\"102.736\">1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(0, 323.622, 397.6)\" x=\"323.622\" y=\"397.6\">Temperature (K)</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(-90, 14.4, 190.157)\" x=\"14.4\" y=\"190.157\">Specific Heat Capacity (meV/K/unitcell)</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  58.6982,372.441 62.2782,292.185 65.8583,234.414 69.4383,193.804 73.0184,164.514 76.5984,142.651 80.1785,125.806 83.7585,112.471 87.3386,101.672 90.9186,92.7588 \n",
       "  94.4987,85.2819 98.0787,78.9235 101.659,73.4518 105.239,68.6947 108.819,64.5216 112.399,60.8316 115.979,57.5458 119.559,54.6014 123.139,51.9481 126.719,49.5449 \n",
       "  130.299,47.358 133.879,45.3597 137.459,43.5265 141.039,41.8389 144.619,40.2802 148.199,38.8362 151.78,37.4948 155.36,36.2453 158.94,35.0787 162.52,33.9869 \n",
       "  166.1,32.9631 169.68,32.001 173.26,31.0953 176.84,30.2412 180.42,29.4343 184,28.6709 187.58,27.9475 191.16,27.2611 194.74,26.6089 198.32,25.9884 \n",
       "  201.9,25.3974 205.48,24.8338 209.06,24.2957 212.64,23.7815 216.22,23.2896 219.801,22.8186 223.381,22.3671 226.961,21.9341 230.541,21.5183 234.121,21.1188 \n",
       "  237.701,20.7346 241.281,20.3649 244.861,20.0089 248.441,19.6658 252.021,19.3349 255.601,19.0156 259.181,18.7073 262.761,18.4094 266.341,18.1215 269.921,17.8429 \n",
       "  273.501,17.5734 277.081,17.3124 280.661,17.0595 284.241,16.8144 287.822,16.5767 291.402,16.3461 294.982,16.1223 298.562,15.905 302.142,15.6938 305.722,15.4886 \n",
       "  309.302,15.2891 312.882,15.0951 316.462,14.9063 320.042,14.7225 323.622,14.5435 327.202,14.3692 330.782,14.1994 334.362,14.0339 337.942,13.8724 341.522,13.715 \n",
       "  345.102,13.5614 348.682,13.4115 352.262,13.2652 355.843,13.1223 359.423,12.9827 363.003,12.8464 366.583,12.7131 370.163,12.5828 373.743,12.4554 377.323,12.3309 \n",
       "  380.903,12.209 384.483,12.0897 388.063,11.973 391.643,11.8587 395.223,11.7468 398.803,11.6372 402.383,11.5299 405.963,11.4247 409.543,11.3216 413.123,11.2206 \n",
       "  416.703,11.1215 420.283,11.0244 423.864,10.9291 427.444,10.8356 431.024,10.7439 434.604,10.654 438.184,10.5656 441.764,10.4789 445.344,10.3938 448.924,10.3102 \n",
       "  452.504,10.2281 456.084,10.1475 459.664,10.0682 463.244,9.99038 466.824,9.91386 470.404,9.83865 473.984,9.76472 477.564,9.69203 481.144,9.62054 484.724,9.55024 \n",
       "  488.304,9.48109 491.885,9.41307 495.465,9.34614 499.045,9.28028 502.625,9.21547 506.205,9.15168 509.785,9.08888 513.365,9.02706 516.945,8.96618 520.525,8.90624 \n",
       "  524.105,8.8472 527.685,8.78905 531.265,8.73177 534.845,8.67534 538.425,8.61973 542.005,8.56494 545.585,8.51094 549.165,8.45772 552.745,8.40526 556.325,8.35354 \n",
       "  559.906,8.30255 563.486,8.25227 567.066,8.20269 570.646,8.1538 574.226,8.10557 577.806,8.058 581.386,8.01107 584.966,7.96477 588.546,7.91909 592.126,7.87402 \n",
       "  \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "470.653,58.994 574.126,58.994 574.126,28.754 470.653,28.754 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  470.653,58.994 574.126,58.994 574.126,28.754 470.653,28.754 470.653,58.994 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  476.653,43.874 512.653,43.874 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 518.653, 48.374)\" x=\"518.653\" y=\"48.374\">BE Sum</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# OK, that's more or less linear in temperature - as the phonon modes become populated, the specific heat capacity --> a constant\n",
    "# So it may be just as easy to calculate specific heat capcity in phonopy and use those figures\n",
    "using Plots\n",
    "plot(Ts,SHCs,lab=\"BE Sum\")\n",
    "plot!(ylab=\"Specific Heat Capacity (meV/K/unitcell)\")\n",
    "plot!(xlab=\"Temperature (K)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "151-element Array{Float64,1}:\n",
       " NaN     \n",
       "  16.7813\n",
       "  40.1088\n",
       "  56.9009\n",
       "  68.7047\n",
       "  77.2185\n",
       "  83.5732\n",
       "  88.4694\n",
       "  92.3454\n",
       "  95.4842\n",
       "  98.0751\n",
       " 100.248 \n",
       " 102.097 \n",
       "   ⋮     \n",
       " 122.594 \n",
       " 122.609 \n",
       " 122.624 \n",
       " 122.638 \n",
       " 122.653 \n",
       " 122.667 \n",
       " 122.681 \n",
       " 122.695 \n",
       " 122.708 \n",
       " 122.722 \n",
       " 122.735 \n",
       " 122.748 "
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# meV / K / unit cell --> J / K / mol\n",
    "# OK; so factor of /1000 for meV -> eV\n",
    "# factor of q for eV -> J\n",
    "# factor of Na + 12(atoms) for moles?\n",
    "Na=6.022E23\n",
    "\n",
    "#1.1*q/1000*Na*12\n",
    "\n",
    "SHCsSI=SHCs.*q/1000*Na*1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 600 400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip00\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"600\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "0,400 600,400 600,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip01\">\n",
       "    <rect x=\"120\" y=\"0\" width=\"421\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "55.1181,372.441 592.126,372.441 592.126,7.87402 55.1181,7.87402 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip02\">\n",
       "    <rect x=\"55\" y=\"7\" width=\"538\" height=\"366\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  55.1181,366.972 55.1181,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  144.619,366.972 144.619,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  234.121,366.972 234.121,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  323.622,366.972 323.622,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  413.123,366.972 413.123,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  502.625,366.972 502.625,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  592.126,366.972 592.126,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  63.1732,372.441 584.071,372.441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  63.1732,298.19 584.071,298.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  63.1732,223.939 584.071,223.939 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  63.1732,149.688 584.071,149.688 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  63.1732,75.4369 584.071,75.4369 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,372.441 592.126,372.441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,372.441 55.1181,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  144.619,372.441 144.619,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  234.121,372.441 234.121,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  323.622,372.441 323.622,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  413.123,372.441 413.123,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  502.625,372.441 502.625,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  592.126,372.441 592.126,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,372.441 55.1181,7.87402 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,372.441 63.1732,372.441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,298.19 63.1732,298.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,223.939 63.1732,223.939 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,149.688 63.1732,149.688 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,75.4369 63.1732,75.4369 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 55.1181, 386.241)\" x=\"55.1181\" y=\"386.241\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 144.619, 386.241)\" x=\"144.619\" y=\"386.241\">500</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 234.121, 386.241)\" x=\"234.121\" y=\"386.241\">1000</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 323.622, 386.241)\" x=\"323.622\" y=\"386.241\">1500</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 413.123, 386.241)\" x=\"413.123\" y=\"386.241\">2000</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 502.625, 386.241)\" x=\"502.625\" y=\"386.241\">2500</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 592.126, 386.241)\" x=\"592.126\" y=\"386.241\">3000</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 49.1181, 376.941)\" x=\"49.1181\" y=\"376.941\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 49.1181, 302.69)\" x=\"49.1181\" y=\"302.69\">25</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 49.1181, 228.439)\" x=\"49.1181\" y=\"228.439\">50</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 49.1181, 154.188)\" x=\"49.1181\" y=\"154.188\">75</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 49.1181, 79.9369)\" x=\"49.1181\" y=\"79.9369\">100</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(0, 323.622, 397.6)\" x=\"323.622\" y=\"397.6\">Temperature (K)</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(-90, 14.4, 190.157)\" x=\"14.4\" y=\"190.157\">Specific Heat Capacity (J/K/mol (unitcell))</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  58.6982,322.6 62.2782,253.316 65.8583,203.443 69.4383,168.385 73.0184,143.099 76.5984,124.225 80.1785,109.683 83.7585,98.1714 87.3386,88.8488 90.9186,81.1539 \n",
       "  94.4987,74.6992 98.0787,69.2101 101.659,64.4865 105.239,60.3797 108.819,56.7771 112.399,53.5916 115.979,50.755 119.559,48.2132 123.139,45.9226 126.719,43.8479 \n",
       "  130.299,41.96 133.879,40.2349 137.459,38.6523 141.039,37.1954 144.619,35.8498 148.199,34.6033 151.78,33.4452 155.36,32.3665 158.94,31.3594 162.52,30.4169 \n",
       "  166.1,29.5331 169.68,28.7026 173.26,27.9207 176.84,27.1833 180.42,26.4867 184,25.8277 187.58,25.2032 191.16,24.6106 194.74,24.0476 198.32,23.5119 \n",
       "  201.9,23.0017 205.48,22.5151 209.06,22.0506 212.64,21.6067 216.22,21.1821 219.801,20.7754 223.381,20.3857 226.961,20.0119 230.541,19.6529 234.121,19.308 \n",
       "  237.701,18.9764 241.281,18.6572 244.861,18.3499 248.441,18.0537 252.021,17.768 255.601,17.4924 259.181,17.2262 262.761,16.9691 266.341,16.7205 269.921,16.4801 \n",
       "  273.501,16.2474 277.081,16.022 280.661,15.8037 284.241,15.5921 287.822,15.387 291.402,15.1879 294.982,14.9947 298.562,14.807 302.142,14.6248 305.722,14.4476 \n",
       "  309.302,14.2754 312.882,14.1079 316.462,13.9449 320.042,13.7862 323.622,13.6317 327.202,13.4813 330.782,13.3346 334.362,13.1917 337.942,13.0524 341.522,12.9165 \n",
       "  345.102,12.7839 348.682,12.6545 352.262,12.5281 355.843,12.4048 359.423,12.2843 363.003,12.1666 366.583,12.0515 370.163,11.9391 373.743,11.8291 377.323,11.7215 \n",
       "  380.903,11.6163 384.483,11.5134 388.063,11.4126 391.643,11.314 395.223,11.2174 398.803,11.1228 402.383,11.0301 405.963,10.9393 409.543,10.8503 413.123,10.7631 \n",
       "  416.703,10.6775 420.283,10.5937 423.864,10.5114 427.444,10.4307 431.024,10.3516 434.604,10.2739 438.184,10.1977 441.764,10.1228 445.344,10.0493 448.924,9.97715 \n",
       "  452.504,9.90628 456.084,9.83666 459.664,9.76825 463.244,9.70104 466.824,9.63499 470.404,9.57006 473.984,9.50623 477.564,9.44348 481.144,9.38177 484.724,9.32108 \n",
       "  488.304,9.26139 491.885,9.20266 495.465,9.14488 499.045,9.08803 502.625,9.03207 506.205,8.977 509.785,8.92279 513.365,8.86942 516.945,8.81687 520.525,8.76512 \n",
       "  524.105,8.71415 527.685,8.66395 531.265,8.6145 534.845,8.56579 538.425,8.51778 542.005,8.47048 545.585,8.42387 549.165,8.37792 552.745,8.33263 556.325,8.28798 \n",
       "  559.906,8.24396 563.486,8.20056 567.066,8.15776 570.646,8.11555 574.226,8.07391 577.806,8.03285 581.386,7.99233 584.966,7.95237 588.546,7.91293 592.126,7.87402 \n",
       "  \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "470.653,58.994 574.126,58.994 574.126,28.754 470.653,28.754 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  470.653,58.994 574.126,58.994 574.126,28.754 470.653,28.754 470.653,58.994 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  476.653,43.874 512.653,43.874 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 518.653, 48.374)\" x=\"518.653\" y=\"48.374\">BE Sum</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Plots\n",
    "plot(Ts,SHCsSI,lab=\"BE Sum\")\n",
    "plot!(ylim=[0,Inf])\n",
    "plot!(ylab=\"Specific Heat Capacity (J/K/mol (unitcell))\")\n",
    "plot!(xlab=\"Temperature (K)\")\n",
    "\n",
    "#savefig(\"SHC_vs_T.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Characteristic phonon energy: 2.25e12 Hz = 0.009305251500657504 eV\n"
     ]
    }
   ],
   "source": [
    "#Characteristic energy\n",
    "f=2.25E12 # 2.25 THz - see Polaron paper, http://arxiv.org/abs/1704.05404 and references therein\n",
    "ω=2*π*f\n",
    "E=ħ*ω\n",
    "println(\"Characteristic phonon energy: $f Hz = \",E/q,\" eV\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial energy loss to scattering into LO modes, 0.0775437625054792 eV/ps\n"
     ]
    }
   ],
   "source": [
    "# optical phonon scattering rate\n",
    "τ=0.12 # in ps\n",
    "println(\"Initial energy loss to scattering into LO modes, \",(E/q)/τ,\" eV/ps\") # in eV/ps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.888888888888889"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For paper --- calculation of polaron temperature, assuming modal occupation as above\n",
    "cells=360 # from volume Polaron, r=27.4 Angstrom\n",
    "SHC=1.25E-3 # meV/K, from rough of above table\n",
    "E=1.75 # maximum of energy in eV\n",
    "\n",
    "T=E/(SHC*cells)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.176679038355128e-7"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Should really use unitful here!\n",
    "\n",
    "ω=2.25E12\n",
    "meff=0.12\n",
    "\n",
    "shm=sqrt((ħ*ω)^2/(meff*MassElectron))\n",
    "#sqrt(ω^2*(meff*me))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "From finite Temperature T=300 variational solution for MAPI (Polaron paper)\n",
      "Franck condon overtone energy: (v=20, f=2.25e12 Hz) = 186.10503001315007 meV\n",
      "Franck condon overtone energy: (v=20, f=1.0e12 Hz) = 82.71334667251114 meV\n",
      "From Feynman small-α approximation, athermal variational solution, for α=2.6765070822960193. \n",
      "Franck condon overtone energy: (v=3.5947793516213373, f=2.25e12 Hz) = 33.45032595620706 meV\n",
      "Franck condon overtone energy: (v=3.5947793516213373, f=1.0e12 Hz) = 14.866811536092026 meV\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "6.626068959988852e-22"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Overtones in MAPI Polarons\n",
    "# These are the Franck-Condon vibration states\n",
    "\n",
    "# Schultz: \n",
    "#phononfreq= 17.1009 #(int units) = 3.8477e+13 [SI, Hz]\n",
    "#electronfreq= 10.3716 #(int units) = 2.33361e+13 [SI, Hz]\n",
    "#combinedfreq= 8.86807 #(int units) = 1.99532e+13 [SI, Hz]\n",
    "#println(\"Energy: (v=$v, f=$f Hz) = \",1000*v*E/q,\" meV\")\n",
    "\n",
    "#Characteristic energy\n",
    "v=20\n",
    "f=2.25E12 # 2.25 THz\n",
    "\n",
    "function FC_overtone(;v=20,f=2.25E12)\n",
    "    ω=2*π*f\n",
    "    E=ħ*ω\n",
    "    println(\"Franck condon overtone energy: (v=$v, f=$f Hz) = \",1000*v*E/q,\" meV\")\n",
    "    return(E)\n",
    "end\n",
    "\n",
    "println(\"From finite Temperature T=300 variational solution for MAPI (Polaron paper)\")\n",
    "FC_overtone(v=20,f=2.25E12)\n",
    "FC_overtone(v=20,f=1.0E12)\n",
    "\n",
    "println(\"From Feynman small-α approximation, athermal variational solution, for α=$(α). \")\n",
    "w=3 ; v=w+(2/9)*α\n",
    "FC_overtone(v=v,f=2.25E12)\n",
    "FC_overtone(v=v,f=1.0E12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "MAPI(meff=0.12): α=2.3939407113094093\n",
      "v=3.5319868247354242 w=3 by the Feynman1955 small-α expansion\n",
      "\n",
      " Schultz1959(2.4): rf= 1.23476 (int units) = 8.66011e-10 m [SI]"
     ]
    }
   ],
   "source": [
    "α=feynmanalpha(ε_Inf, ε_S ,  f,    meff)\n",
    "println(\"\\nMAPI(meff=$meff): α=$α\")\n",
    "\n",
    "scale=sqrt(2*meff*melectron*ω )\n",
    "w=3 ; v=w+(2/9)*α\n",
    "println(\"v=$v w=$w by the Feynman1955 small-α expansion\")\n",
    "\n",
    "\n",
    "#v=20.0; w=1.15\n",
    "\n",
    "mu=μ=((v^2-w^2)/v^2)\n",
    "\n",
    "# (2.4)\n",
    "rf=sqrt(3/(2*mu*v))\n",
    "# (2.4) SI scaling inferred from units in (2.5a) and Table II\n",
    "rfsi=rf*sqrt(2*me*ω )\n",
    "@printf(\"\\n Schultz1959(2.4): rf= %g (int units) = %g m [SI]\",rf,rfsi )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 600 400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip00\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"600\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "0,400 600,400 600,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip01\">\n",
       "    <rect x=\"120\" y=\"0\" width=\"421\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "55.1181,372.441 592.126,372.441 592.126,7.87402 55.1181,7.87402 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip02\">\n",
       "    <rect x=\"55\" y=\"7\" width=\"538\" height=\"366\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  178.652,366.972 178.652,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  323.622,366.972 323.622,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  468.592,366.972 468.592,13.3425 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  63.1732,372.441 584.071,372.441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  63.1732,285.356 584.071,285.356 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  63.1732,198.271 584.071,198.271 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  63.1732,111.185 584.071,111.185 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  63.1732,24.1001 584.071,24.1001 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,372.441 592.126,372.441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  178.652,372.441 178.652,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  323.622,372.441 323.622,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  468.592,372.441 468.592,366.972 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,372.441 55.1181,7.87402 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,372.441 63.1732,372.441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,285.356 63.1732,285.356 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,198.271 63.1732,198.271 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,111.185 63.1732,111.185 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,24.1001 63.1732,24.1001 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 178.652, 386.241)\" x=\"178.652\" y=\"386.241\">-2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 323.622, 386.241)\" x=\"323.622\" y=\"386.241\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 468.592, 386.241)\" x=\"468.592\" y=\"386.241\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 49.1181, 376.941)\" x=\"49.1181\" y=\"376.941\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 49.1181, 289.856)\" x=\"49.1181\" y=\"289.856\">0.1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 49.1181, 202.771)\" x=\"49.1181\" y=\"202.771\">0.2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 49.1181, 115.685)\" x=\"49.1181\" y=\"115.685\">0.3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 49.1181, 28.6001)\" x=\"49.1181\" y=\"28.6001\">0.4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(0, 323.622, 397.6)\" x=\"323.622\" y=\"397.6\">Radius (internal units)</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(-90, 14.4, 190.157)\" x=\"14.4\" y=\"190.157\">Wavefunction</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  55.1181,372.014 56.0131,371.994 56.9081,371.974 57.8031,371.953 58.6982,371.931 59.5932,371.907 60.4882,371.883 61.3832,371.858 62.2782,371.832 63.1732,371.805 \n",
       "  64.0682,371.776 64.9633,371.747 65.8583,371.716 66.7533,371.684 67.6483,371.651 68.5433,371.617 69.4383,371.581 70.3333,371.543 71.2283,371.504 72.1234,371.464 \n",
       "  73.0184,371.422 73.9134,371.378 74.8084,371.333 75.7034,371.286 76.5984,371.237 77.4934,371.186 78.3885,371.134 79.2835,371.079 80.1785,371.022 81.0735,370.963 \n",
       "  81.9685,370.902 82.8635,370.838 83.7585,370.772 84.6535,370.704 85.5486,370.633 86.4436,370.56 87.3386,370.484 88.2336,370.405 89.1286,370.323 90.0236,370.238 \n",
       "  90.9186,370.151 91.8136,370.06 92.7087,369.966 93.6037,369.868 94.4987,369.767 95.3937,369.663 96.2887,369.555 97.1837,369.443 98.0787,369.327 98.9738,369.207 \n",
       "  99.8688,369.083 100.764,368.955 101.659,368.823 102.554,368.686 103.449,368.545 104.344,368.399 105.239,368.248 106.134,368.092 107.029,367.93 107.924,367.764 \n",
       "  108.819,367.592 109.714,367.415 110.609,367.232 111.504,367.043 112.399,366.848 113.294,366.647 114.189,366.439 115.084,366.226 115.979,366.005 116.874,365.778 \n",
       "  117.769,365.543 118.664,365.301 119.559,365.053 120.454,364.796 121.349,364.532 122.244,364.26 123.139,363.98 124.034,363.691 124.929,363.394 125.824,363.089 \n",
       "  126.719,362.774 127.614,362.45 128.509,362.118 129.404,361.775 130.299,361.423 131.194,361.061 132.089,360.689 132.984,360.307 133.879,359.914 134.774,359.51 \n",
       "  135.669,359.095 136.564,358.669 137.459,358.231 138.354,357.782 139.249,357.321 140.144,356.848 141.039,356.362 141.934,355.864 142.829,355.352 143.724,354.828 \n",
       "  144.619,354.29 145.514,353.739 146.409,353.174 147.304,352.594 148.199,352.001 149.094,351.392 149.99,350.769 150.885,350.131 151.78,349.477 152.675,348.808 \n",
       "  153.57,348.123 154.465,347.422 155.36,346.704 156.255,345.97 157.15,345.219 158.045,344.451 158.94,343.666 159.835,342.863 160.73,342.042 161.625,341.203 \n",
       "  162.52,340.346 163.415,339.47 164.31,338.575 165.205,337.661 166.1,336.728 166.995,335.775 167.89,334.803 168.785,333.811 169.68,332.798 170.575,331.765 \n",
       "  171.47,330.711 172.365,329.637 173.26,328.541 174.155,327.424 175.05,326.286 175.945,325.126 176.84,323.944 177.735,322.739 178.63,321.513 179.525,320.264 \n",
       "  180.42,318.993 181.315,317.699 182.21,316.382 183.105,315.041 184,313.678 184.895,312.291 185.79,310.881 186.685,309.447 187.58,307.99 188.475,306.508 \n",
       "  189.37,305.003 190.265,303.473 191.16,301.92 192.055,300.342 192.95,298.74 193.845,297.114 194.74,295.464 195.635,293.789 196.53,292.089 197.425,290.366 \n",
       "  198.32,288.618 199.215,286.845 200.11,285.048 201.005,283.227 201.9,281.382 202.795,279.512 203.69,277.618 204.585,275.7 205.48,273.758 206.375,271.792 \n",
       "  207.27,269.803 208.165,267.79 209.06,265.753 209.955,263.693 210.85,261.61 211.745,259.503 212.64,257.374 213.535,255.223 214.43,253.049 215.325,250.853 \n",
       "  216.22,248.636 217.115,246.396 218.01,244.136 218.906,241.855 219.801,239.552 220.696,237.23 221.591,234.888 222.486,232.526 223.381,230.145 224.276,227.745 \n",
       "  225.171,225.327 226.066,222.89 226.961,220.436 227.856,217.965 228.751,215.478 229.646,212.974 230.541,210.455 231.436,207.92 232.331,205.371 233.226,202.808 \n",
       "  234.121,200.232 235.016,197.642 235.911,195.04 236.806,192.427 237.701,189.802 238.596,187.167 239.491,184.522 240.386,181.868 241.281,179.205 242.176,176.535 \n",
       "  243.071,173.857 243.966,171.173 244.861,168.484 245.756,165.789 246.651,163.09 247.546,160.388 248.441,157.683 249.336,154.977 250.231,152.269 251.126,149.561 \n",
       "  252.021,146.853 252.916,144.147 253.811,141.443 254.706,138.742 255.601,136.045 256.496,133.352 257.391,130.665 258.286,127.985 259.181,125.312 260.076,122.647 \n",
       "  260.971,119.992 261.866,117.346 262.761,114.712 263.656,112.089 264.551,109.479 265.446,106.882 266.341,104.301 267.236,101.734 268.131,99.1847 269.026,96.6522 \n",
       "  269.921,94.138 270.816,91.6431 271.711,89.1682 272.606,86.7144 273.501,84.2826 274.396,81.8737 275.291,79.4885 276.186,77.1281 277.081,74.7934 277.976,72.4851 \n",
       "  278.871,70.2043 279.766,67.9519 280.661,65.7287 281.556,63.5355 282.451,61.3734 283.346,59.2431 284.241,57.1455 285.136,55.0815 286.031,53.0518 286.927,51.0574 \n",
       "  287.822,49.0991 288.717,47.1776 289.612,45.2938 290.507,43.4484 291.402,41.6422 292.297,39.876 293.192,38.1505 294.087,36.4665 294.982,34.8246 295.877,33.2256 \n",
       "  296.772,31.6701 297.667,30.1589 298.562,28.6925 299.457,27.2716 300.352,25.8968 301.247,24.5688 302.142,23.288 303.037,22.055 303.932,20.8705 304.827,19.7348 \n",
       "  305.722,18.6486 306.617,17.6122 307.512,16.6262 308.407,15.691 309.302,14.8069 310.197,13.9745 311.092,13.194 311.987,12.4657 312.882,11.7902 313.777,11.1675 \n",
       "  314.672,10.598 315.567,10.082 316.462,9.61974 317.357,9.21134 318.252,8.85702 319.147,8.55694 320.042,8.31123 320.937,8.12002 321.832,7.98337 322.727,7.90136 \n",
       "  323.622,7.87402 324.517,7.90136 325.412,7.98337 326.307,8.12002 327.202,8.31123 328.097,8.55694 328.992,8.85702 329.887,9.21134 330.782,9.61974 331.677,10.082 \n",
       "  332.572,10.598 333.467,11.1675 334.362,11.7902 335.257,12.4657 336.152,13.194 337.047,13.9745 337.942,14.8069 338.837,15.691 339.732,16.6262 340.627,17.6122 \n",
       "  341.522,18.6486 342.417,19.7348 343.312,20.8705 344.207,22.055 345.102,23.288 345.997,24.5688 346.892,25.8968 347.787,27.2716 348.682,28.6925 349.577,30.1589 \n",
       "  350.472,31.6701 351.367,33.2256 352.262,34.8246 353.157,36.4665 354.052,38.1505 354.948,39.876 355.843,41.6422 356.738,43.4484 357.633,45.2938 358.528,47.1776 \n",
       "  359.423,49.0991 360.318,51.0574 361.213,53.0518 362.108,55.0815 363.003,57.1455 363.898,59.2431 364.793,61.3734 365.688,63.5355 366.583,65.7287 367.478,67.9519 \n",
       "  368.373,70.2043 369.268,72.4851 370.163,74.7934 371.058,77.1281 371.953,79.4885 372.848,81.8737 373.743,84.2826 374.638,86.7144 375.533,89.1682 376.428,91.6431 \n",
       "  377.323,94.138 378.218,96.6522 379.113,99.1847 380.008,101.734 380.903,104.301 381.798,106.882 382.693,109.479 383.588,112.089 384.483,114.712 385.378,117.346 \n",
       "  386.273,119.992 387.168,122.647 388.063,125.312 388.958,127.985 389.853,130.665 390.748,133.352 391.643,136.045 392.538,138.742 393.433,141.443 394.328,144.147 \n",
       "  395.223,146.853 396.118,149.561 397.013,152.269 397.908,154.977 398.803,157.683 399.698,160.388 400.593,163.09 401.488,165.789 402.383,168.484 403.278,171.173 \n",
       "  404.173,173.857 405.068,176.535 405.963,179.205 406.858,181.868 407.753,184.522 408.648,187.167 409.543,189.802 410.438,192.427 411.333,195.04 412.228,197.642 \n",
       "  413.123,200.232 414.018,202.808 414.913,205.371 415.808,207.92 416.703,210.455 417.598,212.974 418.493,215.478 419.388,217.965 420.283,220.436 421.178,222.89 \n",
       "  422.073,225.327 422.969,227.745 423.864,230.145 424.759,232.526 425.654,234.888 426.549,237.23 427.444,239.552 428.339,241.855 429.234,244.136 430.129,246.396 \n",
       "  431.024,248.636 431.919,250.853 432.814,253.049 433.709,255.223 434.604,257.374 435.499,259.503 436.394,261.61 437.289,263.693 438.184,265.753 439.079,267.79 \n",
       "  439.974,269.803 440.869,271.792 441.764,273.758 442.659,275.7 443.554,277.618 444.449,279.512 445.344,281.382 446.239,283.227 447.134,285.048 448.029,286.845 \n",
       "  448.924,288.618 449.819,290.366 450.714,292.089 451.609,293.789 452.504,295.464 453.399,297.114 454.294,298.74 455.189,300.342 456.084,301.92 456.979,303.473 \n",
       "  457.874,305.003 458.769,306.508 459.664,307.99 460.559,309.447 461.454,310.881 462.349,312.291 463.244,313.678 464.139,315.041 465.034,316.382 465.929,317.699 \n",
       "  466.824,318.993 467.719,320.264 468.614,321.513 469.509,322.739 470.404,323.944 471.299,325.126 472.194,326.286 473.089,327.424 473.984,328.541 474.879,329.637 \n",
       "  475.774,330.711 476.669,331.765 477.564,332.798 478.459,333.811 479.354,334.803 480.249,335.775 481.144,336.728 482.039,337.661 482.934,338.575 483.829,339.47 \n",
       "  484.724,340.346 485.619,341.203 486.514,342.042 487.409,342.863 488.304,343.666 489.199,344.451 490.094,345.219 490.99,345.97 491.885,346.704 492.78,347.422 \n",
       "  493.675,348.123 494.57,348.808 495.465,349.477 496.36,350.131 497.255,350.769 498.15,351.392 499.045,352.001 499.94,352.594 500.835,353.174 501.73,353.739 \n",
       "  502.625,354.29 503.52,354.828 504.415,355.352 505.31,355.864 506.205,356.362 507.1,356.848 507.995,357.321 508.89,357.782 509.785,358.231 510.68,358.669 \n",
       "  511.575,359.095 512.47,359.51 513.365,359.914 514.26,360.307 515.155,360.689 516.05,361.061 516.945,361.423 517.84,361.775 518.735,362.118 519.63,362.45 \n",
       "  520.525,362.774 521.42,363.089 522.315,363.394 523.21,363.691 524.105,363.98 525,364.26 525.895,364.532 526.79,364.796 527.685,365.053 528.58,365.301 \n",
       "  529.475,365.543 530.37,365.778 531.265,366.005 532.16,366.226 533.055,366.439 533.95,366.647 534.845,366.848 535.74,367.043 536.635,367.232 537.53,367.415 \n",
       "  538.425,367.592 539.32,367.764 540.215,367.93 541.11,368.092 542.005,368.248 542.9,368.399 543.795,368.545 544.69,368.686 545.585,368.823 546.48,368.955 \n",
       "  547.375,369.083 548.27,369.207 549.165,369.327 550.06,369.443 550.955,369.555 551.85,369.663 552.745,369.767 553.64,369.868 554.535,369.966 555.43,370.06 \n",
       "  556.325,370.151 557.22,370.238 558.115,370.323 559.01,370.405 559.906,370.484 560.801,370.56 561.696,370.633 562.591,370.704 563.486,370.772 564.381,370.838 \n",
       "  565.276,370.902 566.171,370.963 567.066,371.022 567.961,371.079 568.856,371.134 569.751,371.186 570.646,371.237 571.541,371.286 572.436,371.333 573.331,371.378 \n",
       "  574.226,371.422 575.121,371.464 576.016,371.504 576.911,371.543 577.806,371.581 578.701,371.617 579.596,371.651 580.491,371.684 581.386,371.716 582.281,371.747 \n",
       "  583.176,371.776 584.071,371.805 584.966,371.832 585.861,371.858 586.756,371.883 587.651,371.907 588.546,371.931 589.441,371.953 590.336,371.974 591.231,371.994 \n",
       "  592.126,372.014 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip02)\" points=\"\n",
       "234.121,7.87402 413.123,7.87402 413.123,372.441 234.121,372.441 \n",
       "  \" fill=\"#ff0000\" fill-opacity=\"0.2\"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#e26f46; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  234.121,7.87402 413.123,7.87402 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "398.191,74.114 574.126,74.114 574.126,28.754 398.191,28.754 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  398.191,74.114 574.126,74.114 574.126,28.754 398.191,28.754 398.191,74.114 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  404.191,43.874 440.191,43.874 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 446.191, 48.374)\" x=\"446.191\" y=\"48.374\">Polaron Wavefunction</text>\n",
       "</g>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "404.191,65.042 440.191,65.042 440.191,52.946 404.191,52.946 404.191,65.042 \n",
       "  \" fill=\"#ff0000\" fill-opacity=\"0.2\"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#e26f46; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  404.191,52.946 440.191,52.946 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 446.191, 63.494)\" x=\"446.191\" y=\"63.494\">Polaron Radius</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's cross check this by visualising!\n",
    "\n",
    "# Schultz1959 - above (2.4) - SHM wave function\n",
    "ψ(r)=(μ*v/π)^(3/4) * exp(-μ*v*r^2/2)\n",
    "\n",
    "rs=-rf*3:rf/100:rf*3\n",
    "y=[ψ(r) for r in rs]\n",
    "\n",
    "using Plots\n",
    "gr()\n",
    "\n",
    "plot(rs,y,label=\"Polaron Wavefunction\") # fmt=:png)\n",
    "#plot!(x,y.^2,label=\"Polaron Wavefunction ^2\")\n",
    "plot!([0-rf,0+rf],[ψ(0),ψ(0)],label=\"Polaron Radius\",fill=(0,:red),fillalpha=0.2)\n",
    "plot!(xlab=\"Radius (internal units)\")\n",
    "plot!(ylab=\"Wavefunction\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "savefig(\"polaron-radius.png\")\n",
    "\n",
    "# needs to be in a cell by itself as otherwise it breaks the inline output above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 600 400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip00\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"600\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "0,400 600,400 600,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip01\">\n",
       "    <rect x=\"120\" y=\"0\" width=\"421\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "39.3701,376.378 592.126,376.378 592.126,7.87402 39.3701,7.87402 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip02\">\n",
       "    <rect x=\"39\" y=\"7\" width=\"554\" height=\"370\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  166.527,370.85 166.527,13.4016 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  315.748,370.85 315.748,13.4016 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  464.969,370.85 464.969,13.4016 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  47.6614,376.378 583.835,376.378 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  47.6614,271.243 583.835,271.243 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  47.6614,166.108 583.835,166.108 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:0.5; fill:none\" stroke-dasharray=\"1, 2\" points=\"\n",
       "  47.6614,60.9737 583.835,60.9737 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  39.3701,376.378 592.126,376.378 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  166.527,376.378 166.527,370.85 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  315.748,376.378 315.748,370.85 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  464.969,376.378 464.969,370.85 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  39.3701,376.378 39.3701,7.87402 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  39.3701,376.378 47.6614,376.378 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  39.3701,271.243 47.6614,271.243 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  39.3701,166.108 47.6614,166.108 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  39.3701,60.9737 47.6614,60.9737 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 166.527, 390.178)\" x=\"166.527\" y=\"390.178\">-2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 315.748, 390.178)\" x=\"315.748\" y=\"390.178\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 464.969, 390.178)\" x=\"464.969\" y=\"390.178\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 33.3701, 380.878)\" x=\"33.3701\" y=\"380.878\">0.00</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 33.3701, 275.743)\" x=\"33.3701\" y=\"275.743\">0.05</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 33.3701, 170.608)\" x=\"33.3701\" y=\"170.608\">0.10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 33.3701, 65.4737)\" x=\"33.3701\" y=\"65.4737\">0.15</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  39.3701,376.377 40.2913,376.377 41.2126,376.377 42.1339,376.377 43.0551,376.377 43.9764,376.377 44.8976,376.377 45.8189,376.377 46.7402,376.377 47.6614,376.377 \n",
       "  48.5827,376.377 49.5039,376.377 50.4252,376.376 51.3465,376.376 52.2677,376.376 53.189,376.376 54.1102,376.376 55.0315,376.376 55.9528,376.376 56.874,376.375 \n",
       "  57.7953,376.375 58.7165,376.375 59.6378,376.375 60.5591,376.374 61.4803,376.374 62.4016,376.374 63.3228,376.373 64.2441,376.373 65.1654,376.372 66.0866,376.372 \n",
       "  67.0079,376.371 67.9291,376.371 68.8504,376.37 69.7717,376.37 70.6929,376.369 71.6142,376.368 72.5354,376.367 73.4567,376.366 74.378,376.366 75.2992,376.365 \n",
       "  76.2205,376.363 77.1417,376.362 78.063,376.361 78.9843,376.36 79.9055,376.358 80.8268,376.357 81.748,376.355 82.6693,376.353 83.5906,376.351 84.5118,376.349 \n",
       "  85.4331,376.347 86.3543,376.344 87.2756,376.342 88.1969,376.339 89.1181,376.336 90.0394,376.333 90.9606,376.329 91.8819,376.326 92.8031,376.322 93.7244,376.317 \n",
       "  94.6457,376.313 95.5669,376.308 96.4882,376.303 97.4094,376.297 98.3307,376.291 99.252,376.285 100.173,376.278 101.094,376.271 102.016,376.263 102.937,376.255 \n",
       "  103.858,376.246 104.78,376.237 105.701,376.227 106.622,376.216 107.543,376.205 108.465,376.192 109.386,376.179 110.307,376.166 111.228,376.151 112.15,376.135 \n",
       "  113.071,376.119 113.992,376.101 114.913,376.082 115.835,376.063 116.756,376.041 117.677,376.019 118.598,375.995 119.52,375.97 120.441,375.943 121.362,375.914 \n",
       "  122.283,375.884 123.205,375.852 124.126,375.818 125.047,375.782 125.969,375.744 126.89,375.704 127.811,375.661 128.732,375.616 129.654,375.568 130.575,375.518 \n",
       "  131.496,375.465 132.417,375.408 133.339,375.349 134.26,375.286 135.181,375.22 136.102,375.15 137.024,375.076 137.945,374.998 138.866,374.916 139.787,374.829 \n",
       "  140.709,374.738 141.63,374.642 142.551,374.541 143.472,374.435 144.394,374.323 145.315,374.206 146.236,374.082 147.157,373.952 148.079,373.816 149,373.672 \n",
       "  149.921,373.522 150.843,373.364 151.764,373.198 152.685,373.024 153.606,372.842 154.528,372.651 155.449,372.45 156.37,372.24 157.291,372.021 158.213,371.791 \n",
       "  159.134,371.55 160.055,371.298 160.976,371.035 161.898,370.759 162.819,370.471 163.74,370.171 164.661,369.857 165.583,369.529 166.504,369.187 167.425,368.83 \n",
       "  168.346,368.457 169.268,368.069 170.189,367.665 171.11,367.243 172.031,366.804 172.953,366.347 173.874,365.871 174.795,365.376 175.717,364.861 176.638,364.325 \n",
       "  177.559,363.768 178.48,363.19 179.402,362.589 180.323,361.965 181.244,361.318 182.165,360.646 183.087,359.949 184.008,359.226 184.929,358.477 185.85,357.701 \n",
       "  186.772,356.897 187.693,356.064 188.614,355.202 189.535,354.311 190.457,353.388 191.378,352.434 192.299,351.448 193.22,350.43 194.142,349.378 195.063,348.291 \n",
       "  195.984,347.17 196.906,346.013 197.827,344.819 198.748,343.589 199.669,342.32 200.591,341.014 201.512,339.668 202.433,338.282 203.354,336.856 204.276,335.389 \n",
       "  205.197,333.88 206.118,332.329 207.039,330.735 207.961,329.097 208.882,327.416 209.803,325.689 210.724,323.918 211.646,322.101 212.567,320.238 213.488,318.328 \n",
       "  214.409,316.372 215.331,314.368 216.252,312.316 217.173,310.216 218.094,308.068 219.016,305.872 219.937,303.626 220.858,301.332 221.78,298.988 222.701,296.596 \n",
       "  223.622,294.154 224.543,291.662 225.465,289.122 226.386,286.532 227.307,283.892 228.228,281.204 229.15,278.468 230.071,275.682 230.992,272.849 231.913,269.968 \n",
       "  232.835,267.039 233.756,264.064 234.677,261.042 235.598,257.974 236.52,254.861 237.441,251.704 238.362,248.503 239.283,245.259 240.205,241.974 241.126,238.647 \n",
       "  242.047,235.28 242.969,231.875 243.89,228.431 244.811,224.951 245.732,221.436 246.654,217.886 247.575,214.304 248.496,210.691 249.417,207.048 250.339,203.376 \n",
       "  251.26,199.678 252.181,195.955 253.102,192.209 254.024,188.442 254.945,184.655 255.866,180.85 256.787,177.03 257.709,173.196 258.63,169.351 259.551,165.495 \n",
       "  260.472,161.633 261.394,157.765 262.315,153.895 263.236,150.024 264.157,146.154 265.079,142.289 266,138.43 266.921,134.58 267.843,130.742 268.764,126.917 \n",
       "  269.685,123.109 270.606,119.32 271.528,115.553 272.449,111.809 273.37,108.093 274.291,104.405 275.213,100.75 276.134,97.1296 277.055,93.5464 277.976,90.0032 \n",
       "  278.898,86.5025 279.819,83.047 280.74,79.6394 281.661,76.2823 282.583,72.9782 283.504,69.7297 284.425,66.5394 285.346,63.4098 286.268,60.3434 287.189,57.3428 \n",
       "  288.11,54.4102 289.031,51.5481 289.953,48.7589 290.874,46.0449 291.795,43.4083 292.717,40.8513 293.638,38.3761 294.559,35.9848 295.48,33.6793 296.402,31.4618 \n",
       "  297.323,29.334 298.244,27.2979 299.165,25.3551 300.087,23.5074 301.008,21.7563 301.929,20.1035 302.85,18.5503 303.772,17.0982 304.693,15.7484 305.614,14.502 \n",
       "  306.535,13.3603 307.457,12.3242 308.378,11.3947 309.299,10.5726 310.22,9.85857 311.142,9.25332 312.063,8.75736 312.984,8.37116 313.906,8.09505 314.827,7.92929 \n",
       "  315.748,7.87402 316.669,7.92929 317.591,8.09505 318.512,8.37116 319.433,8.75736 320.354,9.25332 321.276,9.85857 322.197,10.5726 323.118,11.3947 324.039,12.3242 \n",
       "  324.961,13.3603 325.882,14.502 326.803,15.7484 327.724,17.0982 328.646,18.5503 329.567,20.1035 330.488,21.7563 331.409,23.5074 332.331,25.3551 333.252,27.2979 \n",
       "  334.173,29.334 335.094,31.4618 336.016,33.6793 336.937,35.9848 337.858,38.3761 338.78,40.8513 339.701,43.4083 340.622,46.0449 341.543,48.7589 342.465,51.5481 \n",
       "  343.386,54.4102 344.307,57.3428 345.228,60.3434 346.15,63.4098 347.071,66.5394 347.992,69.7297 348.913,72.9782 349.835,76.2823 350.756,79.6394 351.677,83.047 \n",
       "  352.598,86.5025 353.52,90.0032 354.441,93.5464 355.362,97.1296 356.283,100.75 357.205,104.405 358.126,108.093 359.047,111.809 359.969,115.553 360.89,119.32 \n",
       "  361.811,123.109 362.732,126.917 363.654,130.742 364.575,134.58 365.496,138.43 366.417,142.289 367.339,146.154 368.26,150.024 369.181,153.895 370.102,157.765 \n",
       "  371.024,161.633 371.945,165.495 372.866,169.351 373.787,173.196 374.709,177.03 375.63,180.85 376.551,184.655 377.472,188.442 378.394,192.209 379.315,195.955 \n",
       "  380.236,199.678 381.157,203.376 382.079,207.048 383,210.691 383.921,214.304 384.843,217.886 385.764,221.436 386.685,224.951 387.606,228.431 388.528,231.875 \n",
       "  389.449,235.28 390.37,238.647 391.291,241.974 392.213,245.259 393.134,248.503 394.055,251.704 394.976,254.861 395.898,257.974 396.819,261.042 397.74,264.064 \n",
       "  398.661,267.039 399.583,269.968 400.504,272.849 401.425,275.682 402.346,278.468 403.268,281.204 404.189,283.892 405.11,286.532 406.031,289.122 406.953,291.662 \n",
       "  407.874,294.154 408.795,296.596 409.717,298.988 410.638,301.332 411.559,303.626 412.48,305.872 413.402,308.068 414.323,310.216 415.244,312.316 416.165,314.368 \n",
       "  417.087,316.372 418.008,318.328 418.929,320.238 419.85,322.101 420.772,323.918 421.693,325.689 422.614,327.416 423.535,329.097 424.457,330.735 425.378,332.329 \n",
       "  426.299,333.88 427.22,335.389 428.142,336.856 429.063,338.282 429.984,339.668 430.906,341.014 431.827,342.32 432.748,343.589 433.669,344.819 434.591,346.013 \n",
       "  435.512,347.17 436.433,348.291 437.354,349.378 438.276,350.43 439.197,351.448 440.118,352.434 441.039,353.388 441.961,354.311 442.882,355.202 443.803,356.064 \n",
       "  444.724,356.897 445.646,357.701 446.567,358.477 447.488,359.226 448.409,359.949 449.331,360.646 450.252,361.318 451.173,361.965 452.094,362.589 453.016,363.19 \n",
       "  453.937,363.768 454.858,364.325 455.78,364.861 456.701,365.376 457.622,365.871 458.543,366.347 459.465,366.804 460.386,367.243 461.307,367.665 462.228,368.069 \n",
       "  463.15,368.457 464.071,368.83 464.992,369.187 465.913,369.529 466.835,369.857 467.756,370.171 468.677,370.471 469.598,370.759 470.52,371.035 471.441,371.298 \n",
       "  472.362,371.55 473.283,371.791 474.205,372.021 475.126,372.24 476.047,372.45 476.969,372.651 477.89,372.842 478.811,373.024 479.732,373.198 480.654,373.364 \n",
       "  481.575,373.522 482.496,373.672 483.417,373.816 484.339,373.952 485.26,374.082 486.181,374.206 487.102,374.323 488.024,374.435 488.945,374.541 489.866,374.642 \n",
       "  490.787,374.738 491.709,374.829 492.63,374.916 493.551,374.998 494.472,375.076 495.394,375.15 496.315,375.22 497.236,375.286 498.157,375.349 499.079,375.408 \n",
       "  500,375.465 500.921,375.518 501.843,375.568 502.764,375.616 503.685,375.661 504.606,375.704 505.528,375.744 506.449,375.782 507.37,375.818 508.291,375.852 \n",
       "  509.213,375.884 510.134,375.914 511.055,375.943 511.976,375.97 512.898,375.995 513.819,376.019 514.74,376.041 515.661,376.063 516.583,376.082 517.504,376.101 \n",
       "  518.425,376.119 519.346,376.135 520.268,376.151 521.189,376.166 522.11,376.179 523.031,376.192 523.953,376.205 524.874,376.216 525.795,376.227 526.717,376.237 \n",
       "  527.638,376.246 528.559,376.255 529.48,376.263 530.402,376.271 531.323,376.278 532.244,376.285 533.165,376.291 534.087,376.297 535.008,376.303 535.929,376.308 \n",
       "  536.85,376.313 537.772,376.317 538.693,376.322 539.614,376.326 540.535,376.329 541.457,376.333 542.378,376.336 543.299,376.339 544.22,376.342 545.142,376.344 \n",
       "  546.063,376.347 546.984,376.349 547.906,376.351 548.827,376.353 549.748,376.355 550.669,376.357 551.591,376.358 552.512,376.36 553.433,376.361 554.354,376.362 \n",
       "  555.276,376.363 556.197,376.365 557.118,376.366 558.039,376.366 558.961,376.367 559.882,376.368 560.803,376.369 561.724,376.37 562.646,376.37 563.567,376.371 \n",
       "  564.488,376.371 565.409,376.372 566.331,376.372 567.252,376.373 568.173,376.373 569.094,376.374 570.016,376.374 570.937,376.374 571.858,376.375 572.78,376.375 \n",
       "  573.701,376.375 574.622,376.375 575.543,376.376 576.465,376.376 577.386,376.376 578.307,376.376 579.228,376.376 580.15,376.376 581.071,376.376 581.992,376.377 \n",
       "  582.913,376.377 583.835,376.377 584.756,376.377 585.677,376.377 586.598,376.377 587.52,376.377 588.441,376.377 589.362,376.377 590.283,376.377 591.205,376.377 \n",
       "  592.126,376.377 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip02)\" points=\"\n",
       "223.622,7.87402 407.874,7.87402 407.874,376.378 223.622,376.378 \n",
       "  \" fill=\"#ff0000\" fill-opacity=\"0.2\"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#e26f46; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  223.622,7.87402 407.874,7.87402 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "421.163,74.114 574.126,74.114 574.126,28.754 421.163,28.754 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  421.163,74.114 574.126,74.114 574.126,28.754 421.163,28.754 421.163,74.114 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  427.163,43.874 463.163,43.874 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 469.163, 49.8059)\" x=\"469.163\" y=\"49.8059\">Polaron Wavefn.</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:10; text-anchor:start;\" transform=\"rotate(0, 556.772, 42.9533)\" x=\"556.772\" y=\"42.9533\">2</text>\n",
       "</g>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "427.163,65.042 463.163,65.042 463.163,52.946 427.163,52.946 427.163,65.042 \n",
       "  \" fill=\"#ff0000\" fill-opacity=\"0.2\"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#e26f46; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  427.163,52.946 463.163,52.946 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 469.163, 63.494)\" x=\"469.163\" y=\"63.494\">Polaron Radius</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(rs,y.^2,label=\"Polaron Wavefn.^2\")\n",
    "plot!([0-rf,0+rf],[ψ(0).^2,ψ(0).^2],label=\"Polaron Radius\",fill=(0,:red),fillalpha=0.2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A note on Schultz\n",
    "\n",
    "In Schultz1959 (https://doi.org/10.1103%2Fphysrev.116.526), a Feynman polaron radius is defined from the Gaussian form of the Simple Harmonic Oscillator specified by the Feynman solution to the polaron problem.\n",
    "\n",
    "The reduced mass is defined in terms of the internal parameters $v$ and $w$,\n",
    "\n",
    "$\\mu=m(v^2-w^2)/v^2$.\n",
    "\n",
    "The (Gaussian) wavefunction is specified (plotted above),\n",
    "\n",
    "$\\psi(r)=(\\mu v/\\pi)^{\\frac{3}{4}} exp(-\\frac{\\mu v r^2}{2})$.\n",
    "\n",
    "And from considering the standard deviation of this wavefunction, a radius is defined: \n",
    "\n",
    "$r_f \\equiv (<p^2>)^{\\frac{1}{2}} = (3/2\\mu v)^{\\frac{1}{2}}$ (2.4, Schultz)\n",
    "\n",
    "Here $p$ is the momentum of the system, $\\mu$ is the reduced (effective) mass of the electron and interacting phonon-cloud, while $v$ and $w$ are internal polaron parameters characterising the harmonic motion of the polaron.\n",
    "The units of $v$ and $w$ are $\\hbar\\omega$.\n",
    "\n",
    "\"Using the weak... coupling expansions given by Feynman for $w$ and $v$, we find\"\n",
    "\n",
    "$r_f (\\alpha\\rightarrow 0) \\approx  (3/0.44\\alpha)^{\\frac{1}{2}} (2m\\omega)^{-\\frac{1}{2}}$\n",
    "\n",
    "To understand where this magic number (0.44) comes from, you must return to Feynman1955 and repeat the derivation used to get the small-$\\alpha$ expansion for energy. This is the part around Eqn 35 to 36 in Feynman 1955. \n",
    "\n",
    "Finding the term $2\\mu v$ is our aim.\n",
    "\n",
    "Following Feynman, assuming $\\alpha$ is small, $v=(1+\\epsilon)w$. \n",
    "You can then go back to the energy integral Feynman(31), substituting your $v$ expansion in. \n",
    "You can expand the square-root on on the denominator, by pulling out terms of $w$ to simplify, before taking a term of $\\tau^{\\frac{1}{2}}$ out to place the root in the form of (1+x)^n. You can then use a binomial expansion for this quare root, and keep only zeroth and first order in $\\epsilon$ terms. The resulting integral in the total energy has two parts, the zeroth-order contribution integrating to 1, and the linear in $\\epsilon$ term (P) apparently analytic and provided by Feynman (35),\n",
    "\n",
    "$2w^{-1} [ (1+w)^{\\frac{1}{2}} - 1] = P$.\n",
    "\n",
    "This makes the total energy,  \n",
    "\n",
    "$E=\\frac{3}{4v}(v-w)^2-A $, (Feynman33) \n",
    "\n",
    "which with our linear-in-$\\epsilon$ form of the integral $A$,\n",
    "\n",
    "$E = \\frac{3}{4v}(v-w)^2 - \\alpha \\frac{v}{w} [1-P] $\n",
    "\n",
    "Again, substituting $v=(1+\\epsilon)w$, \n",
    "\n",
    "$E=\\frac{3}{4}\\frac{((1+\\epsilon)w - w)^2}{(1+\\epsilon)w} - \\alpha\\frac{(1+\\epsilon)w}{w} [1-P]$,\n",
    "\n",
    "and discarding high order (in $\\epsilon$) terms,\n",
    "\n",
    "$E = \\frac{3}{4} w\\epsilon^2 -\\alpha -\\alpha\\epsilon[1-P]$, \n",
    "\n",
    "as given by Feynman. \n",
    "\n",
    "If we reorder terms,\n",
    "\n",
    "$E = \\frac{3}{4} w\\epsilon^2 -\\alpha\\epsilon[1-P]-(\\alpha)$, \n",
    "\n",
    "it becomes clear the total energy has a quadratic ($E =a\\epsilon^2 + b\\epsilon +c $) form in $\\epsilon$. \n",
    "\n",
    "We will want a variational (minimum energy) solution. This found at the minimum, $x=\\frac{-b}{2a}$ in the standard quadratic formula, and so, \n",
    "\n",
    "$\\epsilon = \\frac{2}{3} \\alpha \\frac{(1-P)}{w}$. \n",
    "\n",
    "This is as given by Feynman, valid for small-$\\alpha$ only.\n",
    "\n",
    "We can now return to the key $2\\mu v$ term.\n",
    "\n",
    "Expanding this identity with our definition of $\\mu$,\n",
    "\n",
    "$2\\mu v = 2m_e \\frac{v^2-w^2}{v}$.\n",
    "\n",
    "Substituting in $v=(1+\\epsilon)w$, following through the algebra, and as epsilon is small approximating $\\frac{1}{1+\\epsilon}=\\frac{1}{1}$, we eventually have\n",
    "\n",
    "$2\\mu v = 2m_e[2w\\epsilon]$.\n",
    "\n",
    "Substituting in our form for $\\epsilon$, by simple algebraic rearrangement we get\n",
    "\n",
    "$v=(1+\\frac{2\\alpha(1-P)}{3w})w = w + \\frac{2}{3} \\alpha (1-P)$\n",
    "\n",
    "Feynman states that 'the variational solution is least for $w=3$', and goes on to state that the energy correction is anyway insensitive to the choice of $w$.\n",
    "\n",
    "If $w=3$, $P=\\frac{2}{3}$, and we have that,\n",
    "\n",
    "$v=3+\\frac{2}{9}\\alpha$.\n",
    "\n",
    "Thus, $2\\mu v = 2m_e \\frac{4}{9}\\alpha = 2m_e 0.\\overline{4} \\alpha$.\n",
    "\n",
    "Schultz1959 Eqn. 2.5a approximates this as $0.44$. It is not obvious whether this arrises as a printing error, or was an intentional (1% error) approximation. \n",
    "\n",
    "It is also not clear in Schultz where the table given (reproduced above) comes from. \n",
    "Neither of the provided identities gives it. \n",
    "Possibly these data come from a direct numerical integration of $A$ and a variational solution for the parameter $v$ and $w$, using Feynman(33) as the variational energy. \n",
    "Certainly it would be worth adding this athermal energy to the (more complex) machinery at https://github.com/jarvist/PolaronMobility.jl . "
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Bibliography relevant to the above work.\n",
    "\n",
    "% The original development of Feynman's solution to the polaron problem. \n",
    "% Zero temperature approximate variational solutions developed (in limits w->0,\n",
    "% or w=~v).\n",
    "% Perturbative theories of phonon-drag effective-mass renormalisation given.\n",
    "% (i.e. where the 'me=1+alpha/6' & etc. limits are from. )\n",
    "@article{Feynman1955,\n",
    "  doi = {10.1103/physrev.97.660},\n",
    "  url = {https://doi.org/10.1103%2Fphysrev.97.660},\n",
    "  year  = {1955},\n",
    "  month = {feb},\n",
    "  publisher = {American Physical Society ({APS})},\n",
    "  volume = {97},\n",
    "  number = {3},\n",
    "  pages = {660--665},\n",
    "  author = {R. P. Feynman},\n",
    "  title = {Slow Electrons in a Polar Crystal},\n",
    "  journal = {Physical Review}\n",
    "}\n",
    "\n",
    "% Schultz seemed to spend his PhD solving the Feynman polaron problem with\n",
    "% a digital computer. \n",
    "% Lots of characterisation of the polaron state, and the introduction of an\n",
    "% effective polaron size, from considering the variance of the Gaussian\n",
    "% wavefunction. \n",
    "% Some work towards polaron mobility, but not as developed as in Feynman et al. 1962.\n",
    "% Schultz provides units for some of the quantities - which is useful!\n",
    "@article{Schultz1959,\n",
    "  doi = {10.1103/physrev.116.526},\n",
    "  url = {https://doi.org/10.1103%2Fphysrev.116.526},\n",
    "  year  = {1959},\n",
    "  month = {nov},\n",
    "  publisher = {American Physical Society ({APS})},\n",
    "  volume = {116},\n",
    "  number = {3},\n",
    "  pages = {526--543},\n",
    "  author = {T. D. Schultz},\n",
    "  title = {Slow Electrons in Polar Crystals: Self-Energy,  Mass,  and Mobility},\n",
    "  journal = {Physical Review}\n",
    "}\n",
    "\n",
    "% Original Dresselhaus (G) work on effective mass theory for excitons\n",
    "@article{Dresselhaus1956,\n",
    "  doi = {10.1016/0022-3697(56)90004-x},\n",
    "  url = {https://doi.org/10.1016/0022-3697(56)90004-x},\n",
    "  year  = {1956},\n",
    "  month = {sep},\n",
    "  publisher = {Elsevier {BV}},\n",
    "  volume = {1},\n",
    "  number = {1-2},\n",
    "  pages = {14--22},\n",
    "  author = {G. Dresselhaus},\n",
    "  title = {Effective mass approximation for excitons},\n",
    "  journal = {Journal of Physics and Chemistry of Solids}\n",
    "}"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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": 1
}