10

I am trying to work out the polarisability of the ground state hydrogen atom using second-order Rayleigh-Schrodinger perturbation theory, including the unbound states. I know the angular part can be evaluated easily and that only the $Y_{10}$ states contributes, so the expression can be simplified to, $$ \alpha=\frac{2}{3}\sum_{n=2}^{\infty}\frac{|\langle R_{1 {\rm s}}|r| R_{n {\rm p}}\rangle|^2} {E_n+1/2} + \frac{2}{3} \int_{0}^{\infty} \frac{|\langle R_{1 {\rm s}}|r| R_{E {\rm p}}\rangle|^2} {E+1/2} \mathrm{d} E,$$ where the sum is over the bound states and the integral is over the continuum states and $R$ are the radial wavefunctions. I've evaluated the bound state sum up to $n = 150$ and it converges on $\approx 3.66$, which I believe to be correct. I know that the analytical answer can be obtained from the exact first-order wavefunction and is 4.5, meaning the contribution from the integral term should be $\approx 0.84$, but I can't get that value out.

I've looked up a form for the unbound wavefunctions as $$ R_{El} (r) = \sqrt{\frac{4E}{\pi}} \mathrm{e}^{-\frac{\pi}{2 k}} \left| \Gamma \left(l + 1 - \frac{\mathrm{i}}{k}\right) \right| \frac{1}{\left(2 l + 1\right)!} \left(2 k r\right)^l \mathrm{e}^{- \mathrm{i} k r} ~_1F_1\left( l + 1 - \frac{\mathrm{i}}{k}, 2 l + 2, 2 \mathrm{i} k r \right), $$ where $k = \sqrt{2E}$ is the wavenumber and $_1F_1$ is a confluent hypergeometric function of the first kind. I've got a short C++ program to do the numerical integration, using FLINT Arb for the wavefunction evaluation and GSL for the integration, which seems to converge fine if I do the dipole integral between 0 and 100, and the energy integral from 0.001 to 100 (assuming units of Bohr and Hartree).

I just can't get the right value out. I see from various sources that there seems to be disagreement of the correct form for the unbound states, particularly around the normalisation, so I've tried several powers of $E$ and various factors of $\pi$ and $2$ at the front, and can get values around 0.5 to 2, but none match the analytical value.

Can anyone tell what form I need to use for the unbound states to get the analytical value out, or if I'm making some more serious error?

Related to Scattering states of Hydrogen atom in non-relativistic perturbation theory

Qmechanic
  • 220,844
John
  • 203

1 Answers1

6

There are some weird signs and factors in your expression for the radial function. I used Eq. (4.78) of Jentschura & Adkins: $$ R_{El}(r)=\sqrt{\frac{2k}{\pi}}\frac{(2kr)^l}{(2l+1)!}\exp\left(\frac{\pi}{2k}\right)\left|\Gamma\left(l+1+\frac{i}{k}\right)\right|\exp(-ikr) \, _1F_1\left(l+1+\frac{i}{k},2l+2,2ikr\right) \ , $$ so for $l=1$: $$ R_{E1}(r)=\frac{1}{3}\sqrt{\frac{2k}{\pi}}kr\exp\left(\frac{\pi}{2k}\right)\left|\Gamma\left(2+\frac{i}{k}\right)\right|\exp(-ikr) \, _1F_1\left(2+\frac{i}{k},4,2ikr\right) \ . $$ These scattering states are normalized as $\langle\psi_{Elm}|\psi_{E'l'm'}\rangle=\delta(E-E')\delta_{ll'}\delta_{mm'}$.

I computed both the $r$ integrals: $$ g(E)=2\int_0^\infty\mathrm{d}r \, r^3 \, \exp(-r) R_{E1}(r) \ , $$ and the $E$ integral: $$ S=\frac{4}{3}\int_0^\infty\mathrm{d}E\frac{g^2(E)}{2E+1} $$ numerically with Mathematica, and got $S\approx0.836736$. Based on my previous answer, the sum-over-bound-states part of the polarizability is $-2\times(-1.83163)=3.66326$, so the continuum part should give $9/2-3.66326=0.83674$, and this is exactly what I got from the numerical integration.

Maybe try calculating your integrals with the above radial function. Let me know if that helped.

dennismoore94
  • 3,490