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$


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

  • 1
    $\begingroup$ 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. $\endgroup$ Commented Mar 3, 2020 at 16:15
  • $\begingroup$ 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! $\endgroup$
    – Britzel
    Commented Mar 3, 2020 at 22:42
  • $\begingroup$ FullSimplify[Integrate[s^(2-g)/Sqrt[(r-s)^2],s]] does this help? You can now try Limit and Series with Assuming. $\endgroup$ Commented Mar 4, 2020 at 5:44
  • $\begingroup$ @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/… $\endgroup$ Commented Mar 4, 2020 at 8:58
  • $\begingroup$ @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. $\endgroup$
    – Britzel
    Commented Mar 4, 2020 at 13:16

1 Answer 1


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}$$

which yields your expected result.

  • 1
    $\begingroup$ Fantastic, thanks! I will have to look back into spherical harmonics in order to understand the angular integration you performed. Will do so now! $\endgroup$
    – Britzel
    Commented Mar 4, 2020 at 13:18
  • 1
    $\begingroup$ Soner, can you show how you did this with the Wolfram Language? $\endgroup$ Commented Mar 4, 2020 at 18:47
  • 1
    $\begingroup$ @CATrevillian I am unaware of how to do intermediate steps I wrote above with the Wolfram Language. $\endgroup$ Commented Mar 5, 2020 at 4:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.