Hooke's Atom

This system is of significance as it is, for certain values of the force constant defining the harmonic containment, an exactly solvable ground-state many-electron problem that explicitly includes electron correlation. As such it can provide insight into quantum correlation (albeit in the presence of a non-physical nuclear potential) and can act as a test system for judging the accuracy of approximate quantum chemical methods for solving the Schrödinger equation. The name "Hooke's atom" arises because the harmonic potential used to describe the electron-nucleus interaction is a consequence of Hooke's law.

Employing atomic units, the Hamiltonian defining the Hooke's atom is

${\displaystyle {\hat {H}}=-{\frac {1}{2}}\nabla _{1}^{2}-{\frac {1}{2}}\nabla _{2}^{2}+{\frac {1}{2}}k(r_{1}^{2}+r_{2}^{2})+{\frac {1}{|\mathbf {r} _{1}-\mathbf {r} _{2}|}}.}$ 

As written, the first two terms are the kinetic energy operators of the two electrons, the third term is the harmonic electron-nucleus potential, and the final term the electron-electron interaction potential. The non-relativistic Hamiltonian of the helium atom differs only in the replacement:

${\displaystyle -{\frac {2}{r}}\rightarrow {\frac {1}{2}}kr^{2}.}$ 

The equation to be solved is the two electron Schrödinger equation:

${\displaystyle {\hat {H}}\Psi (\mathbf {r} _{1},\mathbf {r} _{2})=E\Psi (\mathbf {r} _{1},\mathbf {r} _{2}).}$ 

For arbitrary values of the force constant, k, the Schrödinger equation does not have an analytic solution. However, for a countably infinite number of values, such as k=¼, simple closed form solutions can be derived. Given the artificial nature of the system this restriction does not hinder the usefulness of the solution.

To solve, the system is first transformed from the Cartesian electronic coordinates, (r₁,r₂), to the center of mass coordinates, (R,u), defined as

${\displaystyle \mathbf {R} ={\frac {1}{2}}(\mathbf {r} _{1}+\mathbf {r} _{2}),\mathbf {u} =\mathbf {r} _{2}-\mathbf {r} _{1}.}$ 

Under this transformation, the Hamiltonian becomes separable – that is, the |r₁ - r₂| term coupling the two electrons is removed (and not replaced by some other form) allowing the general separation of variables technique to be applied to further a solution for the wave function in the form ${\displaystyle \Psi (\mathbf {r} _{1},\mathbf {r} _{2})=\chi (\mathbf {R} )\Phi (\mathbf {u} )}$ . The original Schrödinger equation is then replaced by:

${\displaystyle \left(-{\frac {1}{4}}\nabla _{\mathbf {R} }^{2}+kR^{2}\right)\chi (\mathbf {R} )=E_{\mathbf {R} }\chi (\mathbf {R} ),}$ 
${\displaystyle \left(-\nabla _{\mathbf {u} }^{2}+{\frac {1}{4}}ku^{2}+{\frac {1}{u}}\right)\Phi (\mathbf {u} )=E_{\mathbf {u} }\Phi (\mathbf {u} ).}$ 

The first equation for ${\displaystyle \chi (\mathbf {R} )}$  is the Schrödinger equation for an isotropic quantum harmonic oscillator with ground-state energy ${\displaystyle E_{\mathbf {R} }=(3/2){\sqrt {k}}E_{\mathrm {h} }}$  and (unnormalized) wave function

${\displaystyle \chi (\mathbf {R} )=e^{-{\sqrt {k}}R^{2}}.}$ 

Asymptotically, the second equation again behaves as a harmonic oscillator of the form ${\displaystyle \exp(-({\sqrt {k}}/4)u^{2})\,}$  and the rotationally invariant ground state can be expressed, in general, as ${\displaystyle \Phi (\mathbf {u} )=f(u)\exp(-({\sqrt {k}}/4)u^{2})\,}$  for some function ${\displaystyle f(u)\,}$ . It was long noted that f(u) is very well approximated by a linear function in u. Thirty years after the proposal of the model an exact solution was discovered for k=¼, and it was seen that f(u)=1+u/2. It was later shown that there are many values of k which lead to an exact solution for the ground state, as will be shown in the following.

Decomposing ${\displaystyle \Phi (\mathbf {u} )=R_{l}(u)Y_{lm}}$  and expressing the Laplacian in spherical coordinates,

${\displaystyle \left(-{\frac {1}{u^{2}}}{\frac {\partial }{\partial u}}\left(u^{2}{\frac {\partial }{\partial u}}\right)+{\frac {{\hat {L}}^{2}}{u^{2}}}+{\frac {1}{4}}ku^{2}+{\frac {1}{u}}\right)R_{l}(u)Y_{lm}({\hat {\mathbf {u} }})=E_{l}R_{l}(u)Y_{lm}({\hat {\mathbf {u} }}),}$ 

one further decomposes the radial wave function as ${\displaystyle R_{l}(u)=S_{l}(u)/u\,}$  which removes the first derivative to yield

${\displaystyle -{\frac {\partial ^{2}S_{l}(u)}{\partial u^{2}}}+\left({\frac {l(l+1)}{u^{2}}}+{\frac {1}{4}}ku^{2}+{\frac {1}{u}}\right)S_{l}(u)=E_{l}S_{l}(u).}$ 

The asymptotic behavior ${\displaystyle S_{l}(u)\sim e^{-{\frac {\sqrt {k}}{4}}u^{2}}\,}$  encourages a solution of the form

${\displaystyle S_{l}(u)=e^{-{\frac {\sqrt {k}}{4}}u^{2}}T_{l}(u).}$ 

