All spherically symmetric systems contain in principle a conserved vector that appears variously under the names Runge, Laplace and Lenz. Since the explicit form of the vector is only simple for an inverse single power of the radial variable, one often hears that it is a special ‘hidden’ symmetry of the Kepler potential alone. On the contrary, it is straightforward to write down the differential equations determining the vector for any spherically symmetric potential. More importantly, there is a relatively simple solution to the equations for the general case, and the solution is not limited to the Kepler potential.

Begin with the Hamiltonian for an *n*-dimensional spherically symmetric system,

$H=\frac{{p}^{2}}{2m}+V\left(r\right)=E$

where the variables have their usual definitions as vector squares:

${r}^{2}=(\mathbf{r}\xb7\mathbf{r})=\sum _{i=1}^{n}{r}_{i}^{2}\phantom{\rule{4em}{0ex}}{p}^{2}=(\mathbf{p}\xb7\mathbf{p})=\sum _{i=1}^{n}{p}_{i}^{2}$

Allowing the coordinate space to have more than the physical three dimensions will not have any effect upon development of the generalized conserved vector and indicates explicitly that the vector is not dependent upon the dimensionality of space. The dot product of these two vectors

$(\mathbf{r}\xb7\mathbf{p})=\sum _{i=1}^{n}{r}_{i}{p}_{i}=r{p}_{r}$

will figure prominently in the development. The last equality provides a definition of the momentum conjugate to the radial variable, which can be verified by evaluating the Poisson bracket of the two quantities:

$[r,{p}_{r}]=\sum _{i=1}^{n}\frac{\partial r}{\partial {r}_{i}}\frac{\partial}{\partial {p}_{i}}\left[\sum _{j=1}^{n}\frac{{r}_{j}{p}_{j}}{r}\right]=\sum _{i=1}^{n}\frac{{r}_{i}^{2}}{{r}^{2}}=1$

Temporal derivatives of conjugate variables are evaluated in the usual manner of Hamiltonian mechanics,

$\frac{d{r}_{i}}{dt}=\frac{\partial H}{\partial {p}_{i}}=\frac{{p}_{i}}{m}\phantom{\rule{5em}{0ex}}\frac{d{p}_{i}}{dt}=-\frac{\partial H}{\partial {r}_{i}}=-\frac{{r}_{i}}{r}\frac{dV}{dr}$

from which follow the temporal derivatives

$\begin{array}{l}\phantom{\rule{2.2em}{0ex}}\frac{dr}{dt}=\frac{(\mathbf{r}\xb7\mathbf{p})}{mr}\\ \phantom{\rule{2.2em}{0ex}}\frac{dp}{dt}=-\frac{(\mathbf{r}\xb7\mathbf{p})}{rp}\frac{dV}{dr}\\ \frac{d}{dt}(\mathbf{r}\xb7\mathbf{p})=\frac{{p}^{2}}{m}-r\phantom{\rule{.2em}{0ex}}\frac{dV}{dr}\end{array}$

Define the square of the angular momentum by

${L}^{2}={r}^{2}{p}^{2}-(\mathbf{r}\xb7\mathbf{p}{)}^{2}$

which is easily seen to be a constant of the motion using the given temporal derivatives. This definition is equivalent to the square of a cross product in three dimensions, and has the advantage that it holds for any number of spatial dimensions greater than one.

The square of the angular momentum can now be used to replace the momentum variables conjugate to spatial angles in the Hamiltonian in the usual manner,

$H=\frac{{p}_{r}^{2}}{2m}+\frac{{L}^{2}}{2m{r}^{2}}+V\left(r\right)=E$

and the temporal derivatives of the two conjugate variables are then

$\begin{array}{l}\frac{dr}{dt}=\frac{\partial H}{\partial {p}_{r}}=\frac{{p}_{r}}{m}=\frac{(\mathbf{r}\xb7\mathbf{p})}{mr}\\ \frac{d{p}_{r}}{dt}=-\frac{\partial H}{\partial r}=\frac{{L}^{2}}{m{r}^{3}}-\frac{dV}{dr}\end{array}$

Conservation of energy implies that there is mathematically only one independent variable in this system, since one can write

$\begin{array}{l}\phantom{\rule{1.5em}{0ex}}{p}^{2}=2m(E-V)\\ \phantom{\rule{1.5em}{0ex}}{p}_{r}^{2}=2m(E-V)-\frac{{L}^{2}}{{r}^{2}}\\ (\mathbf{r}\xb7\mathbf{p})=\sqrt{2m{r}^{2}(E-V)-{L}^{2}}\end{array}$

