newtonian gravity - 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?

Answer

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.