Schroedinger Equation for Hydrogen Atom¶

The Schroedinger equation is:

\begin{eqnarray} (-\frac{\hbar^2}{2m}\nabla^2-\frac{Z e^2}{4\pi\varepsilon_0 r})\psi(\vec{r})=E \psi(\vec{r}) \end{eqnarray}

using ansatz:

$\psi(\vec{r}) = Y_{lm}(\hat{r})\; u(r)/r$

and introducing dimensionless variables:

\begin{eqnarray} x = \frac{r}{r_B}\\ \varepsilon = \frac{E}{E_0} \end{eqnarray}

where \begin{eqnarray} && r_B = \frac{4\pi\varepsilon_0 \hbar^2}{m e^2} \approx 0.529 A\\ && E_0 = \frac{\hbar^2}{2 m r_B^2} == Ry \approx 13.6 eV \end{eqnarray}

we get the differential equation

\begin{eqnarray} u''(x)- \left(\frac{l(l+1)}{x^2}-\frac{2Z}{x}-\varepsilon\right)u(x)=0 \end{eqnarray}

Next we rewrite into the system of first order equations:

\begin{eqnarray} y = \left(u(x),u'(x)\right)\\ \frac{dy}{dx} = \left(u'(x),u''(x)\right) \end{eqnarray}

with boundary conditions \begin{eqnarray} &&u(0) = 0 \rightarrow \psi(0)<\infty\\ &&u(\infty)=0 \rightarrow \int |\psi(r)|^2 r^2 dr \propto \int u^2(r)dr < \infty \end{eqnarray}

Because boundary conditions are given at the two ends, we need so-called shooting method

Shooting algorithm:

Suppose the two boundary condistions are given at $a$ and $b$, i.e., $u(a)=u(b)=0$. Then

• Choose $u(a)=0$ and $u'(a)=c$, with $c$ some constant.
• Solve for $u(x)$ to the other end, and check if $u(b)=0$.
• Using root finding routine find energy $\varepsilon$ for which u(b)=0. This is the bound state.
• Continue with increasing energy $\varepsilon$ until sufficient number of bound states is found

Some remarks

• It turns out that forward integration of the radial Sch. Eq. is unstable. It is better to start integrating from infinity, and then continue down to zero.
• It is better to use logarithmic mesh for radial variable rather than linear. Radial functions need smaller number of points in logarithmic mesh

The implementation will follow these steps

1. call SciPy routine

   integrate.odeint

   to integrate the one-electron Schroedinger equation. Note that the distance is measured in units of bohr radius and energy units is Rydberg ($1 Ry = 13.6058...eV$)
2. The boundary conditions are $u(0)=0$ and $u(\infty)=0$.

   Use shooting method to obtain wave functions:

   1. Use logarithmic mesh of radial points for integration. Start integrating from a large distance ($R_{max} \sim 100$). At $R_{max}$ choose $u=0$ and some nonzero (not too large) derivative.
   2. Integrate the Schroedinger equation down to $r=0$. If your choice for the energy $\varepsilon$ corresponds to the bound state, the wave function at $u(r=0)$ will be zero.

3. Start searching for the first bound state at sufficiently negative energy (for example $\sim -1.2 Z^2$) and increase energy in sufficiently small steps to bracket all necessary bound states. Ones the wave function at $r=0$ changes sign, use root finding routine, for example

   optimize.brentq,

   to compute zero to very high precision. Store the index and the energy of the bound state for further processing.
4. Ones bound state energies are found, recompute $u(r)$ for all bound states. Normalize $u(r)$ and plot them.
5. Compute electron density for various atoms (for example He, Li, ..) neglecting Coulomb repulsion:

   Populate first $Z$ lowest laying electron states and compute $\rho = \sum_{lm\in occupied} u_{lm}^2(r)/(4\pi r^2)$. Each state with quantum number $l$ can take $2(2l+1)$ electrons. Be carefull, if atom is not one of the Nobel gases (He, Ne, ...) the last orbital is only partially filled.

Recall:

\begin{eqnarray} && y = (u(r), u'(r) )\\ && dy/dr = (u'(r), u''(r)) \end{eqnarray}\begin{eqnarray} u''(r)= \left(\frac{l(l+1)}{r^2}-\frac{2Z}{r}-\varepsilon\right)u(r) \end{eqnarray}

First we try linear mesh and forward integration. It is supposed to be unstable. We know the ground state has energy $E_0=-1 Ry$ and we should get $1s$ state with integrating Scroedinger equation.

Recal Euler's method and Runge Kutta method

Indeed the integration is unstable, and needs to be done in opposite direction. Let's try from large R.

Clearly the integration from infinity is stable, and we will use it here.

Logarithmic mesh is better suited for higher excited states, as they extend far away.

Lets create a subroutine of what we learned up to now:

Shooting algorithm:

The boundary condistions are given at two points $a$ and $b$, i.e., $u(a)=u(b)=0$.

• Choose $u(a)=0$ and $u'(a)=c$, with $c$ some constant.

• Solve for $u(x)$ to the other end, and evaluate $u(b)$.

• Using root finding routine find energy $\varepsilon$ for which u(b)=0. This is the bound state.

• Continue with increasing energy $\varepsilon$ until sufficient number of bound states is found

Shooting algorithm:

The boundary condistions are given at two points $a$ and $b$, i.e., $u(a)=u(b)=0$.

• Choose $u(a)=0$ and $u'(a)=c$, with $c$ some constant.

• Solve for $u(x)$ to the other end, and evaluate $u(b)$.