and these quantities can be treated as explicit functions of the radial variable and implicit functions of time in differential equations and integrals.

Now introduce a general *n*-dimensional vector of the form

$\mathbf{R}=A\left(r\right)\mathbf{r}-B\left(r\right)\mathbf{p}$

The minus sign within the vector simplifies the subsequent development. This vector will be constrained to be constant by an appropriate choice of the coefficient functions *A* and *B*. Evaluate the temporal derivative of this vector and set it equal to zero:

$\frac{d\mathbf{R}}{dt}=[\frac{(\mathbf{r}\xb7\mathbf{p})}{mr}\frac{dA}{dr}+\frac{B}{r}\frac{dV}{dr}]\mathbf{r}+[\frac{A}{m}-\frac{(\mathbf{r}\xb7\mathbf{p})}{mr}\frac{dB}{dr}]\mathbf{p}=0$

This is equivalent to the coupled equations

$\begin{array}{l}\sqrt{2m{r}^{2}(E-V)-{L}^{2}}\frac{dA}{dr}=-m\frac{dV}{dr}B\\ \sqrt{2m{r}^{2}(E-V)-{L}^{2}}\frac{dB}{dr}=r\phantom{\rule{.2em}{0ex}}A\end{array}$

These equations can be converted into linear second-order equations for each coefficient function by taking another derivative with respect to the radial variable and substituting for known quantities. The resulting equations have regular singularities at infinity and all finite points, indicating that a relatively simple solution is likely to exist. This solution, however, appears more cleanly without the introduction of second-order equations.

First write the square of the Runge vector explicitly,

$\begin{array}{l}{R}^{2}={r}^{2}{A}^{2}-2(\mathbf{r}\xb7\mathbf{p})AB+{p}^{2}{B}^{2}\\ \phantom{{R}^{2}}={r}^{2}{A}^{2}-2\sqrt{2m{r}^{2}(E-V)-{L}^{2}}AB+2m(E-V){B}^{2}\end{array}$

This quadratic equation can be easily solved for either coefficient function:

$\begin{array}{l}A=\frac{1}{{r}^{2}}[\sqrt{2m{r}^{2}(E-V)-{L}^{2}}B\pm \sqrt{{r}^{2}{R}^{2}-{L}^{2}{B}^{2}}\phantom{\rule{.2em}{0ex}}]\\ B=\frac{1}{2m(E-V)}[\sqrt{2m{r}^{2}(E-V)-{L}^{2}}A\pm \sqrt{2m(E-V){R}^{2}-{L}^{2}{A}^{2}}\phantom{\rule{.2em}{0ex}}]\end{array}$

These two results can be substituted in the right-hand sides of the coupled equations to produce separated first-order equations for each coefficient function. The first-order equations look nonlinear at first sight, but both can be solved with simple substitutions. The differential equation for *B* is

$\sqrt{2m{r}^{2}(E-V)-{L}^{2}}\phantom{\rule{.2em}{0ex}}[\frac{dB}{dr}-\frac{B}{r}]=\pm \sqrt{{R}^{2}-(\frac{LB}{r}{)}^{2}}$

but under the substitution $B=r\mathcal{B}$ this becomes simply

$\pm \frac{d\mathcal{B}}{\sqrt{{R}^{2}-{L}^{2}{\mathcal{B}}^{2}}}=\frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{2m{r}^{2}(E-V)-{L}^{2}}}$

so that the corresponding coefficient function is either of

