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 $\Psi=\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{\Psi \mathrm d E}{\sqrt{U_{M}-E}}$.

$\Psi$ may have energy dependance, however we require $\Psi(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 \Psi}{\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$