# 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. Commented Mar 3, 2020 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! Commented Mar 3, 2020 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. Commented Mar 4, 2020 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/… Commented Mar 4, 2020 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. Commented Mar 4, 2020 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}