Hooke's atom, also known as harmonium or hookium, refers to an artificial helium-like atom where the Coulombic electron-nucleus interaction potential is replaced by a harmonic potential.^[1][2] This system is of significance as it is, for certain values of the force constant defining the harmonic containment, an exactly solvable^[3] 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.^[4][5] 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.^[5] 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.^[2] Thirty years after the proposal of the model an exact solution was discovered for k=¼,^[3] 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,^[5] 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^[4]

${\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

• Cioslowski, Jerzy; Pernal, Katarzyna (2000). "The Ground State of Harmonium". The Journal of Chemical Physics. 113 (19): 8434–8443. Bibcode:2000JChPh.113.8434C. doi:10.1063/1.1318767.
• O’Neill, Darragh P.; Gill, Peter M. W. (2003). "Wave functions and two-electron probability distributions of the Hooke's-law atom and helium" (PDF). Physical Review A. 68 (2): 022505. Bibcode:2003PhRvA..68b2505O. doi:10.1103/PhysRevA.68.022505.