# Is it possible to prove that planets should be approximately spherical using the calculus of variations?

Is it possible to use the Lagrangian formalism involving physical terms to answer the question of why all planets are approximately spherical?

Let's assume that a planet is 'born' when lots of particles of uniform density are piled together, and constrained only by the Newtonian gravity of the $N$ particles with no external gravity.

My question is: Is it possible to use the variational calculus starting from this situation to prove that planets should be spherical?

• @Jim: Although he uses the word "approximately", I think what the OP is after is a proof that a sphere is indeed the minimum energy shape of an object being acted on by its own gravity, ie, that a perturbation to a nonspherical shape along any conformational degree of freedom will be energetically unfavorable. Apr 9, 2014 at 16:27

I) Let us work in units where $$G=1=\rho$$.

In this answer, we will additionally assume that the massive object is star-shaped. We can then use spherical coordinates $$(r,\theta,\varphi)$$. The surface profile is then given as

$${\bf r}~=~r{\bf n}~=~f({\bf n}){\bf n}, \qquad r~=~f({\bf n})~=~\sum_{\ell m} c_{\ell m} Y_{\ell m}({\bf n})~\geq~0 ,$$ $$c_{\ell m}~=~\int_{S^2}\!d^2n ~Y^{\ast}_{\ell m}({\bf n})f({\bf n}),\tag{1}$$

where $${\bf n}\in S^2$$ is a unit vector. In eq. (1) we have used DumpsterDoofus' idea to expand in spherical harmonics. The volume

$$V[f]~:=~ \int_V \! d^3r ~=~\frac{1}{3}\int_{S^2}\!d^2n ~ f({\bf n})^3 .\tag{2}$$

is a functional of the surface profile $$f({\bf n})$$.

$$\delta V ~=~\int_{S^2}\!d^2n ~ f({\bf n})^2\delta f({\bf n}),\qquad \delta^2 V ~=~2\int_{S^2}\!d^2n ~ f({\bf n})\delta f({\bf n})^2.\tag{3}$$

The volume is kept fixed by a constraint

$$V[f]~=~V_0.\tag{4}$$

Potential:

$$-\Phi({\bf r})~:=~\int_{V} \! \frac{d^3r^{\prime}}{|{\bf r}-{\bf r}^{\prime}|} ~=~\int_{S^2}\!d^2n^{\prime}\int_0^{f({\bf n}^{\prime})}\! \frac{r^{\prime 2}dr^{\prime}}{|{\bf r}-r^{\prime}{\bf n}^{\prime}|} .\tag{5}$$

$$-\delta \Phi({\bf r}) ~=~\int_{S^2}\!d^2n^{\prime} \frac{f({\bf n}^{\prime})^2\delta f({\bf n}^{\prime})}{|{\bf r}-f({\bf n}^{\prime}){\bf n}^{\prime}|}.\tag{6}$$

Field:

$$-{\bf g}({\bf r})~:=~ \nabla \Phi({\bf r}) ~=~\int_{V} \!d^3r^{\prime} \frac{{\bf r}-{\bf r}^{\prime}}{|{\bf r}-{\bf r}^{\prime}|^3} .\tag{7}$$

\begin{align}-g_r({\bf r})~=~&\frac{\partial\Phi({\bf r})}{\partial r} ~=~{\bf n}\cdot \nabla \Phi({\bf r}) ~=~\int_{V} \!d^3r^{\prime} \frac{r-{\bf n}\cdot{\bf r}^{\prime}}{|{\bf r}-{\bf r}^{\prime}|^3}\cr ~=~&\int_{S^2}\!d^2n^{\prime}\int_0^{f({\bf n}^{\prime})}\!r^{\prime 2}dr^{\prime} \frac{r-{\bf n}\cdot{\bf n}^{\prime}r^{\prime}}{|{\bf r}-r^{\prime}{\bf n}^{\prime}|^3}.\end{align}\tag{8}

II) Potential energy:

\begin{align} U[f]~:=~& -\iint_{V\times V} \! \frac{d^3r~ d^3r^{\prime}}{2|{\bf r}-{\bf r}^{\prime}|} ~=~ \frac{1}{2}\int_V \! d^3r ~\Phi({\bf r}) \cr ~=~& -\int_{S^2}\!d^2n\int_{S^2}\!d^2n^{\prime} \int_0^{f({\bf n})}\int_0^{f({\bf n}^{\prime})} \frac{r^2dr~ r^{\prime 2}dr^{\prime}}{2|r{\bf n}-r^{\prime}{\bf n}^{\prime}|}.\end{align}\tag{9}

$$\delta U ~=~\int_{S^2}\!d^2n ~ f({\bf n})^2\delta f({\bf n})~\Phi(f({\bf n}){\bf n}).\tag{10}$$

