# I am failing to verify an analytical integral with Mathematica

I am trying to use Mathematica to check an integral from a book, and cannot get it done. The book talks about a spherically symmetric mass distribution in $$R^3$$ of the form

$$\rho(r)=r^{-\gamma}, \gamma<3$$

and it goes on by stating that the Newtonian potential to that distribution is given by

$$\phi(r) = 4π\ln r + \mathrm{const.}, \gamma=2$$

and

$$\phi(r) = \frac{4π}{(2-\gamma)(3-\gamma)}r^{2-\gamma}+\mathrm{const.}, \gamma\neq2$$ .

(There are a couple more constants in the book which I set to 1 here for simplicity.)

I wanted to verify this by directly evaluating the integral for the potential

$$\phi(\mathbf x) = -\int_{R^3}\frac{\rho(\mathbf x')}{|\mathbf x-\mathbf x'|}\,d^3x' =-4π\int_0^\infty\frac{\rho(r')}{|r-r'|}r'^2\,dr'$$ ,

and I am failing.

Here is my code to check the $$\gamma=2$$ case:

Assuming[{r > 0, g == 2}, -4 π Integrate[
s^(-g) / Abs[r - s] s^2, {s, 0, ∞}]]


Mathematica gives me the message: "Integrate: Integral of 1/Abs[r-s] does not converge on {0,[Infinity]}"

My code for the other case is

Assuming[{r > 0, g != 2, g < 3}, -4 π Integrate[
s^(-g) / Abs[r - s] s^2, {s, 0, ∞}]]


which also yields "Integrate: Integral of s^(2-g)/Abs[r-s] does not converge on {0, ∞}."

What am I missing?

• One does problems of this sort with the help of the Gauss theorem. The way you try it, is much too difficult, and directly it will not work. Mar 3 '20 at 16:15
• Hi @AlexeiBoulbitch , thanks for the answer! Sure, there might be simpler ways than directly evaluating the integral. I wanted to try this, since later I want to generalise this to more complex distributions and integrate numerically. But I found a mistake now: I cannot simply replace $|\mathbf x-\mathbf x'|$ in the 3D integrand by $|r-r'|$ in the 1D integrand! I will correct this and try again tomorrow and then report back! Mar 3 '20 at 22:42
• FullSimplify[Integrate[s^(2-g)/Sqrt[(r-s)^2],s]] does this help? You can now try Limit and Series with Assuming. Mar 4 '20 at 5:44
• @Britzel If your (more complex) distributions are spherically- or cylindrically symmetric, the Gauss theorem approach is the way to implement. Do you know it? If you have an asymmetric problem you will have a bad time. I am just struggling against a comparable problem, where I need to integrate a Green function over a distribution. [Here][1] you can see the level of mathematical difficulties arising there. Have fun! [1]: mathematica.stackexchange.com/questions/215614/… Mar 4 '20 at 8:58
• @AlexeiBoulbitch For the start they will be spherically symmetric, yes. Perhaps it is a good idea to go with a Gauß theorem approach for now, as you suggest. Eventually though they will become axisymmetric. But not for now. Mar 4 '20 at 13:16

The problem is with the expansion: to integrate $$\frac{1}{\lvert \mathbf{x'} - \mathbf{x}\rvert}$$ over a volume with a spherically invariant factor (the charge density), the standard (and the most sensible trick) is to apply spherical harmonic decomposition:

$$\frac{1}{\lvert \mathbf{x'} - \mathbf{x}\rvert}=\sum\limits_{l=0}^\infty\frac{r_<^l}{r_>^{l+1}}\left(\frac{4\pi}{2l+1}\right)\sum\limits_{m=-l}^{l}Y_{lm}(\theta,\phi)Y_{lm}^*(\theta',\phi')$$

With this expansion, the integral will be much simpler; in particular, the angular integral will immediately fix radial dependence. In fact, for spherically symmetric overall factor, the angular integration simply reads

\begin{align} \int d\Omega'\rho(r')\frac{1}{\lvert \mathbf{x'} - \mathbf{x}\rvert}=&\sum\limits_{l=0}^\infty\frac{r_<^l}{r_>^{l+1}}\left(\frac{4\pi}{2l+1}\right)\sum\limits_{m=-l}^{l}Y_{lm}(\theta,\phi)\rho(r')\int d\Omega'Y_{lm}^*(\theta',\phi')\\ =&\sum\limits_{l=0}^\infty\frac{r_<^l}{r_>^{l+1}}\left(\frac{4\pi}{2l+1}\right)\sum\limits_{m=-l}^{l}Y_{lm}(\theta,\phi)\rho(r')\sqrt{4\pi}\delta_l^0\delta_m^0\\ =&\frac{4\pi}{r_>}\rho(r') \end{align}

The radial integral is now straightforward:

\begin{align} \int d^3x'\frac{\rho(r')}{\lvert \mathbf{x'} - \mathbf{x}\rvert}=&\int r'^2 dr'\int d\Omega'\rho(r')\frac{1}{\lvert \mathbf{x'} - \mathbf{x}\rvert}\\ =&4\pi\int \rho(r')r' dr'\\ =&4\pi\int r'^{1-\gamma} dr' \end{align}