Binding energy for a uniform sphere bound by a Yukawa potential rather than an inverse-square-law potential?

Question

The binding energy of a uniform sphere that interacts through a force that follows the inverse square law is proportional to $a\frac{C^2}{R}$, where $C$ is whatever charge determines the strength of the force on the particle, such as mass or electric charge, $R$ is the radius, and $a$. For gravity, this looks like $\frac{3GM^2}{5R}$. What would the equation be for a Yukawa potential though? I thought the equation might be $a\frac{C^2e^{-bR}}R$, where $a$ and $b$ are constants, but when I predicted the binding energy by using an integral to calculate potential for a given radius and plotting the results, I couldn't get $a\frac{C^2e^{-bR}}R$ to fit the data.

The equation I used to predict the binding energy of the sphere for a given radius was this: $$\frac{\int_0^R4\pi r^2\left(\int_{-R}^R\left(\int_0^\sqrt{1-x^2}\left(2\pi y \frac{e^{-\sqrt{(x-r)^2+y^2}}}{\sqrt{(x-r)^2+y^2}}\right)dy\right)dx\right)dr}{\left(\frac43\pi R^3\right)^2}$$ The outermost integral variable, $r$, represents a radius from the center of the sphere that is less than the radius $R$. The $4\pi r^2$ bit is because each radius $r$ also represents an equipotential shell. The rest of the equation treats $r$ like a distance to a certain point, and determines the potential at that point. $\int_{-R}^R\int_0^\sqrt{1-x^2}2\pi y...dydx$ maps each point in the sphere and calculates each point's "contribution" to the potential at point $r$. I divide the whole thing by the volume squared because I want the total "charge" to remain constant.

Answer 1

Per this Math.SE post, given a Yukawa potential with range $1/a$, a spherical shell of mass $m$ with radius $r$ behaves like a point mass potential at the origin, except with mass $m~\sinh(ar)/(ar).$

If we nest these then the binding energy can be found as follows: bring a small mass from infinity to the origin for free, forming a little ball of density $\rho$ and radius $\delta r$, bring another small mass from infinity to the radius $\delta r$ to grow the sphere a little more, paying some energy, and so on. To keep a constant density the mass we have to move each time is $\mathrm dm=\rho~4\pi r^2\mathrm dr,$ while the mass already moved is $m(r)=\frac43\pi r^3~\rho$. So we need to know two things:

First, what is the effective mass of a uniform ball, instead of a spherical shell? That's a straightforward integration by parts:$$m_{\text{eff}}(R)=\int_0^R\mathrm dr ~4\pi r^2\rho~\frac{\sinh(ar)}{ar} = \frac{4\pi\rho}{a^3}~\big(aR\cosh(aR)-\sinh(aR)\big) $$ Then assuming your potential looks like $$U(r)=-G\frac{m_1m_2}{r}~e^{-ar}$$ the binding energy is some integral $$E= G\int_0^R\mathrm dr~4\pi r^2\rho~m_{\text{eff}}(r)\frac{e^{-ar}}{r}$$ Calling $4\pi\rho/a^3=\mu$ this gives $$E= G\mu^2~a\int_0^{aR}\mathrm dx~x~e^{-x}(x \cosh x-\sinh x)$$This is technically not a hard integral, it is just a couple more integrations by parts and the integrand actually simplifies when you substitute in the definition of cosh/sinh... But it is a rather big expression so I trust Wolfram Alpha more than I trust my pen and paper, and Wolfram Alpha says it is $$E= {G\mu^2a \over 12}\left( 2 (aR)^3 - 3 (aR)^2 - 3 e^{-2 aR} (aR + 1)^2 + 3 \right)$$ Wolfram Alpha is also kind enough to give us a check on this expression via the Taylor series, which it says is$$E=G\mu^2a\left(\frac{(aR)^5}{15} -\frac{(aR)^6}{18}+\frac{(aR)^7}{35}-\dots\right)$$ which indeed limits in the infinite range $a=0$ case to $(3/5) GM^2/R$ (remember that $\mu\sim1/a^3$).

So yeah, I am not surprised that you could not model it that easily, it has a rather complicated form!