\begin{align}\delta^2 U~=~&-\int_{S^2}\!d^2n\int_{S^2}\!d^2n^{\prime}~ \frac{f({\bf n})^2\delta f({\bf n})~ f({\bf n}^{\prime})^2\delta f({\bf n}^{\prime})}{|f({\bf n}){\bf n}-f({\bf n}^{\prime}){\bf n}^{\prime}|}\cr &-\int_{S^2}\!d^2n~ f({\bf n})^2 ~\delta f({\bf n})^2g_r(f({\bf n}){\bf n})\cr &+2\int_{S^2}\!d^2n~f({\bf n})\delta f({\bf n})^2 \Phi(f({\bf n}){\bf n}).\end{align}\tag{11}

III) To treat the constraint (4) we use Lagrange multiplier method. We should minimize the functional

$$E[f]~:=~U[f] +\lambda (V_0-V[f]).\tag{12}$$

The first variation is $$\delta E ~=~\int_{S^2}\!d^2n ~ f({\bf n})^2\delta f({\bf n}) \{\Phi(f({\bf n}){\bf n})-\lambda\}.\tag{13}$$

We conclude that

The surface potential $$\Phi(f({\bf n}){\bf n})=\lambda<0$$ of a stationary shape $$f$$ is a constant, i.e. independent of the unit vector $${\bf n}$$.

The second variation is

\begin{align}\delta^2 E~=~&-\int_{S^2}\!d^2n\int_{S^2}\!d^2n^{\prime}~ \frac{f({\bf n})^2\delta f({\bf n})~ f({\bf n}^{\prime})^2\delta f({\bf n}^{\prime})}{|f({\bf n}){\bf n}-f({\bf n}^{\prime}){\bf n}^{\prime}|}\cr &-\int_{S^2}\!d^2n~ f({\bf n})^2 \delta f({\bf n})^2g_r(f({\bf n}){\bf n})\cr &+2\int_{S^2}\!d^2n~f({\bf n})\delta f({\bf n})^2 \underbrace{\{\Phi(f({\bf n}){\bf n})-\lambda\}}_{=0}.\end{align}\tag{14}

For systematic reasons, the second variations (14) should be consistent with the constraint (4) to first order.

IV) We finally consider a ball $$f({\bf n})=R$$. We then have

$$V~=~\frac{4\pi}{3}R^3, \qquad {\bf g}({\bf r})~=~-\frac{4\pi}{3}{\bf r}, \qquad \Phi({\bf r})~=~2\pi\left(\frac{r^2}{3}-R^2\right),$$ $$\Phi(R)~=~-\frac{4\pi}{3}R^2, \qquad U~=~-\frac{16\pi^2}{15}R^5.\tag{15}$$

We want to show that the ball is a stable stationary shape. We have to show that the second variation (14) is semipositive definite. Eq. (3) simplifies to

$$\frac{\delta V}{R^3} ~=~\int_{S^2}\!d^2n ~\delta f({\bf n}) ~=~\sqrt{4\pi}\delta c_{00},\tag{16}$$

$$\frac{\delta^2 V}{2R^3}~=~\int_{S^2}\!d^2n ~|\delta f({\bf n})|^2 ~=~\sum_{\ell m}|\delta c_{\ell m}|^2.\tag{17}$$

Since the constraint (4) should be maintained to first (but not necessarily second) order, we must demand that the zeromode vanishes

$$\delta c_{00}~=~0.\tag{18}$$

The second variation (14) simplifies to (see e.g. Ref. 1)

\begin{align}\frac{\delta^2 E}{R^3} ~=~&-\int_{S^2}\!d^2n \int_{S^2}\!d^2n^{\prime}~\frac{\delta f({\bf n})~\delta f({\bf n}^{\prime})^{\ast}}{|{\bf n}-{\bf n}^{\prime}|} -\frac{g_r(R)}{R}\int_{S^2}\!d^2n~\delta f({\bf n})^2 \cr ~=~&-\int_{S^2}\!d^2n\int_{S^2}\!d^2n^{\prime}~\delta f({\bf n})~\delta f({\bf n}^{\prime})^{\ast}\sum_{\ell m} \frac{4\pi}{2\ell+1}Y^{\ast}_{\ell m}({\bf n})Y_{\ell m}({\bf n}^{\prime})\cr &+\frac{4\pi}{3}\int_{S^2}\!d^2n~|\delta f({\bf n})|^2 \cr ~=~&4\pi\sum_{\ell m}\left(\frac{1}{3}-\frac{1}{2\ell+1}\right)|\delta c_{\ell m}|^2 ~\geq~0,\end{align}\tag{19}

which is non-negative since the zeromode (18) is absent. So the ball is a stable stationary shape.