• Using root finding routine find energy $\varepsilon$ for which u(b)=0. This is the bound state.

• Continue with increasing energy $\varepsilon$ until sufficient number of bound states is found

Resulting charge density for a Ni-like Hydrogen atom

It seems that error accumulation still prevents one to calculate f-states near the nucleous. Only if we start sufficiently away from nucleout (0.1$R_b$) numerics works.

Can we do better?

Numerov algorithm¶

The general purpose integration routine is not the best method for solving the Schroedinger equation, which does not have first derivative terms.

Numerov algorithm is better fit for such equations, and its algorithm is summarized below.

The second order linear differential equation (DE) of the form

\begin{eqnarray} x''(t) = f(t) x(t) + u(t) \end{eqnarray}

is a target of Numerov algorithm.

Due to a special structure of the DE, the fourth order error cancels and leads to sixth order algorithm using second order integration scheme.

If we expand x(t) to some higher power and take into account the time reversal symmetry of the equation, all odd term cancel

\begin{eqnarray} x(h) = x(0)+h x'(0)+\frac{1}{2}h^2 x''(0)+\frac{1}{3!}h^3 x^{(3)}(0)+\frac{1}{4!}h^4 x^{(4)}(0)+\frac{1}{5!}h^5 x^{(5)}(0)+...\\ x(-h) = x(0)-h x'(0)+\frac{1}{2}h^2 x''(0)-\frac{1}{3!}h^3 x^{(3)}(0)+\frac{1}{4!}h^4 x^{(4)}(0)-\frac{1}{5!}h^5 x^{(5)}(0)+... \end{eqnarray}

hence

\begin{eqnarray} x(h)+x(-h) = 2x(0)+h^2 (f(0)x(0)+u(0))+\frac{2}{4!}h^4 x^{(4)}(0)+O(h^6) \end{eqnarray}

If we are happy with $O(h^4)$ algorithm, we can neglect $x^{(4)}$ term and get the following recursion relation

\begin{equation} x_{i+1} = 2 x_i - x_{i-1} + h^2 (f_i x_i+u_i) +O(h^4). \end{equation}

where we renaimed \begin{eqnarray} &x_{i-1}&= x(-h)\\ &x_i &= x(0)\\ &x_{i+1} &= x(h) \end{eqnarray}

But we know from the differential equation that

\begin{equation} x^{(4)} = \frac{d^2 x''(t)}{dt^2} = \frac{d^2}{dt^2}(f(t) x(t)+u(t)) \end{equation}

and we will use the well known descrete expression for the second order derivative \begin{eqnarray} g''(t) = \frac{g(t+h) - 2 g(t) + g(t-h)}{h^2} + O(h^2) \end{eqnarray} which can be approximated by

\begin{equation} x^{(4)}= \frac{f_{i+1}x_{i+1}+u_{i+1}-2 f_i x_i -2 u_i+ f_{i-1}x_{i-1}+u_{i-1}}{h^2} + O(h^2) \end{equation}

Inserting the fourth order derivative into the above recursive equation (forth equation in his chapter), we get

\begin{equation} x_{i+1}-2 x_i+x_{i-1}=h^2(f_i x_i+u_i)+\frac{h^2}{12}(f_{i+1}x_{i+1}+u_{i+1}-2 f_i x_i -2 u_i+ f_{i-1}x_{i-1}+u_{i-1}) + O(h^6) \end{equation}

If we switch to a new variable $w_i=x_i(1-\frac{h^2}{12} f_i)-\frac{h^2}{12}u_i$ we are left with the following equation

\begin{equation} w_{i+1} = 2 w_i - w_{i-1} + h^2 (f_i x_i + u_i)+O(h^6) \end{equation}

The variable $x$ needs to be recomputed at each step with \begin{equation} x_i=\frac{w_i+\frac{h^2}{12}u_i}{1-\frac{h^2}{12}f_i}. \end{equation}

For Numerov algorithm we can evaluate derivative part $f(r)$ for all points at once:

Because Schroedinger Eq is: \begin{eqnarray} u''(r)= \left(\frac{l(l+1)}{r^2}-\frac{2Z}{r}-\varepsilon\right)u(r) \end{eqnarray} the $f$ function is \begin{eqnarray} f(r)= \left(\frac{l(l+1)}{r^2}-\frac{2Z}{r}-\varepsilon\right) \end{eqnarray}

The Numerov algorithm is much faster, but the price we pay is the mesh has to be linear. We can not use logarithmic mesh in combination with Numerov algorithm.

Numerov seems much more precise than odeint, and avoids numerical problems we had before

Put it all together

It seems with Numerov we are getting substantial error-bar for the energy of 1s state. We could increase the number of points in the mesh, but the error decreases only linearly with the number of points used.

Where is the problem? What should be done?

Check that approximate solution gives smaller wave function at zero than exact energy, which confirms that root finding routine works fine.

Let's check how the function looks like near zero

Idea: The mesh is very sparse near zero, and in the range of the first few points, the curve is not linear enough. Linear extrapolation gives the error.

Can we do better?

Let's use cubic extrapolation with first 4 points.

Indeed we get $10^{-8}$ error as compared to $10^{-5}$ error before. So, the extrapolation must be improved to reduce the error.

Increasing the number of points for 10-times reduces the error factor of 1000, i.e., $O(1/N^3)$.

Better cubic extrapolation reduces the error for the same factor of 1000, equivalent to 10-times more points.

It also helps to increase the number of points. 10-times denser grid gives roughly $10^{-3}$ smaller error.