The differential equation satisfied by ${\displaystyle T_{l}(u)\,}$  is

${\displaystyle -{\frac {\partial ^{2}T_{l}(u)}{\partial u^{2}}}+{\sqrt {k}}u{\frac {\partial T_{l}(u)}{\partial u}}+\left({\frac {l(l+1)}{u^{2}}}+{\frac {1}{u}}+\left({\frac {\sqrt {k}}{2}}-E_{l}\right)\right)T_{l}(u)=0.}$ 

This equation lends itself to a solution by way of the Frobenius method. That is, ${\displaystyle T_{l}(u)\,}$  is expressed as

${\displaystyle T_{l}(u)=u^{m}\sum _{k=0}^{\infty }\ a_{k}u^{k}.}$ 

for some ${\displaystyle m\,}$  and ${\displaystyle \{a_{k}\}_{k=0}^{k=\infty }\,}$  which satisfy:

${\displaystyle m(m-1)=l(l+1)\,,}$ 

${\displaystyle a_{0}\neq 0\,}$ 

${\displaystyle a_{1}={\frac {a_{0}}{2(l+1)}},}$ 

${\displaystyle a_{2}={\frac {a_{1}+\left({\sqrt {k}}(l+{\frac {3}{2}})-E_{l}\right)a_{0}}{2(2l+3)}}={\frac {a_{0}}{2(2l+3)}}\left({\frac {1}{2(l+1)}}+{\sqrt {k}}\left(l+{\frac {3}{2}}\right)-E_{l}\right),}$ 

${\displaystyle a_{3}={\frac {a_{2}+\left({\sqrt {k}}(l+{\frac {5}{2}})-E_{l}\right)a_{1}}{6(l+2)}},}$ 

${\displaystyle a_{n+1}={\frac {a_{n}+\left({\sqrt {k}}(l+{\frac {1}{2}}+n)-E_{l}\right)a_{n-1}}{(n+1)(2l+2+n)}}.}$ 

The two solutions to the indicial equation are ${\displaystyle m=l+1}$  and ${\displaystyle m=-l}$  of which the former is taken as it yields the regular (bounded, normalizable) wave function. For a simple solution to exist, the infinite series is sought to terminate and it is here where particular values of k are exploited for an exact closed-form solution. Terminating the polynomial at any particular order can be accomplished with different values of k defining the Hamiltonian. As such there exists an infinite number of systems, differing only in the strength of the harmonic containment, with exact ground-state solutions. Most simply, to impose a_k = 0 for k ≥ 2, two conditions must be satisfied:

${\displaystyle {\frac {1}{2(l+1)}}+{\sqrt {k}}\left(l+{\frac {3}{2}}\right)-E_{l}=0,}$ 
${\displaystyle {\sqrt {k}}(l+{\frac {5}{2}})=E_{l}.}$ 

These directly force a₂ = 0 and a₃ = 0 respectively, and as a consequence of the three term recession, all higher coefficients also vanish. Solving for ${\displaystyle {\sqrt {k}}\,}$  and ${\displaystyle E_{l}\,}$  yields

${\displaystyle {\sqrt {k}}={\frac {1}{2(l+1)}},}$ 
${\displaystyle E_{l}={\frac {2l+5}{4(l+1)}},}$ 

and the radial wave function

${\displaystyle T_{l}=u^{l+1}\left(a_{0}+{\frac {a_{0}}{2(l+1)}}u\right).}$ 

Transforming back to ${\displaystyle R_{l}(u)\,}$ 

${\displaystyle R_{l}(u)={\frac {T_{l}(u)e^{-{\frac {\sqrt {k}}{4}}u^{2}}}{u}}=u^{l}\left(1+{\frac {1}{2(l+1)}}u\right)e^{-{\frac {\sqrt {k}}{4}}u^{2}},}$ 

the ground-state (with ${\displaystyle l=0\,}$  and energy ${\displaystyle 5/4E_{\mathrm {h} }\,}$ ) is finally

${\displaystyle \Phi (\mathbf {u} )=\left(1+{\frac {u}{2}}\right)e^{-u^{2}/8}.}$ 

Combining, normalizing, and transforming back to the original coordinates yields the ground state wave function:

${\displaystyle \Psi (\mathbf {r} _{1},\mathbf {r} _{2})={\frac {1}{2{\sqrt {8\pi ^{5/2}+5\pi ^{3}}}}}\left(1+{\frac {1}{2}}|\mathbf {r} _{1}-\mathbf {r} _{2}|\right)\exp \left(-{\frac {1}{4}}{\big (}r_{1}^{2}+r_{2}^{2}{\big )}\right).}$ 

The corresponding ground-state total energy is then ${\displaystyle E=E_{R}+E_{u}={\frac {3}{4}}+{\frac {5}{4}}=2E_{\mathrm {h} }}$ .

The exact ground state electronic density of the Hooke atom for the special case ${\displaystyle k=1/4}$  is

${\displaystyle \rho (\mathbf {r} )={\frac {2}{\pi ^{3/2}(8+5{\sqrt {\pi }})}}e^{-(1/2)r^{2}}\left(\left({\frac {\pi }{2}}\right)^{1/2}\left({\frac {7}{4}}+{\frac {1}{4}}r^{2}+\left(r+{\frac {1}{r}}\right)\mathrm {erf} \left({\frac {r}{\sqrt {2}}}\right)\right)+e^{-(1/2)r^{2}}\right).}$ 

From this we see that the radial derivative of the density vanishes at the nucleus. This is in stark contrast to the real (non-relativistic) helium atom where the density displays a cusp at the nucleus as a result of the unbounded Coulomb potential.

• List of quantum-mechanical systems with analytical solutions