A discussion by John Baez prompts the question of whether Sundman transformations can give insight into systems other than the Kepler potential. A Sundman transformation is a reparametrization of the temporal variable that was first used as part of an analytic solution to the three-body problem. Such transformations provide an intriguing perspective on the two-body problem for any spherically symmetric 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}$

The dot product of these two vectors

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

is the third rotationally invariant quantity available in a spherically symmetric system. The square of the angular momentum tensor is given by

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

which can 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$

The equation of motion for the radial variable is

$\frac{dr}{dt}=\frac{\partial H}{\partial {p}_{r}}=\frac{{p}_{r}}{m}=\frac{1}{m}\sqrt{2m(E-V)-\frac{{L}^{2}}{{r}^{2}}}$

To determine completely the motion in the invariant plane one needs an angle with respect to a fixed direction. This is given by angular momentum conservation,

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

which is easily integrated for

$\phi =\int dt\frac{L}{m{r}^{2}}=\int dr\frac{L}{r\sqrt{2m{r}^{2}(E-V)-{L}^{2}}}$

With a known functional dependence of the radial variable on the temporal variable this in principle provides a complete solution to the orbital motion for any spherically symmetric potential:

$\mathbf{r}=\left[\phantom{\rule{.3em}{0ex}}r\right(t)cos\phi [r\left(t\right)]\phantom{\rule{.3em}{0ex}},\phantom{\rule{.3em}{0ex}}r(t)sin\phi [r\left(t\right)\left]\phantom{\rule{.3em}{0ex}}\right]$

The problem of course is that the equation of motion for the radial variable, the energy equation

$t=\int dt\frac{mr}{\sqrt{2m{r}^{2}(E-V)-{L}^{2}}}$

is rarely integrable and invertible. The reparametrization of the temporal variable via a Sundman transformation allows one to choose a simple form for the dependence of the radial variable on the new parameter, obviating the issue of integration and inversion.

The solution presented here will be determined in parametric form up to the two quadratures directly above. There will remain the issue of inverting the dependence of the temporal variable on the new parameter, but this same issue occurs in the Kepler potential in the form of the Kepler equation linking time with the eccentric anomaly.

The Sundman transformation will be written

${t}^{\prime}=\frac{dt}{ds}=f\left(r\right)$

where following Baez a prime will be used to denote differentiation with respect to the transformed temporal variable *s*. A dot will be used to denote differentiation with respect to the original temporal variable *t* and derivatives with respect to the radial variable will be written explicitly.

The first derivative of the radial vector with respect to the transformed temporal variable is

${\mathbf{r}}^{\prime}=\frac{d}{ds}\mathbf{r}=\frac{dt}{ds}\stackrel{\xb7}{\mathbf{r}}=f\stackrel{\xb7}{\mathbf{r}}$

The first derivative of the radial variable with respect to the transformed temporal variable is

${r}^{\prime}=\frac{dr}{ds}=\sum _{i=1}^{n}\frac{{r}_{i}}{r}\frac{d{r}_{i}}{ds}=\frac{(\mathbf{r}\xb7{\mathbf{r}}^{\prime})}{r}$

which in terms of the original dynamic variables can be variously written

${r}^{\prime}=\frac{f}{r}(\mathbf{r}\xb7\stackrel{\xb7}{\mathbf{r}})=\frac{f}{mr}(\mathbf{r}\xb7\mathbf{p})=\frac{f}{m}\sqrt{2m(E-V)-\frac{{L}^{2}}{{r}^{2}}}$

The second derivative of the radial vector with respect to the transformed temporal variable is

${\mathbf{r}}^{\prime \prime}=\frac{d}{ds}{\mathbf{r}}^{\prime}=\frac{d}{ds}\left(f\stackrel{\xb7}{\mathbf{r}}\right)={f}^{\prime}\stackrel{\xb7}{\mathbf{r}}+{f}^{2}\stackrel{\xb7\xb7}{\mathbf{r}}$

The second derivative of the radial vector with respect to the original temporal variable is given by the equations of motion

$\stackrel{\xb7\xb7}{\mathbf{r}}=\frac{\stackrel{\xb7}{\mathbf{p}}}{m}=-\frac{1}{mr}\frac{dV}{dr}\mathbf{r}=-g\left(r\right)\mathbf{r}$

so that the second derivative of the radial vector with respect to the transformed temporal variable becomes