$B=\{\phantom{\rule{.5em}{0ex}}\begin{array}{l}\frac{rR}{L}sin[\int \frac{L\phantom{\rule{.2em}{0ex}}dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{2m{r}^{2}(E-V)-{L}^{2}}}]\\ \frac{rR}{L}cos[\int \frac{L\phantom{\rule{.2em}{0ex}}dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{2m{r}^{2}(E-V)-{L}^{2}}}]\end{array}$

The differential equation for *A* is a bit more complicated,

$\sqrt{2m{r}^{2}(E-V)-{L}^{2}}\phantom{\rule{.2em}{0ex}}[\frac{dA}{dr}+\frac{m{V}^{\prime}A}{2m(E-V)}]=\pm \frac{m{V}^{\prime}}{2m(E-V)}\sqrt{2m(E-V){R}^{2}-{L}^{2}{A}^{2}}$

but under the substitution $A=\sqrt{2m(E-V)}\mathcal{A}$ the equation becomes

$\pm \frac{d\mathcal{A}}{\sqrt{{R}^{2}-{L}^{2}{\mathcal{A}}^{2}}}=\frac{m{V}^{\prime}dr}{2m(E-V)\sqrt{2m{r}^{2}(E-V)-{L}^{2}}}$

so that the corresponding coefficient function is either of

$A=\{\phantom{\rule{.5em}{0ex}}\begin{array}{l}\frac{R}{L}\sqrt{2m(E-V)}sin[\int \frac{mL{V}^{\prime}dr}{2m(E-V)\sqrt{2m{r}^{2}(E-V)-{L}^{2}}}]\\ \frac{R}{L}\sqrt{2m(E-V)}cos[\int \frac{mL{V}^{\prime}dr}{2m(E-V)\sqrt{2m{r}^{2}(E-V)-{L}^{2}}}]\end{array}$

This last result could be simplified slightly, but is left in this form for a reason that will soon become clear. Note that the quantity under the radical is proportional to squared quantities and so is always positive.

An immediate question at this juncture is why there should be two possible choices for each coefficient function, leading to two invariant vectors. A trivial mathematical answer is that the coupled first-order equations are equivalent to second-order equations, each of which has two linearly independent solutions. A more physically meaningful answer is that in three dimensions angular momentum can be represented as a vector, and there will be two independent directions orthogonal to the angular momentum vector. For the Kepler problem, one can use the Runge vector as a defining constant of the motion, or equally its cross product with the angular momentum vector. The existence of two possible constant vectors continues to hold for higher dimensions, where the orbital motion occurs in an invariant hyperplane spanned by two independent directions.

This leaves the question of how the two possible choices should be correlated for a single invariant vector, and this is best presented at this stage using the constituent imaginary exponentials of the circular functions. First noting that

$\frac{dp}{dr}=\frac{d}{dr}\sqrt{2m(E-V)}=-\frac{m}{p}\frac{dV}{dr}$

the coupled equations can be written symmetrically:

$\sqrt{(\frac{rp}{L}{)}^{2}-1}\frac{dA}{dp}=\frac{p}{L}B\phantom{\rule{4em}{0ex}}\sqrt{(\frac{rp}{L}{)}^{2}-1}\frac{dB}{dr}=\frac{r}{L}A$

Apart from an overall constant factor $\frac{R}{L}$ , the behavior of the coefficient functions can be written symmetrically as well:

$A\sim pexp[\pm i\phantom{\rule{.2em}{0ex}}\int \frac{dp}{p\phantom{\rule{.2em}{0ex}}\sqrt{(\frac{rp}{L}{)}^{2}-1}}]\phantom{\rule{4em}{0ex}}B\sim rexp[\pm i\phantom{\rule{.2em}{0ex}}\int \frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{(\frac{rp}{L}{)}^{2}-1}}]$

If these exponentials are combined with the same sign, the integral is easily evaluated

$\begin{array}{l}\int \frac{dp}{p\phantom{\rule{.2em}{0ex}}\sqrt{(\frac{rp}{L}{)}^{2}-1}}+\int \frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{(\frac{rp}{L}{)}^{2}-1}}=\int \frac{d\left(\frac{rp}{L}\right)}{\frac{rp}{L}\sqrt{(\frac{rp}{L}{)}^{2}-1}}\\ \phantom{\rule{4em}{0ex}}=\frac{i}{2}ln\left[\frac{\sqrt{(\frac{rp}{L}{)}^{2}-1}+i}{\sqrt{(\frac{rp}{L}{)}^{2}-1}-i}\right]\end{array}$

and since the two coefficient functions appear on opposite sides of the differential equations, the exponentials should be anticorrelated in sign for the invariant vector

$\frac{R}{L}[\phantom{\rule{1em}{0ex}}pexp[i\phantom{\rule{.2em}{0ex}}\int \frac{dp}{p\phantom{\rule{.2em}{0ex}}\sqrt{(\frac{rp}{L}{)}^{2}-1}}]\phantom{\rule{.3em}{0ex}}\mathbf{r}-rexp[-i\phantom{\rule{.2em}{0ex}}\int \frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{(\frac{rp}{L}{)}^{2}-1}}]\phantom{\rule{.3em}{0ex}}\mathbf{p}\phantom{\rule{1em}{0ex}}]$

Since this is a classical problem, the corresponding real vectors

