Further information about the possible orbits in the Schwarzschild geometry can be obtained by analyzing the potential more carefully. Using the same conventions as before, we have \begin{equation} \frac12 \dot{r}^2 = \frac{e^2-1}{2} - V \label{rdot} \end{equation} with \begin{equation} V = \frac12\left(1-\frac{2m}{r}\right) \left(1+\frac{\ell^2}{r^2}\right) - \frac12 \end{equation} Differentiating $V$ with respect to $r$ yields \begin{equation} V' = \frac{dV}{dr} = \frac{mr^2-\ell^2r+3m\ell^2}{r^4} \label{Vprime} \end{equation} so that $V'=0$ when \begin{equation} \frac{r}{m} = \frac{\ell}{2m} \left( \frac{\ell}{m}\pm\sqrt{\frac{\ell^2}{m^2}-12} \right) \end{equation}
We have argued graphically that circular orbits can only occur when $V'=0$, and this can only occur if \begin{equation} \ell^2 \ge 12 m^2 \end{equation} so that the roots are real. Circular orbits do not exist for all values of the parameters $\ell$ and $e$! But why must $V'$ vanish for circular orbits?
Return to expression ($\ref{rdot}$) for $\dot{r}$, and differentiate with respect to $r$, leading to \begin{equation} \dot{r}\ddot{r} = -V' \dot{r} \label{rddot0} \end{equation} We would like to conclude that $\ddot{r}=0$ if and only if $V'=0$, but this appears to require us to assume $\dot{r}\ne0$. One could argue that well-behaved solutions must depend continuously on the parameters, so that it is legitimate to divide ($\ref{rddot0}$) by $\dot{r}$. A better argument is to begin by solving ($\ref{rdot}$) for $\dot{r}$, yielding \begin{equation} \dot{r} = \sqrt{e^2-1-2V} \end{equation} so that \begin{equation} \ddot{r} = \frac{-V'\dot{r}}{\sqrt{e^2-1-2V}} = \frac{-V'\dot{r}}{\dot{r}} = -V' \end{equation} and now the last equality can be justified in the usual way by computing the derivative using limits.
What's going on here? The use of symmetry to find constants of the motion effectively turns the second-order geodesic equation into a first-order differential equation, and this process can introduce spurious solutions. We must in principle check that our claimed solutions actually solve the full second-order equation. We will postpone that computation for now, as it first requires computing the connection in the Schwarzschild geometry, but the ($r$-component of the) full geodesic equation does indeed imply that \begin{equation} \ddot{r} = -V' = - \frac{mr^2-\ell^2r+3m\ell^2}{r^4} \end{equation} for all geodesics, including circular geodesics. 1)
Now that we have established that $V'$ vanishes on circular orbits, it is straightforward to use ($\ref{Vprime}$) to express $\ell$ in terms of $r$, yielding \begin{equation} \frac{\ell^2}{r^2} = \frac{\frac{m}{r}}{1-\frac{3m}{r}} \label{l2eq} \end{equation} and ($\ref{rdot}$) can then be used to express $e$ in terms of $r$, yielding \begin{equation} e^2 = \frac{\left(1-\frac{2m}{r}\right)^2}{1-\frac{3m}{r}} \label{e2eq} \end{equation} Not surprisingly, specifying the radius of a circular orbit completely determines the orbiting object's energy and angular momentum (per unit mass).
We can use ($\ref{l2eq}$) and ($\ref{e2eq}$) to compute the apparent angular velocity $\Omega$ of an orbiting object as seen by an observer far away, since \begin{equation} \Omega^2 = \frac{d\phi^2}{dt^2} = \frac{\dot{\phi}^2}{\dot{t}^2} = \frac{\ell^2/r^4}{e^2/\left(1-\frac{2m}{r}\right)^2} = \frac{m}{r^3} \end{equation} Thus, \begin{equation} m^1 = \Omega^2 r^3 \end{equation} which is the same “1–2–3” law as in Kepler's third law (for circular orbits).