${\mathbf{r}}^{\prime \prime}=\frac{{f}^{\prime}}{f}{\mathbf{r}}^{\prime}-g{f}^{2}\mathbf{r}$

The second derivative of the radial variable with respect to the transformed temporal variable can be written

$\begin{array}{l}{r}^{\prime \prime}=\frac{({\mathbf{r}}^{\prime}\xb7{\mathbf{r}}^{\prime})}{r}+\frac{(\mathbf{r}\xb7{\mathbf{r}}^{\prime \prime})}{r}-\frac{{r}^{\prime}}{{r}^{2}}(\mathbf{r}\xb7{\mathbf{r}}^{\prime})\\ \phantom{{r}^{\prime \prime}}=\frac{({\mathbf{r}}^{\prime}\xb7{\mathbf{r}}^{\prime})}{r}-g{f}^{2}r+{r}^{\prime}\phantom{\rule{.3em}{0ex}}[\frac{{f}^{\prime}}{f}-\frac{{r}^{\prime}}{r}]\end{array}$

but will be more useful in the form

${r}^{\prime \prime}=\frac{d}{ds}{r}^{\prime}={r}^{\prime}\frac{d}{dr}{r}^{\prime}=\frac{1}{2}\frac{d}{dr}({r}^{\prime}{)}^{2}$

These two forms are easily shown to be identical using the original dynamic variables.

The third derivative of the radial vector with respect to the transformed temporal variable is

$\begin{array}{l}{\mathbf{r}}^{\prime \prime \prime}=(\frac{{f}^{\prime}}{f}{)}^{\prime}{\mathbf{r}}^{\prime}+\frac{{f}^{\prime}}{f}{\mathbf{r}}^{\prime \prime}-(g{f}^{2}{)}^{\prime}\mathbf{r}-g{f}^{2}{\mathbf{r}}^{\prime}\\ \phantom{{\mathbf{r}}^{\prime \prime \prime}}=[\frac{{f}^{\prime \prime}}{f}-g{f}^{2}]{\mathbf{r}}^{\prime}-[\frac{{f}^{\prime}}{f}g{f}^{2}+(g{f}^{2}{)}^{\prime}]\mathbf{r}\\ {\mathbf{r}}^{\prime \prime \prime}=[\frac{{f}^{\prime \prime}}{f}-g{f}^{2}]{\mathbf{r}}^{\prime}-\frac{(g{f}^{3}{)}^{\prime}}{f}\mathbf{r}\end{array}$

where the second derivative of the reparametrization function is

${f}^{\prime \prime}=\frac{d}{ds}\left[{r}^{\prime}\frac{df}{dr}\right]={r}^{\prime \prime}\frac{df}{dr}+({r}^{\prime}{)}^{2}\frac{{d}^{2}f}{d{r}^{2}}$

One way to determine the reparametrization function $f\left(r\right)$ is to remove the last term in the third derivative of the radial vector with respect to the transformed temporal variable by choosing

$g{f}^{3}=\mathrm{constant}={c}_{1}$

so that the third derivative of the radial vector becomes

${\mathbf{r}}^{\prime \prime \prime}=\frac{1}{f}[{r}^{\prime \prime}\frac{df}{dr}+({r}^{\prime}{)}^{2}\frac{{d}^{2}f}{d{r}^{2}}-{c}_{1}]{\mathbf{r}}^{\prime}$

Since the reparametrization function is here a function of the first derivative of the potential, the quantity in brackets is in general complicated, but there are two very simple particular solutions available. If the the reparametrization function is constant, so that both derivatives in the brackets disappear, then

${\mathbf{r}}^{\prime \prime \prime}=-{\mathbf{r}}^{\prime}\phantom{\rule{5em}{0ex}}f={c}_{1}$

which corresponds to the potential

$\begin{array}{c}\frac{1}{mr}\frac{dV}{dr}=g\left(r\right)=\frac{1}{{c}_{1}^{2}}\\ V\left(r\right)=\frac{m}{2{c}_{1}^{2}}{r}^{2}\equiv \frac{1}{2}m{\omega}^{2}{r}^{2}\end{array}$

This is the potential for the simple harmonic oscillator, and the Sundman transformation in this case is merely a rescaling of the temporal variable. If the reparametrization function is linear in the radial variable, so that the second derivative in the brackets disappears, one has

