Potential energy of a sphere in the field of itself

[Orodruin]

So you're saying that post #30 is nonsense because we are trying to solve Poisson's equation here? (I dread the answer).

No, I attempted "humor". I may have erroneously assumed that many people would be at least vaguely aware that those are among the contents of my textbook. Hence, saying that there is no good textbook on the subject would be along the lines of self-derogative humor.

To be on topic though ...

The Spherical harmonics $Y_{\ell m}(\theta,\varphi)$ are eigenfunctions of the angular part of the negative of the Laplace operator $-\nabla^2$, i.e., in spherical coordinates
$$
-\nabla^2 = -\partial_r^2 - \frac 2r \partial_r + \frac 1{r^2} \hat \Lambda,
$$
where $\hat \Lambda Y_{\ell m}(\theta,\varphi) = \ell(\ell+1) Y_{\ell m}(\theta,\varphi)$. Being constructed from individual Sturm-Liouville problems in the $\theta$ and $\varphi$ directions, it follows that any function on $\mathbb R^3$ may be written on the form
$$
f(\vec r) = \sum_{\ell = 0}^\infty \sum_{m = -\ell}^\ell f_{\ell m}(r) Y_{\ell m}(\theta, \varphi)
$$
and that the $Y_{\ell m}$ are linearly independent so the expansion is unique. (From now on, let me just write $\sum_{\ell, m}$ for the sums to save typing.)

So, armed with this knowledge, we can diagonalise any linear problem involving the Laplace operator, in particular on the form of Poisson's equation
$$
-\nabla^2 u(\vec r) = \rho(\vec r)
$$
where I have used convenient constant normalisation. As both $u$ and $\rho$ are functions of position, they can both be written in terms of spherical coordinates and therefore expanded in terms of the spherical harmonics as described above with expansion coefficients $u_{\ell m}(r)$ and $\rho_{\ell m}(r)$, respectively. (Note that the spherical harmonics form the expansion basis for functions on a sphere. We are therefore expanding each function on the separate spheres of radius $r$, which is why the expansion coefficients depend on $r$.)