References:

1. J.D. Jackson, Classical Electrodynamics; chapter 3.

I interpreted the question as asking whether it's possible to show that a sphere is the minimum energy shape of an object being acted on by its own gravity. I attempted to do this using spherical harmonics, but got stuck part of the way there; I'm posting this anyways just in case someone can figure out how to complete the last bit.

The gravitational self-energy of a matter distribution $\rho(\mathbf{r})$ which gives rise to a gravitational potential $V(\mathbf{r})$ is given by $$E=\frac{1}{2}\left\langle V,\rho\right\rangle=\frac{1}{2}\left\langle \nabla^{-2}\rho,\rho\right\rangle$$ where $\langle,\rangle$ represents the inner product and $\nabla^{-2}$ is the inverse Laplacian. If there is a change in the matter distribution $\rho=\rho_0+\delta\rho$, then the self-energy changes to $$E=\frac{1}{2}\left\langle \nabla^{-2}(\rho_0+\delta\rho),\rho_0+\delta\rho\right\rangle=\frac{1}{2}\left\langle \nabla^{-2}\rho_0,\rho_0\right\rangle+\left\langle \nabla^{-2}\rho_0,\delta\rho\right\rangle+\frac{1}{2}\left\langle \nabla^{-2}\delta\rho,\delta\rho\right\rangle \\ =E_0+\delta E+O(\delta^2)$$ where the self-adjointness of the Laplacian was used in the second to last step.

Thus it suffices to determine the sign of $\delta E=\left\langle V_0,\delta\rho\right\rangle$. Since $V_0$ is just the gravitational potential of a uniform-density sphere (which has simple closed form) this becomes tractable.

Let the matter density be $\rho_d$ and assume the planet is incompressible. The radius of the planet $R(\Omega)$ can be expanded in the real harmonics as $$R(\Omega)=R_0+\sum_{L=0}^\infty\sum_{m=-L}^L\delta C_{Lm}Y_{Lm}(\Omega)=R_0+\delta R(\Omega).$$ To conserve volume, note that $$V=\frac{4}{3}\pi R_0^3=\frac{1}{3}\int R(\Omega)^3\,d\Omega=\frac{1}{3}\int \left(R_0+\sum_{L=0}^\infty\sum_{m=-L}^L\delta C_{Lm}Y_{Lm}(\Omega)\right)^3\,d\Omega \\ =\frac{4}{3}\pi R_0^3+2\sqrt{\pi}R_0^2\delta C_{00}+R_0\sum_{L=0}^\infty\sum_{m=-L}^L\delta C_{Lm}^2+O(\delta^3) \\ \Rightarrow \delta C_{00}=-\frac{1}{2\sqrt{\pi}R_0}\sum_{L=1}^\infty\sum_{m=-L}^L\delta C_{Lm}^2+O(\delta^3).$$ Since the volume conservation term $\delta C_{00}$ scales as $\delta^2$ let's rewrite it as $\delta^2C_{00}$, giving $$R(\Omega)=R_0+\sum_{L=1}^\infty\sum_{m=-L}^L\delta C_{Lm}Y_{Lm}(\Omega)+Y_{00}\delta^2C_{00}=R_0+\delta R_1+\delta^2 R_2.$$ With a bit of visualization one can see that $$\left\langle V_0,\delta\rho\right\rangle=\int d\Omega R(\Omega)^2\int_{R_0}^{R(\Omega)}V_0(r)\rho_d\,dr \\ =\int d\Omega R(\Omega)^2\left[-\frac{G M R_1 \rho_d}{R_0}\delta+\frac{GM(R_1^2-2R_0R_2)\rho_d}{2R_0^2}\delta^2+\:...\right] \\ =\int d\Omega\left[-GMR_0\delta R_1\rho_d-\frac{1}{2}GM\rho_d\left(3(\delta R_1)^2+2R_0\delta^2R_2\right)\:...\right]$$ where in the last step terms in the integral have been arranged according to their order in $\delta$.

The first term vanishes due to symmetry of the real harmonics, and the second term simplifies due to orthonormality to become $$=-\frac{1}{2}GM\rho_d\left[3\sum_{L=1}^\infty\sum_{m=-L}^L\delta C_{Lm}^2-2\sum_{L=1}^\infty\sum_{m=-L}^L\delta C_{Lm}^2\right] \\ =-\frac{1}{2}GM\rho_d\sum_{L=1}^\infty\sum_{m=-L}^L\delta C_{Lm}^2<0.$$ Unfortunately, this predicts that the energy change is always favorable to deformation, rather than always unfavorable! So I must have an arithmetic error somewhere, but I am too tired to find out where it is (the units at least work out correctly to Joules). However, this gives a general idea of how one can proceed down this line of proof.