${\mathbf{r}}^{\prime \prime \prime}=\frac{1}{r}[{r}^{\prime \prime}-{c}_{1}]{\mathbf{r}}^{\prime}\phantom{\rule{5em}{0ex}}f=r$

The potential in this case is

$\begin{array}{c}\frac{1}{mr}\frac{dV}{dr}=g\left(r\right)=\frac{{c}_{1}}{{r}^{3}}\\ V\left(r\right)=-\frac{m{c}_{1}}{r}\equiv -\frac{k}{r}\end{array}$

which is the Kepler potential. The second derivative of the radial variable with respect to the transformed temporal variable is now

$\begin{array}{l}{r}^{\prime \prime}=\frac{1}{2}\frac{d}{dr}({r}^{\prime}{)}^{2}\\ \phantom{{r}^{\prime \prime}}=\frac{d}{dr}{r}^{2}\phantom{\rule{.2em}{0ex}}[\frac{1}{m}(E+\frac{k}{r})-\frac{{L}^{2}}{2{m}^{2}{r}^{2}}]\\ {r}^{\prime \prime}=\frac{1}{m}(2Er+k)\end{array}$

and the third derivative of the radial vector with respect to the transformed temporal variable is

${\mathbf{r}}^{\prime \prime \prime}=\frac{2E}{m}{\mathbf{r}}^{\prime}\phantom{\rule{5em}{0ex}}f=r$

For both the simple harmonic oscillator and the Kepler potential, the Sundman transformation here will produce harmonic motion for the radial vector. This will occur as well for the potential determined by the general solution of the nonlinear differential equation

$\frac{1}{f}[{r}^{\prime \prime}\frac{df}{dr}+({r}^{\prime}{)}^{2}\frac{{d}^{2}f}{d{r}^{2}}-{c}_{1}]={c}_{2}$

which is third order with respect to derivatives of the potential. It will not, however, be true for an arbitrary spherically symmetric potential, so that this Sundman transformation is inadequate for a general solution.

To find a Sundman transformation applicable to all spherically symmetric potentials, note that for the Kepler potential the previous transformation led to the behavior

${r}^{\prime \prime}=-\frac{2(-E)}{m}r+\frac{k}{m}$

for the radial variable. Since energy is negative for bound states of the Kepler potential, this represents harmonic motion with a constant offset. Taking a cue from this, consider a Sundman transformation that produces harmonic motion in the radial variable for all potentials. This can be done by setting the second derivative of the radial variable with respect to the transformed temporal variable in its more useful form equal to a linear with constant offset:

${r}^{\prime \prime}=\frac{1}{2}\frac{d}{dr}({r}^{\prime}{)}^{2}=-r+{c}_{1}$

With the expression above for the first derivative of the radial variable with respect to the transformed temporal variable, the integration is simple:

$\begin{array}{c}\frac{1}{2{m}^{2}}\frac{d}{dr}\left[{f}^{2}[2m(E-V)-\frac{{L}^{2}}{{r}^{2}}]\phantom{\rule{.3em}{0ex}}\right]=-r+{c}_{1}\\ f\left(r\right)=mr\sqrt{\frac{{c}_{0}+2{c}_{1}r-{r}^{2}}{2m{r}^{2}(E-V)-{L}^{2}}}\end{array}$

Since the radial variable is harmonic by construction, one can immediately write

$r={c}_{1}+{c}_{2}coss\equiv a(1-ecoss)$

where parameters analogous to the semimajor axis and eccentricity of the Kepler problem have been introduced. The turning points of the motion are

${r}_{\mathrm{inner}}=a(1-e)\phantom{\rule{5em}{0ex}}{r}_{\mathrm{outer}}=a(1+e)$

and correspond to zeroes of the momentum ${p}_{r}$ conjugate to the radial variable. Evaluating the energy at both end points

$E=\frac{{L}^{2}}{2m{a}^{2}(1-e{)}^{2}}+V\left[a\right(1-e\left)\right]=\frac{{L}^{2}}{2m{a}^{2}(1+e{)}^{2}}+V\left[a\right(1+e\left)\right]$

one can easily determine both the angular momentum

${L}^{2}=\frac{m{a}^{2}(1-{e}^{2}{)}^{2}}{2e}[\phantom{\rule{.3em}{0ex}}V\left[a\right(1+e\left)\right]-V\left[a\right(1-e\left)\right]\phantom{\rule{.3em}{0ex}}]$