Inserting the expansion into Poisson's equation directly leads to
$$
\sum_{\ell, m} \left[ - u_{\ell m}''(r) - \frac 2r u_{\ell m}'(r) + \frac 1{r^2} \ell(\ell + 1) u_{\ell m}(r)\right] Y_{\ell m}(\theta, \varphi) = \sum_{\ell, m} \rho_{\ell m}(r) Y_{\ell m}(\theta, \varphi).
$$
Now, since the spherical harmonics are linearly independent, it directly follows that
$$
- u_{\ell m}''(r) - \frac 2r u_{\ell m}'(r) + \frac 1{r^2} \ell(\ell + 1) u_{\ell m}(r) - \rho_{\ell m}(r) = 0,
$$
which is an ordinary inhomogeneous differential equation of second order that can be solved. Boundary conditions are typically found from $u|_{r = 0}$ being regular (although in actuality this will follow from the differential equation itself in coordinates other than spherical) and $\lim_{r \to \infty} u = 0$. Solving for each $u_{\ell m}(r)$ solves Poisson's equation for an arbitrary source term $\rho(\vec r)$.

Note that the functions $\rho_{\ell m}$ may be found through the inner product (restricting to real functions for simplicity)
$$
\langle f, g \rangle = \int_{\theta = 0}^\pi \int_{\varphi = 0}^{2\pi} f(\theta,\varphi) g(\theta,\varphi) \sin(\theta) d\theta \, d\varphi
$$
as
$$
\rho_{\ell m}(r) = \frac{\langle Y_{\ell m},\rho\rangle }{\langle Y_{\ell m}, Y_{\ell m}\rangle},
$$
where I am leaving the normalisation term $\langle Y_{\ell m}, Y_{\ell m}\rangle$ as written since there are a relatively large number of different normalisation conventions for spherical harmonics and there will always be someone complaining if you pick a particular one ...

Regardless of normalisation, the numerator is on the form
$$
\langle Y_{\ell m},\rho\rangle = \int_{\theta = 0}^\pi \int_{\varphi = 0}^{2\pi} Y_{\ell m}(\theta,\varphi) \rho(r,\theta,\varphi) \sin(\theta) d\theta \, d\varphi,
$$
making it clear that $\rho_{\ell m}$ is indeed a function of $r$.

Now, looking at a source distribution such that $\rho(r) = 0$ for all $r > R$ for some $R$, the differential equations become homogeneous
$$
- u_{\ell m}''(r) - \frac 2r u_{\ell m}'(r) + \frac 1{r^2} \ell(\ell + 1) u_{\ell m}(r)= 0
$$
for $r > R$. This has the general solution
$$
u_{\ell m}(r) = \frac{A_{\ell m}}{r^{\ell+1}} + B_{\ell m} r^{\ell}.
$$
With the boundary condition at infinity, this leads to $B_{\ell m}= 0$ and therefore $u_{\ell m}(r) = A_{\ell m}/r^{\ell + 1}$ outside of the source distribution. This leads to the multipole expansion
$$
u = \sum_{\ell, m} \frac{A_{\ell m}}{r^{\ell + 1}} Y_{\ell m}(\theta,\varphi)
$$
for $r > R$.

(I had planned to go a bit further but it is getting late ...)

Edit: Minor typos.

 

[Orodruin]

Oh wait, this is the introductory physics homework forum …

 

[Delta2]

Oh wait, this is the introductory physics homework forum …

Well this "self" potential energy is kind of a bit advanced concept since it is not cover in introductory physics book, where they usually talk about potential energy of a body in the field of another body. However I decided to put it here since the solution seemed to me a bit easy and straightforward, without the use of heavy math.

I didn't expect of course that the thread would evolve to spherical harmonics and multipole expansion of a field which are of course advanced physics theory.

 

[vanhees71]

Really? I'd also say that spherical harmonics are at a higher level, but to understand the energy in the gravitational field (Newtonian theory of gravitation of course) of a mass distribution is not so difficult.

You start with a gedanken experiment, asking what energy it takes (or in this case how much energy you get out) when bringing a set of $N$ point masses from infinity into a given static configuration, i.e., particle number $i$ to position $s\vec{x}_i$ given Newton's law for the gravitational force between two point particles:
$$\vec{F}_{21}=-\gamma m_1 m_2 \frac{\vec{x}_2-\vec{x}_1}{|\vec{x}_2-\vec{x}_1|}$$
is the force on particle 2 at presence of particle 1.

For the following it's good to know that this interaction force has a potential, i.e., it's easy to show that
$$\vec{F}_{21}=-\vec{\nabla} \left [-\frac{\gamma m_1 m_2}{|\vec{x}_2-\vec{x}_1|} \right].$$

Now start with bringing particle 1 to $\vec{x}_1$. Since there are no other particles present yet, there's no force acting on it and thus you neither need nor gain energy.

Now bring particle 2 at its position. This is very much easier, because there's the attractive gravitational force due to the presence of particle 1, i.e., you get energy out. Since we count the energy contained in the mass configuration that means this configuration looses energy, i.e., you count this energy as negative:
$$E_{2}=-\frac{\gamma m_1 m_2}{|\vec{x}_1-\vec{x}_2|}.$$
Now you bring particle 3 to its position. Now the attractive gravitational force from both particles 1 and 2 act on it, and according to Newton's theory there is no generic 3-body force involved in gravity, i.e., the total force is simply the sum of the two pair forces. This implies that for the three particles the total energy is
$$E_3=E_2 - \frac{\gamma m_1 m_3}{|\vec{x}_1-\vec{x}_3|} - \frac{\gamma m_2 m_3}{|\vec{x}_2-\vec{x}_3|} = -\gamma \sum_{j<k} \frac{m_j m_k}{|\vec{x}_j-\vec{x}_k|} = -\frac{\gamma}{2} \sum_{j \neq k} \frac{m_j m_k}{|\vec{x}_j-\vec{x}_k|},$$
and here you see the factor $1/2$ which compensates for counting the interaction energies between all possible pairs twice. Here $j$ and $k$ each run from $1$ to $3$, but in the sum you leave out the divergent "self-energy contributions" by contraining the sum to $j \neq k$.

Of course the same kind of sum you have for all $N$ particles, i.e., the $j,k$ each run from 1 to $N$ but leaving out the undefined contributions for $j=k$.

Now go to the continuum limit and think of it as composed of very many infinitesimal volume elements $\mathrm{d}^3 x$. Then instead of the sum you get the integral
$$E_{\text{grav}}=-\frac{\gamma}{2} \int_{\mathbb{R}^3} \mathrm{d}^3 x_1 \int_{\mathbb{R}^3} \mathrm{d}^3 x_2 \frac{\rho(\vec{x}_1) \rho(\vec{x}_2)}{|\vec{x}_1-\vec{x}_2|}.$$
Now tacitly I omitted the important assumption that I don't consider the self-energies, i.e., I tacitly didn't take into account that in the sum over the discrete point particles I've put $j \neq k$. Of course in the latter model I'd get an infinite self-energy due to the singularities of the gravitational interaction potential. In the continuum limit these self-energy contributions are finite, because one integrates out the singularity. So in the continuum model the gravitational field-energy also contains these self-energy contributions.

Often it is simpler to bring this integral in other forms. You can introduce the gravitational field as
$$\vec{g}(\vec{x})=-\gamma \int_{\mathbb{R}^3} \mathrm{d}^3 x' \rho(\vec{x}') \frac{\vec{x}-\vec{x}'}{|\vec{x}-\vec{x}'|^3},$$
which you also can derive from a potential $\Phi$,
$$\vec{g}(\vec{x})=-\vec{\nabla} \Phi \quad \text{with} \quad \Phi(\vec{x})=-\gamma \int_{\mathbb{R}^3} \mathrm{d}^3 x' \frac{\rho(\vec{x}')}{|\vec{x}-\vec{x}'|}.$$
Then the gravitational-field energy obviously can be written as
$$E_{\text{grav}}=\frac{1}{2} \int_{\mathbb{R}^3} \mathrm{d}^3 x \rho(\vec{x}) \Phi(\vec{x}).$$
So often it's simpler to first calculate the potential $\Phi$ and then evaluate the field energy using this formula.

It's a good exercise to calculate $\Phi$ and $E_{\text{grav}}$ for a homogeneous spherical mass distribution,
$$\rho(\vec{x})=\rho_0 \Theta(a-|\vec{x}|).$$

 

[PhDeezNutz]

@Delta2

There is a very nice discussion from slides 6 to 9: https://home.iitk.ac.in/~akjha/PHY103_Notes_HW_Solutions/PHY103_Lec_6.pdf (albeit in an analogous discussion of electrostatics) That talks about the factor of $\frac{1}{2}$ and summing interaction terms while not double counting. It generalizes from the discrete case to the continuous case.

I'm going to try and derive an analogous "field squared" equation for gravitational potentials. Although their example is for a spherical shell I think we can apply it to a solid sphere. (Others feel free to correct me if I am wrong).

$W_{grav} = \frac{1}{2} \int_{all \text{ } space} \vec{g} \cdot \vec{g} \, d \tau = \frac{1}{2} \int_{inside} \vec{g}_{inside} \cdot \vec{g}_{inside} \, d \tau + \frac{1}{2} \int_{outside} \vec{g}_{outside} \cdot \vec{g}_{outside} \, d \tau$

$\vec{g}_{inside} = -\frac{GMr}{R^3}$ where $r \lt R$ (Can be derived through the integral form of Gauss' Law of Gravitation)
$\vec{g}_{outside} = -\frac{GM}{r^2}$ where $r \gt R$

$\vec{g}_{inside} \cdot \vec{g}_{inside} = \frac{G^2 M^2 r^2 }{R^6}$

$\vec{g}_{outside} \cdot \vec{g}_{outside} = \frac{G^2 M^2}{r^4}$

Dividing the integral into 2 parts one for the inside and one for the outside while realizing that this is a spherically symmetric distribution/field so $d \tau = 4 \pi r^2 dr$. We integrate $\vec{g}_{inside} \cdot \vec{g}_{inside}$ from $0$ to $R$ and $\vec{g}_{outside} \cdot \vec{g}_{outside}$ from $R$ to $\infty$

Here we go

$W_1 = \frac{4 \pi G^2 M^2}{R^6} \int_{r=0}^{4} r^4 \,dr = \frac{4 \pi G^2 M^2}{5 R}$

$W_2 = 4 \pi G^2 M^2 \int_{r=R}^{\infty} \frac{1}{r^2} \,dr = \frac{4 \pi G^2 M^2}{R}$

Adding them together $W_{total} = W_{1} + W_{2} = \frac{24 \pi G^2 M^2}{5 R}$

Or something like that. I might be completely wrong.

 

[PhDeezNutz]

Just realized @Orodruin used a similar approach before me. I must have skimmed past it. Seems like we have slightly different answers. Gonna look at mine more closely.

 

[Orodruin]

Just realized @Orodruin used a similar approach before me. I must have skimmed past it. Seems like we have slightly different answers. Gonna look at mine more closely.

You are missing the crucial differences between the gravitational and electrostatic cases in terms of the sign and the constants appearing in Gauss’ law for each scenario when deriving the energy density. I discussed this towards the end of my post.

Edit:
*: Compare this with the electrostatic potential energy stored in an electric field
ue(x→)=12ε0E→(x→)2.
The constants are a bit different because of how ε0 and constants appear in Gauss' law, but the essential argument is exactly the same. (Just replace g→→E→ and −1/4πG→ε0.)

 

[vanhees71]

Let's check it with the other method using the potential. It's easier to work it out with differential equations but you can also do it by using the integrals given in my previous posting. The equation for the gravitational potential is
$$\Delta \Phi=4 \pi \gamma \rho.$$
With $\rho=\rho_0 \Theta(a-|\vec{x}|)$ and in spherical coordinates, using the spherical symmetry of the problem by making the ansatz $\Phi(\vec{x}) \equiv \Phi(r)$, where $r=|\vec{x}|$, you get
$$\Delta \Phi(r)=\frac{1}{r} \partial_r^2 (r \Phi)=4 \pi \gamma \rho(r).$$
For $r<a$ we have $\rho(r)=\rho_0$ and thus by integrating you get
$$\partial_r (r \Phi)=2 \pi \gamma \rho_0 r^2 + A_< \; \Rightarrow \; r \Phi=\frac{2 \pi}{3} \gamma \rho_0 r^3 + A_< r + B_> \; \Rightarrow \; \Phi=\frac{2 \pi}{3} \gamma \rho_0 r^2 + A_<+\frac{B_<}{r}.$$
Since there's no singularity at the origin, we must have $B_<=0$, while $A_<$ is an arbitrary integration constant.

For $r>a$ you get in the same way, because of $\rho=0$,
$$\Phi(r)=\frac{B_>}{r},$$
where I've set the other constant to $0$, because we define the potential such that it goes to 0 at infinity (that's just choosing an arbitrary physically irrelevant additive constant to the total energy).

Now both the potential as well as the gravitational field must be continuous at $r=a$, i.e.,
$$\Phi(a-0^+)=\Phi(a+0^+) \; \Rightarrow \; \frac{2 \pi}{3} \gamma \rho_0 a^2 + A_<=\frac{B_>}{a}$$
and
$$\Phi'(a-0^+)=\Phi'(a-0^+) \; \Rightarrow \; \frac{4 \pi}{3} \gamma \rho_0 a = -\frac{B_>}{a^2}.$$
So we have
$$B_>=-\frac{4 \pi}{3} \gamma \rho_0 a^3=-M \gamma,$$
where $M$ is the total mass of the sphere and from this
$$A_<=\frac{B_>}{a}-\frac{2 \pi}{3} \gamma \rho^0 a^2 = -2 \pi \gamma \rho_0^2 a^2.$$
So we have
$$\Phi(r)=\begin{cases} \frac{2 \pi}{3} \gamma \rho_0 (r^2-3 a^2) & \text{for} \quad r<a \\
-\frac{4 \pi \gamma \rho_0 a^3}{3 r}=-\frac{M}{r} &\text{for} \quad r>a. \end{cases}$$
The total energy of the gravitational field thus is
$$E_{\text{grav}}=\frac{1}{2} \int_{\mathbb{R}^3} \mathrm{d}^3 x \rho(r) \Phi(r) = \frac{1}{2} 4 \pi \int_{0}^a \mathrm{d} r r^2 \rho_0^2 (r^2-3a^2)=-\frac{16}{15} \gamma \pi^2 \rho_0 a^5=-\frac{3}{5a} \gamma M^2.$$
That's different from the result above.

Now let's see how to express the field energy in terms of $\vec{g}=-\vec{\nabla} \Phi$. Note that you have
$$\vec{\nabla} \cdot \vec{g}=-4 \pi \rho$$
and thus
$$E=\frac{1}{2} \int_{\mathbb{R}^3} \mathrm{d}^3 x \rho \Phi=-\frac{1}{8 \pi \gamma} \int_{\mathbb{R}} \mathrm{d}^3 x (\vec{\nabla} \cdot \vec{g}) \Phi = +\frac{1}{8 \pi \gamma} \int_{\mathbb{R}} \mathrm{d}^3 x \vec{g} \cdot \vec{\nabla} \Phi = -\frac{1}{8 \pi \gamma} \int_{\mathbb{R}} \vec{g}^2.$$
For our case, obviously you have $\vec{g}=g(r) \vec{e}_r$ with
$$g(r)=-\Phi'(r) = \begin{cases} -\frac{4 \pi}{3} \gamma \rho_0 r &\text{for} \quad r<a \\
-\frac{4 \pi a^3}{3r} \gamma \rho_0 & \text{for} \quad r>a. \end{cases}$$
Plugging this into the said integral you get the same as above with the potential.

 

[Delta2]

Well @vanhees71 I am Greek and you are German and I like it that you use the Greek letter "gamma" $\gamma$ for the gravitational constant G , however it looks like it has something to do with the gamma of special relativity...