$\begin{array}{c}\frac{R}{L}[\phantom{\rule{1em}{0ex}}pcos[\int \frac{dp}{p\phantom{\rule{.2em}{0ex}}\sqrt{(\frac{rp}{L}{)}^{2}-1}}]\phantom{\rule{.3em}{0ex}}\mathbf{r}-rcos[\int \frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{(\frac{rp}{L}{)}^{2}-1}}]\phantom{\rule{.3em}{0ex}}\mathbf{p}\phantom{\rule{1em}{0ex}}]\\ \frac{R}{L}[\phantom{\rule{1em}{0ex}}psin[\int \frac{dp}{p\phantom{\rule{.2em}{0ex}}\sqrt{(\frac{rp}{L}{)}^{2}-1}}]\phantom{\rule{.3em}{0ex}}\mathbf{r}+rsin[\int \frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{(\frac{rp}{L}{)}^{2}-1}}]\phantom{\rule{.3em}{0ex}}\mathbf{p}\phantom{\rule{1em}{0ex}}]\end{array}$

are necessarily constants of the motion. One could at this point rewrite the integrals entirely in terms of the radial variable, but it is more stimulating to leave them in this symmetric form. Since there is only one independent variable of integration in the system, it can be chosen as either the radial variable or the total linear momentum variable by inverting the defining relations:

$p=\sqrt{2m[E-V(r\left)\right]}\phantom{\rule{5em}{0ex}}r={V}^{-1}\phantom{\rule{.2em}{0ex}}[E-\frac{{p}^{2}}{2m}]$

The symmetry between these two variables helps to explain a remarkable way of writing the differential equations for the two coefficient functions *A* and *B*. With the symmetric form of the coupled equations, dividing one equation by the other to eliminate the radical leads to the simple statement

$\frac{d\phantom{\rule{.3em}{0ex}}{A}^{2}}{d\left({p}^{2}\right)}=\frac{d\phantom{\rule{.3em}{0ex}}{B}^{2}}{d\left({r}^{2}\right)}$

If the functional form of *B*^{2} is known, then up to a constant the form of *A*^{2} is identical with *r*^{2} replaced by *p*^{2} or its equivalent function of the radial variable.

It now remains to find a physical meaning for the integrals appearing within the circular functions. The Runge vector is used in the Kepler problem to indicate the direction of either perihelion or aphelion, and its dot product with the radius vector defines the cosine of the angular position of the body. Since dot products are independent of dimension, let the same definition hold here:

$cos\phi \equiv \frac{\mathbf{r}\xb7\mathbf{R}}{rR}=\frac{A{r}^{2}-B(\mathbf{r}\xb7\mathbf{p})}{rR}$

The complementary angle is found with simple trigonometry:

$sin\phi =\sqrt{1-{cos}^{2}\phi}=\frac{BL}{rR}$

Since the Runge vector is identically constant, one can write

$\begin{array}{c}\mathbf{p}\xb7\mathbf{R}=m\stackrel{\xb7}{\mathbf{r}}\xb7\mathbf{R}=m\frac{d}{dt}\mathbf{r}\xb7\mathbf{R}\\ A(\mathbf{r}\xb7\mathbf{p})-B{p}^{2}=m(\stackrel{\xb7}{r}Rcos\phi -rR\stackrel{\xb7}{\phi}sin\phi )\end{array}$

Using previously given quantities, this simplifies to

$\stackrel{\xb7}{\phi}=\frac{L}{m{r}^{2}}$

which expresses the conservation of angular momentum for the motion. Integrating this statement with a derivative above,

$\phi =\int dt\frac{L}{m{r}^{2}}=\int dr\frac{L}{r(\mathbf{r}\xb7\mathbf{p})}=\int \frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{(\frac{rp}{L}{)}^{2}-1}}$

which is precisely the integral appearing in the coefficient function *B*: it is the angular position of the body in motion.

Repeat the procedure with the linear momentum vector replacing the radial vector. Define

$cos\alpha \equiv \frac{\mathbf{p}\xb7\mathbf{R}}{pR}=\frac{A(\mathbf{r}\xb7\mathbf{p})-B{p}^{2}}{pR}$

The complementary angle is again found with simple trigonometry,

$sin\alpha =\sqrt{1-{cos}^{2}\alpha}=\frac{AL}{pR}$

Since the Runge vector is identically constant, one can write

$\begin{array}{c}m\stackrel{\xb7\xb7}{\mathbf{r}}\xb7\mathbf{R}=\stackrel{\xb7}{\mathbf{p}}\xb7\mathbf{R}=\frac{d}{dt}\mathbf{p}\xb7\mathbf{R}\\ -\frac{1}{r}\frac{dV}{dr}[A{r}^{2}-B(\mathbf{r}\xb7\mathbf{p})]=\stackrel{\xb7}{p}Rcos\alpha -pR\stackrel{\xb7}{\alpha}sin\alpha \end{array}$