and the energy

$E=\frac{1}{4e}[\phantom{\rule{.3em}{0ex}}(1+e{)}^{2}V\left[a\right(1+e\left)\right]-(1-e{)}^{2}V\left[a\right(1-e\left)\right]\phantom{\rule{.3em}{0ex}}]$

in terms of the orbit parameters introduced. The transformation as it now stands determines the original temporal variable as

$t=\int ds\phantom{\rule{.3em}{0ex}}mr\sqrt{\frac{{c}_{0}+2ar-{r}^{2}}{2m{r}^{2}(E-V)-{L}^{2}}}$

with dynamic constants of the motion as given directly above. While this holds for any constant of integration ${c}_{0}$ , this constant can be chosen to simplify the integral. Evaluating part of the numerator under the radical,

$2ar-{r}^{2}={a}^{2}(1-{e}^{2}{cos}^{2}s)={a}^{2}(1-{e}^{2})+(aesins{)}^{2}$

the choice ${c}_{0}={a}^{2}({e}^{2}-1)$ cancels the constant term, with the remainder being the square of the derivative of the radial variable with respect to the transformed temporal variable. The integral for the temporal variable becomes

$t=\int dr\frac{mr}{\sqrt{2m{r}^{2}(E-V)-{L}^{2}}}$

which surprisingly is the original energy equation, but now parametrically dependent on the transformed temporal variable. While at first sight it appears that nothing has really been accomplished, first write out the second derivative of the radial variable with respect to the original temporal variable:

$\stackrel{\xb7\xb7}{r}=\frac{d}{dt}\stackrel{\xb7}{r}=\stackrel{\xb7}{r}\frac{d}{dr}\stackrel{\xb7}{r}=\frac{1}{2}{\stackrel{\xb7}{r}}^{2}=\frac{1}{2{m}^{2}}\frac{d}{dr}[2m(E-V)-\frac{{L}^{2}}{{r}^{2}}]$

The final right-hand side is a fully determined function of the radial variable, with no leeway to assign a convenient parametric form for the radial variable. The introduction of the Sundman transformation provides this leeway in a consistent manner, allowing any dependence of the radial variable to be chosen. The simple harmonic motion has a significant advantage in that it characterizes the turning points efficiently for evaluation of the energy and angular momentum in terms of parameters introduced.

The complete solution includes the integral for the angle in the invariant plane

$\phi =\int dr\frac{L}{r\sqrt{2m{r}^{2}(E-V)-{L}^{2}}}$

as a function parametrically dependent on the transformed temporal variable as well. The integrals can be evaluated with respect to either the radial variable or the transformed temporal variable to determine the motion parametrically. This can be done explicitly for only a handful of potentials, but represents a complete solution in integral form.

The parametric determination of the motion is an extension of the classical eccentric anomaly, which is normally described by circumscribing the elliptical orbit of the Kepler potential. It should more accurately be described with a circular orbit halfway between the inner and outer limits of the actual, and generally precessing, orbit. The normal description holds because the eccentric anomaly sweeps on the same angle on all circles centered with the elliptical orbit. The extension of the eccentric anomaly is labeled “Keplerization” in a conference paper by Vladimir Martinusi and Pini Gurfil.

In illustration of the method, orbits in terms of the transformed temporal variable will be evaluated for two simple cases, the Kepler potential and the simple harmonic oscillator. Using expressions above, the energy and angular momentum for the Kepler potential are

$E=-\frac{k}{2a}\phantom{\rule{5em}{0ex}}{L}^{2}=mak(1-{e}^{2})$

The radical for denominators of integrals is

$\sqrt{2m{r}^{2}(E-V)-{L}^{2}}=\sqrt{mak}\phantom{\rule{.3em}{0ex}}esins$

The integral for the temporal variable is

$\begin{array}{l}t=\sqrt{\frac{m{a}^{3}}{k}}\int ds(1-ecoss)=\sqrt{\frac{m{a}^{3}}{k}}(s-esins)\end{array}$

which is the Kepler equation as expected. The integral for the angle in the invariant plane is

$\phi =\sqrt{1-{e}^{2}}\int \frac{ds}{1-ecoss}=2{tan}^{-1}\left[\sqrt{\frac{1+e}{1-e}}\frac{sins}{1+coss}\right]$

