Bertrand's Theorem

Theorem
Let $U:\mathbb{R}^+\to\mathbb{R}$ be analytic for $r>0$, and let $M>0$ be a nonvanishing angular momentum such that a stable circular orbit exists. If every orbit sufficiently close to the circular orbit is closed, then $U$ is either $k r^2$ or $-k/r$ (for $k>0$) up to an additive constant.

Proof
For simplicity we set $m=1$, so that the effective potential becomes $U_{M}=U+\frac{M^2}{2 r^2}$. Recall that an orbit is closed, if the angle between adjacent apocenters (pericenters) is commensurable with $\pi$, and that this angle, called the apsidial angle, is given by $\Phi=\sqrt{2}\int_{r_{min}}^{r_{max}} \frac{M \mathrm d r}{r^2 \sqrt{E-U_{M}}}$, where $E$ is the energy, and $r_{min},r_{max}$ are solutions to $U_{M}(r)=E$.

Non-perturbative proof
In general $U_{M}$ is not monotonic on $(r_{min},r_{max})$, so a unique inverse $r(U_{M})$ does not exist. However suppose it is possible to construct separate inverse functions $r_{1,2}$ for the intervals $(r_0,r_{min})$ and $(r_0,r_{max})$, where $r_0$ is the minimum of $U_{M}$ on $(r_{min},r_{max})$. Note that the orbit corresponding to $r_0$ is the stable circular orbit, so this will be possible for orbits sufficiently close to it. Now write $\Phi=\int_{U_0}^E F(U_{M}) \frac{\mathrm d U_{eff}}{\sqrt{E-U_{M}}}$, where $F=\sqrt2 M \frac{\mathrm d}{\mathrm d U_{M}}\big[\frac{1}{r_1}-\frac{1}{r_2}\big]$ and $U_0\equiv U_M(r_0)$. This is Abel's integral equation which can be solved for $F(U_M)$, giving:
 * $F(U) = \frac{1}{\pi} \int_{U_0}^{U_M} \frac{\Phi \mathrm d E}{\sqrt{U_{M}-E}}$.

$\Phi$ may have energy dependance, however we require $\Phi(E)=2\pi q(E)$ such that $q:\mathbb{R}\to \mathbb{Q}$ is continuous, and a rational continuous function must be constant. Finally, integrating gives:
 * $\frac{1}{r_1}-\frac{1}{r_2} = \gamma \sqrt{U_{M}-U_0}, \gamma = \frac{\sqrt2 \Phi}{\pi M}$.

Now consider defining the left hand side to be a single function $1/r(U)$, meromorphic with a simple pole at $r=0$, and write $\frac{1}{r} = \gamma \sqrt{U_{M}-U_0} + \Omega(U_0,U)$ where $\Omega$ is an analytical function in $U$ in an open neighbourhood of $U_0$ such that $\Omega(U_0,U_0)=1/r_0$. Explicitly, $\sqrt{U_{M}-U_0}=\sqrt{U+\frac{M^2}{2 r^2}-U_0}=\frac{1}{r} \sqrt{M^2/2 +(U-U_0) r^2}$. We see that the RHS also has a simple pole at $r=0$, however, it also has a branching point at $r=r_0$. To avoid this, the radicand must be the square of an analytic function with a zero at $r_0$. The only possibilities are $U\propto r^2, -1/r^2$

Proof based on asymptotic behaviour
The substitution $x=M/r$ gives $\Phi=\sqrt{2}\int_{x_{min}}^{x_{max}} \frac{ \mathrm d x}{ \sqrt{E-W(x)}}$, where $W(x)\equiv U(M/x)+\frac{1}{2}x^2$. In general, the frequency of oscillations around a stable equilibrium at $x=x_0$ for a particle of mass $m$ in a potential $V$ is given by $\omega = \sqrt{\frac{V(x_0)}{m}}$, which means that $\Phi\approx 2\pi\frac{M}{x_0^2}\sqrt{W(x_0)}=2\pi\sqrt{\frac{U'}{3U'+x_0U''}}$. If $\Phi=const.$, the solutions are $U(r)=k r^\alpha, \alpha\geq-2,\alpha\neq0$ and $U(r)=b\log r$, and therefore $\Phi = 2\pi/\sqrt{\alpha+2}$, with $\alpha=0$ corresponding to the logarithmic solution - which we discard since in that case $\Phi$ is not commensurable with $\pi$.

Now consider $\alpha>0$, so that $U(r)\rightarrow\infty$ as $r\rightarrow\infty$. The substitution $x=y x_{max}$ reduces $\Phi$ to $\Phi=\sqrt{2}\int_{y_{min}}^{1} \frac{ \mathrm d y}{\sqrt{W(1)-W(y)} }$ where $W(y)\equiv \frac{y^2}{2}+\frac{1}{x_{max}^2}U(\frac{M}{y x_{max}})$, so when $E\rightarrow\infty$, $x_{max}\rightarrow\infty$ and $y_{min}\rightarrow 0$. This means that $\lim_{E\rightarrow\infty}\Phi=\pi$. Equating this with the previous result gives $\frac{2\pi}{\sqrt{\alpha+2}}=\pi$, therefore $\alpha=2$.

Now let $U(r)=-k r^{-\beta},0<\beta<2$. A similar approach gives $\lim_{E\rightarrow -0}\Phi=\frac{2\pi}{2-\beta}$, so that $\frac{2\pi}{2+\alpha}=\frac{2\pi}{\sqrt{2+\alpha}}$ from which it follows that $\alpha=-1$.