I'm studying many-electron atoms without e-e repulsion, in particular, an atom with two electrons and Z protons. Once we get the wave function for the ground state, which (in atomic units) is $$\psi (r_1,r_2)= {Z^3 \over \pi } e^{-Z(r_1+r_2)},$$
we want to approximate how much the e-e interaction would contribute to the energy. For that purpuse, we calculate the expected value of the potential energy $ ~ U_{ee} ={1 \over | \bf r _1 – r_2 |} $ (again in atomic units), in the state given by the previous wave function. The integral to be calculated is: $$ \int d ^3 {\bf r } _1 d ^3 {\bf r } _2 ~ \psi ^* (r_1,r_2) {1 \over | \bf r _1 – r_2 |} \psi (r_1,r_2) = {Z^6 \over \pi^2} \int d ^3 {\bf r } _1 d ^3 {\bf r } _2 ~ { e^{-2Z(r_1+r_2)} \over | \bf r _1 – r_2 |} $$
Now, it's been a while since I took a course in multivariable calculus, but I gave it a try… and couldn't do it. I know that the integral in cartesian coordinates is $${Z^6 \over \pi^2} \int _{-\infty} ^\infty \int _{-\infty} ^\infty\int _{-\infty} ^\infty\int _{-\infty} ^\infty\int _{-\infty} ^\infty\int _{-\infty} ^\infty { e^{-2Z \left(\sqrt{x_1^2+y_1^2+z_1^2}+\sqrt{x_2^2+y_2^2+z_2^2} \right) } \over \sqrt{(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2}}dx_1dy_1dz_1dx_2dy_2dz_2 , $$
and that, in order to compute it, you'd need to do the six iterated integrals. Mathematica can't do that in a reasonable time, and neither can I (I doubt there is a closed form for most of them). So instead I tried two sets of spherical coordinates: $$ x_1=r_1 \sin \theta _1 \cos \phi _1 ~~~~~~~~ x_2=r_2 \sin \theta _2 \cos \phi _2 \\ y_1=r_1 \sin \theta _1 \sin \phi _1 ~~~~~~~~~ y_2=r_2 \sin \theta _2 \sin \phi _2 \\ z_1=r_1 \cos \theta _1 ~~~~~~~~~ z_2=r_2 \cos \theta _2 $$
I used mathematica to calculate the Jacobian determinant and it turns out it's just the product of the two individual Jacobians (which should have been obvious but I'm a bit dumb). After calculating that denominator, the integral that needs to be done is:
$$ {Z^6 \over \pi ^2} \int ^{2\pi} _0 \int ^{2\pi} _0 \int ^\pi _0 \int ^\pi _0 \int ^\infty _0 \int ^\infty _0 {e^{-2Z(r_1+r_2)} r_1^2r_2^2 \sin \theta _1 \sin \theta _2 \over \sqrt{r_1^2+r_2^2-2r_1r_2 (\sin \theta _1 \sin \theta _2 \cos (\phi _1 – \phi _2)+\cos \theta _1 \cos \theta _2)}} dr_1dr_2d\theta_1d\theta_2d\phi_1d\phi_2$$
…which looks… better? Mathematica still hangs up if I try to do the whole integral. Integrating first in $r_1$ or $r_2$ seems impossible analytically. Trying to integrate first in the rest of variables gives a closed answer, but it's in terms of hypergeometric functions or elliptic integrals (and it's absurdly complicated). And the thing is: the answer to this integral is suposedly ${5 \over 8}Z$. Is it really necessary to do six integrals involving hypergeometric functions and elliptic integrals to get such a simple answer? I wouldn't know how to continue anyways, because well… I don't know how to integrate eliptic integrals composed with other complicated functions and it seems mathematica doesn't know either. What should I do here?
Edit: I swear the integrals fit in one line when I was writing.
Best Answer
One way of computing it is by using the expansion \begin{equation} \frac{1}{|\mathbf{r_1} - \mathbf{r_2}|} = 4 \pi \sum_{l=0}^{\infty} \sum_{m=-l}^l \frac{1}{2l+1} \frac{r_<^l}{r_>^{l+1}} Y^*_{lm}(\theta_2,\phi_2) Y_{lm}(\theta_1,\phi_1) \, , \end{equation} where $r_>$ ($r_<$) is the maximum (minimum) of $r_1$ and $r_2$. You can find it as equation $(3.70)$ in Jackson's Classical Electrodynamics. Plugging it into the integral makes the angular integrations trivial, thanks to the orthonormalization of spherical harmonics when integrated over the whole solid angle, \begin{equation} \int \mathrm{d} S \, Y^*_{l m}(\theta,\phi) Y_{l' m'}(\theta,\phi) = \delta_{l l'} \delta_{m m'} \, . \end{equation} Hence we have, \begin{eqnarray} \langle \frac{1}{|\mathbf{r_1} - \mathbf{r_2}|} \rangle &=& \left(\frac{Z^3}{\pi}\right)^2 4 \pi \sum_{l=0}^{\infty} \sum_{m=-l}^l \frac{1}{2l+1} \int_0^\infty \mathrm{d} r_1 \int_0^\infty \mathrm{d} r_2 \frac{r_<^l}{r_>^{l+1}} \mathrm{e}^{-2 Z(r_1 + r_2)} r_1^2 r_2^2 \int \mathrm{d} S_1 Y_{lm}(\theta_1,\phi_1) \int \mathrm{d} S_2 Y^*_{lm}(\theta_2,\phi_2) \\ &=& 4 \frac{Z^6}{\pi} \int_0^\infty \mathrm{d} r_1 \int_0^\infty \mathrm{d} r_2 \frac{1}{r_>} \mathrm{e}^{-2 Z(r_1 + r_2)} r_1^2 r_2^2 \times \sqrt{4 \pi} \times \sqrt{4 \pi} \\ &=& 16 Z^6 \int_0^\infty \mathrm{d} r_2 \left[ \int_0^{r_2} \mathrm{d} r_1 \frac{1}{r_2} \mathrm{e}^{-2 Z(r_1 + r_2)} r_1^2 r_2^2 + \int_{r_2}^\infty \mathrm{d} r_1 \frac{1}{r_1} \mathrm{e}^{-2 Z(r_1 + r_2)} r_1^2 r_2^2 \right] \\ &=& \frac{5}{8} Z \,, \end{eqnarray} where we have used the fact that only the $l = 0$ and $m = 0$ term contributes.