which can be verified by differentiation. Circular functions in the invariant plane are

$cos\phi =\frac{coss-e}{1-ecoss}\phantom{\rule{5em}{0ex}}sin\phi =\frac{\sqrt{1-{e}^{2}}sins}{1-ecoss}$

Multiplying by the parametrization of the radial variable, the radial vector for the Kepler potential is

${\mathbf{r}}_{\mathrm{Kepler}}=[\phantom{\rule{.3em}{0ex}}a(coss-e)\phantom{\rule{.3em}{0ex}},\phantom{\rule{.3em}{0ex}}a\sqrt{1-{e}^{2}}sins\phantom{\rule{.3em}{0ex}}]$

For the two independent coordinates in the invariant plane, the orbit is

$\frac{(x+ae{)}^{2}}{{a}^{2}}+\frac{{y}^{2}}{{a}^{2}(1-{e}^{2})}=1$

which is an ellipse displaced along the major axis as expected. When thinking about this system in terms of a fourth dimension, the eccentric circle is rotated into the fourth dimension at an angle determined by the eccentricity.

The energy and angular momentum for the simple harmonic oscillator are

$E=m{\omega}^{2}{a}^{2}(1+{e}^{2})\phantom{\rule{4em}{0ex}}{L}^{2}=\left[m\omega {a}^{2}\right(1-{e}^{2}){]}^{2}$

Integrals for this systems are more easily evaluated in terms of the radial variable to avoid elliptic functions. The radical for denominators of integrals needs to be written in two different ways:

$\begin{array}{l}\sqrt{2m{r}^{2}(E-V)-{L}^{2}}=m\omega \sqrt{4{e}^{2}{a}^{4}-[{r}^{2}-{a}^{2}(1+{e}^{2}){]}^{2}}\\ \phantom{\rule{4em}{0ex}}=m\omega {r}^{2}\sqrt{\frac{4{e}^{2}}{(1-{e}^{2}{)}^{2}}-[\frac{{a}^{2}(1-{e}^{2})}{{r}^{2}}-\frac{1+{e}^{2}}{1-{e}^{2}}{]}^{2}}\end{array}$

Using the first form of the radical, the integral for the temporal variable is

$\begin{array}{l}t=\frac{1}{2\omega}\int \frac{d\left({r}^{2}\right)}{\sqrt{4{e}^{2}{a}^{4}-[{r}^{2}-{a}^{2}(1+{e}^{2}){]}^{2}}}\\ \phantom{t}=\frac{1}{2\omega}{sin}^{-1}\left[\frac{{r}^{2}-{a}^{2}(1+{e}^{2})}{2e{a}^{2}}\right]+{c}_{0}\\ t=\frac{1}{2\omega}{sin}^{-1}\left[\frac{(1-ecoss{)}^{2}-(1+{e}^{2})}{2e}\right]+{c}_{0}\end{array}$

where a constant of integration has been included for later use. Using the second form of the radical, the integral for the angle in the invariant plane is

$\begin{array}{l}\phi =\int dr\frac{{a}^{2}(1-{e}^{2})}{{r}^{3}\sqrt{\frac{4{e}^{2}}{(1-{e}^{2}{)}^{2}}-[\frac{{a}^{2}(1-{e}^{2})}{{r}^{2}}-\frac{1+{e}^{2}}{1-{e}^{2}}{]}^{2}}}\\ \phantom{\phi}=-\frac{1}{2}\int d\left[\frac{{a}^{2}(1-{e}^{2})}{{r}^{2}}\right]\frac{1}{\sqrt{\frac{4{e}^{2}}{(1-{e}^{2}{)}^{2}}-[\frac{{a}^{2}(1-{e}^{2})}{{r}^{2}}-\frac{1+{e}^{2}}{1-{e}^{2}}{]}^{2}}}\\ \phantom{\phi}=\frac{1}{2}{cos}^{-1}\left[\frac{{a}^{2}(1-{e}^{2}{)}^{2}-(1+{e}^{2}){r}^{2}}{2e{r}^{2}}\right]\\ \phi =\frac{1}{2}{cos}^{-1}\left[\frac{(1-{e}^{2}{)}^{2}-(1+{e}^{2}\left)\right(1-ecoss{)}^{2}}{2e(1-ecoss{)}^{2}}\right]\end{array}$