Using previously given quantities, this time it simplifies to

$\stackrel{\xb7}{\alpha}=\frac{L}{r{p}^{2}}\frac{dV}{dr}$

which is the statement of conservation of angular momentum in terms of the position of the linear momentum vector. Integrating this statement with a derivative above,

$\alpha =\int dt\frac{L}{r{p}^{2}}\frac{dV}{dr}=-\int dp\frac{L}{p(\mathbf{r}\xb7\mathbf{p})}=-\int \frac{dp}{p\phantom{\rule{.2em}{0ex}}\sqrt{(\frac{rp}{L}{)}^{2}-1}}$

which is the negative of the integral appearing in the coefficient function *A*. This integral is also a relative angle, but this time between the Runge vector and the linear momentum vector rather than the radius vector.

Restoring the overall constant factors, the two possible choices for the generalized constant vector can now be written

$\mathbf{R}=\{\phantom{\rule{.5em}{0ex}}\begin{array}{l}\frac{R}{L}\phantom{\rule{.4em}{0ex}}[\phantom{\rule{.4em}{0ex}}pcos\alpha \phantom{\rule{.5em}{0ex}}\mathbf{r}-rcos\phi \phantom{\rule{.5em}{0ex}}\mathbf{p}\phantom{\rule{.3em}{0ex}}]\\ \frac{R}{L}\phantom{\rule{.4em}{0ex}}[-psin\alpha \phantom{\rule{.5em}{0ex}}\mathbf{r}+rsin\phi \phantom{\rule{.5em}{0ex}}\mathbf{p}\phantom{\rule{.3em}{0ex}}]\end{array}$

Using quantities given, it is straightforward to verify that these two vectors are both constants of the motion and are orthogonal. The angles appearing in the vectors have explicit significance with regard to the dynamical variables of the system. And most importantly, this is a general solution for any spherically symmetric potential, not just the Kepler potential. So much for ‘hidden’ symmetries.

As a final note, squaring each vector leads to the symmetric equations

$\begin{array}{r}{cos}^{2}\alpha +{cos}^{2}\phi -2\frac{(\mathbf{r}\xb7\mathbf{p})}{rp}cos\alpha cos\phi -(\frac{L}{rp}{)}^{2}=0\\ {sin}^{2}\alpha +{sin}^{2}\phi -2\frac{(\mathbf{r}\xb7\mathbf{p})}{rp}sin\alpha sin\phi -(\frac{L}{rp}{)}^{2}=0\end{array}$

These are quadratic equations in the circular functions, with symmetric solutions

$\begin{array}{l}cos\alpha =\frac{(\mathbf{r}\xb7\mathbf{p})cos\phi \pm Lsin\phi}{rp}\\ cos\phi =\frac{(\mathbf{r}\xb7\mathbf{p})cos\alpha \pm Lsin\alpha}{rp}\\ sin\alpha =\frac{(\mathbf{r}\xb7\mathbf{p})sin\phi \pm Lcos\phi}{rp}\\ sin\phi =\frac{(\mathbf{r}\xb7\mathbf{p})sin\alpha \pm Lcos\alpha}{rp}\end{array}$

which provide relations between the dynamic angles of the system. These relations simplify the process of correlating circular functions in the invariant vectors. Designating the angle between the radius and linear momentum vectors with γ, one has

$(\mathbf{r}\xb7\mathbf{p})=rpcos\gamma \phantom{\rule{5em}{0ex}}L=\sqrt{{r}^{2}{p}^{2}-(\mathbf{r}\xb7\mathbf{p}{)}^{2}}=rpsin\gamma $

and the relations between the dynamic angles are seen to be simple addition formulae

$\begin{array}{lll}cos\alpha =cos(\phi \pm \gamma )& & cos\phi =cos(\alpha \pm \gamma )\\ sin\alpha =sin(\phi \pm \gamma )& & sin\phi =sin(\alpha \pm \gamma )\end{array}$

where the sign chosen depends on the geometrical arrangement of the three vectors.

Future presentations will consider the relationship of this solution to the angular momentum tensor, visualization of the constancy of the solution relative to the dynamical system variables, and conditions for closed orbits.

*Uploaded 2013.11.22 — Updated 2015.03.31*
analyticphysics.com