Solving 3D Kepler Problem substitution goes wrong I'm trying to arrive at the effective potential equation in Kepler Problem using Routh reduction method. We can procede in two ways, either using polar coordinates in the plane where the orbit happens or using spherical coordinates. I'm having trouble with this last one. I'm gonna follow steps taken in this Wikipedia page. Recalling,
\begin{gather*}
\mathcal{L}(r, \dot{r}, \theta, \dot{\theta}, \dot{\phi}) = \frac{1}{2}m \left(\dot{r}^2 + r^2 \dot{\theta}^2 + r^2 \sin^2 (\theta) \dot{\phi}^2\right) - V(r).\tag{1}
\end{gather*}
Because $\phi$ is cyclic, its momentum conjugate is conserved
\begin{gather*}
p_{\phi} = \frac{\partial \mathcal{L}}{\partial \dot{\phi}} = mr^2 \sin^2 (\theta) \dot{\phi} = L_z = cte.\tag{2}
\end{gather*}
Now, consider the Routhian
\begin{gather*}
\mathcal{R}(r, \dot{r}, \theta, \dot{\theta}) = \frac{1}{2} \frac{p_{\phi}^2}{mr^2\sin^2 (\theta)} - \frac{1}{2}m \left(\dot{r}^2 + r^2 \dot{\theta}^2 \right) + V(r).\tag{3}
\end{gather*}
Now consider $\theta$ Lagrange equation, which is equivalent to the conservation of the modulus of the momentum
\begin{align*}
m(2r\dot{r}\dot{\theta} + r^2\ddot{\theta}) &= \frac{p_{\phi}^2\cos(\theta)}{mr^2\sin^3(\theta)} \tag{4} \\
m^2r^4\dot{\theta}^2 + \frac{p_{\phi}^2}{\sin^2(\theta)} &= L^2 = cte.\tag{5}
\end{align*}
However, if I substitute in the Routhian, it does not work properly
\begin{gather*}
\mathcal{R}(r, \dot{r}, \dot{\theta}) = \frac{1}{2} \frac{L^2}{mr^2} - \frac{1}{2}m \left(\dot{r}^2 + 2r^2 \dot{\theta}^2 \right) + V(r).\tag{6}
\end{gather*}
 A: *

*The underlying reason for OP's flawed argument is, that a premature use of EOMs in the stationary action principle
$$\begin{align} S~=~-&\int\!dt ~R(r,\dot{r};\theta,\dot{\theta}), \cr 
-R(r,\dot{r};\theta,\dot{\theta})~=~&\frac{1}{2}m(\dot{r}^2 +r^2\dot{\theta}^2) -\frac{p_{\phi}^2}{2mr^2\sin^2\theta} -V(r),
 \end{align}\tag{A}$$
invalidates the variational principle, cf. e.g. this Phys.SE post.


*In OP's case the EOM happens to be the integrated EOM for $\theta$:
$$\begin{align}{\rm const}~\approx~\vec{L}^2~=~&(\vec{r}\times \vec{p})^2~=~r^2\vec{p}_{\perp}^2~=~m^2r^2\dot{\vec{r}}_{\perp}^2\cr
~=~&m^2r^4\left(\dot{\theta}^2+\sin^2\theta ~\dot{\phi}^2\right)~=~p_{\theta}^2 + \frac{p_{\phi}^2}{\sin^2\theta}.\end{align}\tag{B}$$
(An on-shell conserved quantity is always the result of some EOM.)


*Specifically, OP plugs EOM (B) into the Routhian (A) to obtain
$$-R(r,\dot{r};\theta,\dot{\theta})~\approx~\frac{1}{2}m(\dot{r}^2 +\color{red}{2}r^2\dot{\theta}^2) -\frac{\color{red}{\vec{L}^2}}{2mr^2} -V(r)\tag{C}$$
on-shell. However the on-shell Routhian (C) cannot be used to deduce the EOMs.
A: Substituting some conserved quantity is not how the Routhian mechanics works. If your Lagrangian is $\mathcal{L}(r,\dot{r},\theta,\dot{\theta},\phi,\dot{\phi})$ (we do not exclude $\phi$ here as  generalization for other possible potential), we always have a Routhian $R(r,\dot{r},\theta,\dot{\theta})$ given by
$$R(r,\dot{r},\theta,\dot{\theta})=p_{\phi}\dot{\phi}-\mathcal{L}$$
which satisfies the equations of motion
\begin{align}
{d \over dt}\bigg({\partial R \over \partial \dot{r}}\bigg) & = {\partial R \over \partial r} \\
{d \over dt}\bigg({\partial R \over \partial \dot{\theta}}\bigg) & = {\partial R \over \partial \theta}
\end{align}
You can also choose some other Routhians like $R(r,\dot{r},\phi,\dot{\phi})=p_{\theta}\dot{\theta}-\mathcal{L}$, and it will satisfy
\begin{align}
{d \over dt}\bigg({\partial R \over \partial \dot{r}}\bigg) & = {\partial R \over \partial r} \\
{d \over dt}\bigg({\partial R \over \partial \dot{\phi}}\bigg) & = {\partial R \over \partial \phi}
\end{align}
So as pointed out in the wiki page, the choice of generalized coordinates to obtaining the corresponding Ruthian is arbitrary. However, some can probably make solving the problem easier. In the example given in your question, since $p_\phi$ is conserved, there will be $2$ equations of motion $R(r,\dot{r},\theta,\dot{\theta})$ have to satisfy with only $2$ degrees of freedom left. This is not the case if we choose $R(r,\dot{r},\phi,\dot{\phi})$. However, you cannot have the Routhian with the substitution for the total angular momentum $L$ since we do not have $p_\phi\dot{\phi}={L^2 \over mr^2}$.