Circular functions in the invariant plane are

$\begin{array}{l}cos\phi =\frac{(1-e)}{2\sqrt{e}(1-ecoss)}\sqrt{(1+e{)}^{2}-(1-ecoss{)}^{2}}\\ sin\phi =\frac{(1+e)}{2\sqrt{e}(1-ecoss)}\sqrt{(1-ecoss{)}^{2}-(1-e{)}^{2}}\end{array}$

Multiplying by the parametrization of the radial variable, the radial vector for the simple harmonic oscillator is

$\begin{array}{l}{\mathbf{r}}_{\mathrm{SHO}}=[\phantom{\rule{.3em}{0ex}}\frac{a(1-e)}{2\sqrt{e}}\sqrt{(1+e{)}^{2}-(1-ecoss{)}^{2}}\phantom{\rule{.3em}{0ex}},\\ \phantom{\rule{6em}{0ex}}\frac{a(1+e)}{2\sqrt{e}}\sqrt{(1-ecoss{)}^{2}-(1-e{)}^{2}}\phantom{\rule{.3em}{0ex}}]\end{array}$

For the two independent coordinates in the invariant plane, the orbit is

$\frac{{x}^{2}}{{a}^{2}(1-e{)}^{2}}+\frac{{y}^{2}}{{a}^{2}(1+e{)}^{2}}=1$

which is an ellipse centered at the origin. As the parameter *s* progresses from zero to π only one quadrant of the ellipse is traversed, so that the entire ellipse requires the parameter to vary from zero to 4π. This differs from a period of 2π for the Kepler orbit and explains the difference between them in the location of the center of force.

Unlike the Kepler equation, the transformed temporal variable for the simple harmonic oscillator is easily inverted. Adjusting the constant of integration so that zeroes of both variables coincide, one has

$(1-ecoss{)}^{2}=1+{e}^{2}-2ecos2\omega t$

so that the radial vector becomes

$\begin{array}{l}{\mathbf{r}}_{\mathrm{SHO}}=[\phantom{\rule{.3em}{0ex}}a(1-e)\sqrt{\frac{1+cos2\omega t}{2}}\phantom{\rule{.3em}{0ex}},\phantom{\rule{.3em}{0ex}}a(1+e)\sqrt{\frac{1-cos2\omega t}{2}}\phantom{\rule{.3em}{0ex}}]\\ \phantom{{\mathbf{r}}_{\mathrm{SHO}}}=[\phantom{\rule{.3em}{0ex}}a(1-e)cos\omega t\phantom{\rule{.3em}{0ex}},\phantom{\rule{.3em}{0ex}}a(1+e)sin\omega t\phantom{\rule{.3em}{0ex}}]\end{array}$

and the moderately complicated intermediate parametrization leads directly back to the expected simple circular functions.

As an aside, it is easy to show from the expression for energy that the Kepler potential is the only power potential for which energy is independent of eccentricity.

To recap, the Sundman transformation

$\frac{dt}{ds}=mr\sqrt{\frac{{a}^{2}({e}^{2}-1)+2ar-{r}^{2}}{2m{r}^{2}(E-V)-{L}^{2}}}$

along with expressions for the energy

$E=\frac{1}{4e}[\phantom{\rule{.3em}{0ex}}(1+e{)}^{2}V\left[a\right(1+e\left)\right]-(1-e{)}^{2}V\left[a\right(1-e\left)\right]\phantom{\rule{.3em}{0ex}}]$

and the angular momentum

${L}^{2}=\frac{m{a}^{2}(1-{e}^{2}{)}^{2}}{2e}[\phantom{\rule{.3em}{0ex}}V\left[a\right(1+e\left)\right]-V\left[a\right(1-e\left)\right]\phantom{\rule{.3em}{0ex}}]$

determines harmonic motion in the radial variable

$r=a(1-ecoss)$

for any spherically symmetric potential $V\left(r\right)$ in conjunction with the two quadratures

$t=\int dr\frac{mr}{\sqrt{2m{r}^{2}(E-V)-{L}^{2}}}$

$\phi =\int dr\frac{L}{r\sqrt{2m{r}^{2}(E-V)-{L}^{2}}}$

The motion is determined parametrically in an extension of the classical eccentric anomaly.

*Uploaded 2015.08.15 — Updated 2015.08.17*
analyticphysics.com