Consider a body moving in a bound state under the influence of a spherically symmetric power potential. The classical mechanical Hamiltonian for this motion is

$H=\frac{{p}_{r}^{2}}{2m}+\frac{{L}^{2}}{2m{r}^{2}}+a{r}^{k}=E$

with coupling constant *a* and arbitrary real exponent *k*. The motion is independent of the number of spatial dimensions, as long as there are at least two spatial dimensions to allow a definition of angular momentum.

The Kepler potential, of exponent minus one, is particularly distinguished among power potentials. For a body in orbit around a Kepler potential, movement from the point of closest approach to the center of force (perihelion, perigee) to the point of farthest separation from the center of force (aphelion, apogee) corresponds to a change in angular position of exactly π radians. This means that a full circuit of radial values, from closet to farthest back to closest, corresponds to exactly 2π radians, a complete angular cycle, and the orbiting body returns to the same radial and angular positions at the same time. The orbit retraces itself on subsequent periods and is said to be closed.

For exponents differing slightly from minus one, there will be a slight deficit or excess compared to π in moving from closest approach to farthest separation, according to whether the exponent is slightly more positive than minus one or slightly more negative. As the difference from this reference point increases, so does the deficit or excess. At any point where the angular change relative to π is a rational number there will again be closed orbits, but these will have multiple loops in physical space and take multiple circuits to fully retrace.

A similar behavior will occur for positive exponents. The reference point here is the simple harmonic oscillator, of exponent two, for which the movement from closest approach to farthest separation corresponds to a change in angular position of exactly $\frac{\pi}{2}$. The oscillator takes two circuits from closest to farthest back to closest for a complete angular cycle of 2π radians, after which the orbits begin to retrace themselves. The main difference in this behavior from the Kepler potential is that the center of force for the oscillator lies at the center of the elliptical orbits, not at a focus, so that there are two points of closest approach and farthest separation in each complete orbit.

For positive exponents slightly smaller or larger than two, there will again be a deficit or excess compared to $\frac{\pi}{2}$ in moving from closest approach to farthest separation, according to whether the exponent is slightly larger or smaller than two. As the difference from the reference point of two increases, so does the deficit or excess up to a point, and when the angular change relative to π is a rational number there will again be closed orbits taking multiple circuits to fully retrace.

Now to make the description concrete with mathematics. The angular position of the orbiting body is a simple result of angular momentum conservation:

$m{r}^{2}\stackrel{\xb7}{\phi}=L\phantom{\rule{3em}{0ex}}\to \phantom{\rule{3em}{0ex}}\phi =\int dt\frac{L}{m{r}^{2}}$

The integral can be written entirely in terms of the radial variable by evaluating the temporal derivative of the radius from the Hamiltonian,

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

so that the change in the angular position of an orbiting body in moving from closest approach to farthest separation is

$\Delta \phi =\underset{{r}_{-}}{\overset{{r}_{+}}{\int}}\frac{L\phantom{\rule{.3em}{0ex}}dr}{r\phantom{\rule{.3em}{0ex}}\sqrt{2mE{r}^{2}-2ma{r}^{k+2}-{L}^{2}}}$

This same integral appears in the presentation of a generalized Runge vector for spherically symmetric potentials. Introducing the combinations

$\epsilon =\frac{2mE}{{L}^{2}}\phantom{\rule{5em}{0ex}}\alpha =\frac{2ma}{{L}^{2}}$

which will be referred to as reduced constants, the change in angular position is

$\Delta \phi (k;\epsilon ,\alpha )=\underset{{r}_{-}}{\overset{{r}_{+}}{\int}}\frac{dr}{r\phantom{\rule{.3em}{0ex}}\sqrt{\epsilon {r}^{2}-\alpha {r}^{k+2}-1}}$

where the exponent of the power potential has been separated from the other two constants in the notation for a reason that will soon be apparent.

The endpoints of the integration are given by the roots of the polynomial under the radical, which are equivalently the turning points of the motion in the effective potential of the attractive center plus the angular momentum term. For bound states the polynomial will have two real positive roots, but may have additional negative real roots and pairs of imaginary roots. The roots cannot be determined algebraically except for very special cases: this presents a major challenge in describing the analytic behavior of this change as a function of the exponent *k* of the power potential. The bulk of possible real values of the exponent correspond to changes not expressible in terms of known special functions.

For the Kepler potential *k* = −1

$\int \frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{\left|\alpha \right|r-\left|\epsilon \right|{r}^{2}-1}}={tan}^{-1}\left[\frac{\left|\alpha \right|r-2}{2\phantom{\rule{.2em}{0ex}}\sqrt{\left|\alpha \right|r-\left|\epsilon \right|{r}^{2}-1}}\right]$

The endpoints of the integration are

${r}_{\pm}=\frac{\left|\alpha \right|\pm \sqrt{|\alpha {|}^{2}-4|\epsilon |}}{2\left|\epsilon \right|}$

and the total change in angular position in moving from closest approach to farthest separation is

$\begin{array}{l}\Delta \phi =\underset{{r}_{-}}{\overset{{r}_{+}}{\int}}\frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{\left|\alpha \right|r-\left|\epsilon \right|{r}^{2}-1}}\\ \phantom{\Delta \phi}={tan}^{-1}\left[\frac{\left|\alpha \right|{r}_{+}-2}{2\times 0}\right]-{tan}^{-1}\left[\frac{\left|\alpha \right|{r}_{-}-2}{2\times 0}\right]\\ \phantom{\Delta \phi}={tan}^{-1}\left(\infty \right)-{tan}^{-1}(-\infty )\\ \Delta \phi =\pi \end{array}$

which is precisely the exact result described above. This change is independent of the two combinations of constants, depending only on the value of the exponent in the power potential, and thus holds for any body moving in an attractive inverse potential.

For the simple harmonic oscillator *k* = 2

$\int \frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{\epsilon {r}^{2}-\alpha {r}^{4}-1}}=\frac{1}{2}\phantom{\rule{.2em}{0ex}}\int \frac{du}{u\phantom{\rule{.2em}{0ex}}\sqrt{\epsilon u-\alpha {u}^{2}-1}}$

This last integral is exactly the same form as for the Kepler potential, apart from the factor of one half. For the total change in angular position in moving from closest approach to farthest separation, one can thus immediately write

$\Delta \phi =\frac{1}{2}\phantom{\rule{.2em}{0ex}}\underset{{u}_{-}}{\overset{{u}_{+}}{\int}}\frac{du}{u\phantom{\rule{.2em}{0ex}}\sqrt{\epsilon u-\alpha {u}^{2}-1}}=\frac{\pi}{2}$

which is the exact result described above. This result is again independent of constants of the system apart from the exponent, and thus holds for all bodies moving in an attractive squared potential. The only difference in the evaluation from the Kepler potential is the factor of one half resulting from the change of variable of integration.

There are other cases of rational values of the exponent of the power potential that can be evaluated exactly but all involve elliptic functions. For example when *k* = 1

There are, however, three limiting cases that are very useful. As *k* → ∞

$\epsilon {r}_{-}^{2}-1\approx 0\phantom{\rule{2em}{0ex}}\to \phantom{\rule{2em}{0ex}}{r}_{-}\approx \frac{1}{\sqrt{\epsilon}}$

For a root with value greater than unity the second term overwhelms the other two terms, and the upper limiting root is determined by

$\epsilon {r}_{+}^{2}-\alpha {r}_{+}^{k+2}\approx 0\phantom{\rule{2em}{0ex}}\to \phantom{\rule{2em}{0ex}}{r}_{+}\approx \underset{k\to \infty}{lim}(\frac{\epsilon}{\alpha}{)}^{\frac{1}{k}}=(\frac{\epsilon}{\alpha}{)}^{0}=1$

Since the variable of integration takes values less than unity except at the upper endpoint, the asymptotic change of angular position can be estimated by ignoring the second term in the denominator polynomial. The integral is evaluated as

$\begin{array}{l}\Delta \phi (k\to \infty )\approx \underset{\frac{1}{\sqrt{\epsilon}}}{\overset{1}{\int}}\frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{\epsilon {r}^{2}-1}}=\underset{1}{\overset{\sqrt{\epsilon}}{\int}}\frac{du}{u\phantom{\rule{.2em}{0ex}}\sqrt{{u}^{2}-1}}\\ \phantom{\Delta \phi (k\to \infty )}\approx \frac{\pi}{2}-{tan}^{-1}\left[\frac{1}{\sqrt{\epsilon -1}}\right]\end{array}$

This asymptotic angular change ranges from zero to $\frac{\pi}{2}$ as the parameter combination ε varies from one to infinity.

As *k* → 0

$\Delta \phi (k\to {0}^{+})\approx \underset{\frac{1}{\sqrt{\epsilon -\alpha}}}{\overset{\infty}{\int}}\frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{(\epsilon -\alpha ){r}^{2}-1}}=\underset{1}{\overset{\infty}{\int}}\frac{du}{u\phantom{\rule{.2em}{0ex}}\sqrt{{u}^{2}-1}}=\frac{\pi}{2}$

Since this is the same value as at *k* = 2*k*) between these two values of the exponent.

For negative values of the exponent, the energy and coupling constant are also negative in value for bound states. That means that as *k* → 0

$\Delta \phi (k\to {0}^{-})\approx \underset{\frac{1}{\sqrt{\left|\alpha \right|-\left|\epsilon \right|}}}{\overset{\infty}{\int}}\frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{\left(\right|\alpha |-|\epsilon \left|\right){r}^{2}-1}}=\underset{1}{\overset{\infty}{\int}}\frac{du}{u\phantom{\rule{.2em}{0ex}}\sqrt{{u}^{2}-1}}=\frac{\pi}{2}$

for exactly the same result. The function Δφ(*k*) is thus continuous at the origin, but its differentiability there must be ascertained.

As *k* → −2

$\begin{array}{c}\left(\right|\alpha |-1){r}^{2}-\left|\epsilon \right|{r}^{4}=0\\ {r}_{-}=0\phantom{\rule{5em}{0ex}}{r}_{+}=\sqrt{\frac{\left|\alpha \right|-1}{\left|\epsilon \right|}}\end{array}$

so that the evaluation of the angular change is

$\begin{array}{l}\Delta \phi (k\to -2)\approx \underset{0}{\overset{\sqrt{\frac{\left|\alpha \right|-1}{\left|\epsilon \right|}}}{\int}}\frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{\left|\alpha \right|-1-\left|\epsilon \right|{r}^{2}}}\\ \phantom{\Delta \phi (k\to -2)}=\frac{1}{\sqrt{\left|\alpha \right|-1}}\phantom{\rule{.2em}{0ex}}\underset{0}{\overset{1}{\int}}\frac{du}{u\phantom{\rule{.2em}{0ex}}\sqrt{1-{u}^{2}}}\\ \Delta \phi (k\to -2)=\infty \end{array}$

There are now four values of the exponent, *k* = −2, −1, 0 and 2,

Before exhibiting numerical evaluations of the change in angular position as a function of the exponent of the power potential, consider a polynomial estimate to this function. The vertical asymptote at *k* = −2*k* = ∞*k*. In units of π a good candidate is

$f\left(k\right)=\frac{1}{k+2}+\frac{{c}_{1}k+{c}_{2}}{{k}^{2}+{c}_{3}}+c$

where the constants are determined by

$\begin{array}{c}f(-1)=1\\ f\left(0\right)=f\left(2\right)=\frac{1}{2}\end{array}\phantom{\rule{2em}{0ex}}\to \phantom{\rule{2em}{0ex}}\begin{array}{l}{c}_{1}=c\\ {c}_{2}=-4c(6c-1)\\ {c}_{3}=4(6c-1)\end{array}$

In order to avoid a pole in the middle term, *c* must be restricted to range from one sixth to one half. The behavior of this inverse cubic polynomial approximation is given in an interactive graphic using SageMathCell:

This behavior is to be compared with an accurate numerical evaluation. It will turn out to have generally the same appearance, but will differ in the details. This inverse cubic polynomial approximation should be taken in the same spirit as the van der Waals equation of state: qualitatively instructive but quantitatively inaccurate.

Numerical evaluation requires having values for the zeroes of the denominator polynomial as a function of all three arguments. For positive values of the exponent, typical contours for the zeroes are

For negative values of the exponent, typical contours for the zeroes are

Manipulation of the graphics indicates that for a well-posed system with positive exponents, the reduced energy parameter must exceed the reduced coupling constant, but that for negative exponents the relationship between their absolute values is reversed. This is easy enough to understand from the discriminant of the polynomial for two of the exactly soluble cases. For the simple harmonic oscillator *k* = 2*k* = −1

The limitations imposed by discriminant values do not hold for the entire analytic domain of exponents, but one can determine a decent limitation by graphic inspection. Plotting the mathematical expression for the contours of positive exponent in the three-dimensional space of the radial variable and the two reduced constants,

it becomes apparent that there is a separatrix in the system in the vicinity of *r* = 1

This approximate separatrix has the following application: setting *r* = 1

$\begin{array}{llll}\hfill \epsilon -\alpha -1>0& \phantom{\rule{.5em}{0ex}}\to \phantom{\rule{.5em}{0ex}}& \phantom{\rule{.5em}{0ex}}\alpha <\epsilon -1\hfill & \phantom{\rule{.5em}{0ex}},\phantom{\rule{.5em}{0ex}}k>0\\ \left|\alpha \right|-\left|\epsilon \right|-1>0& \phantom{\rule{.5em}{0ex}}\to \phantom{\rule{.5em}{0ex}}& \left|\epsilon \right|<\left|\alpha \right|-1& \phantom{\rule{.5em}{0ex}},\phantom{\rule{.5em}{0ex}}k<0\end{array}$

It should be stressed that these linear inequalities on the constant combinations are not exact, but are useful in ensuring interactive graphics function.

The reversal of the relationship between reduced constants for negative exponents implies that there will be a symmetry relation between pairs of exponents. Explicitly separating negative signs,

$\begin{array}{l}\Delta \phi [-|k|;-|\epsilon |,-|\alpha \left|\right]=\underset{{r}_{-}}{\overset{{r}_{+}}{\int}}\frac{dr}{r\phantom{\rule{.2em}{0ex}}\sqrt{\left|\alpha \right|{r}^{2-\left|k\right|}-\left|\epsilon \right|{r}^{2}-1}}\\ \phantom{\rule{15em}{0ex}}=\frac{2}{2-\left|k\right|}\underset{{u}_{-}}{\overset{{u}_{+}}{\int}}\frac{du}{u\phantom{\rule{.2em}{0ex}}\sqrt{\left|\alpha \right|{u}^{2}-\left|\epsilon \right|{u}^{4/(2-|k\left|\right)}-1}}\\ \Delta \phi [-|k|;-|\epsilon |,-|\alpha \left|\right]=\frac{2}{2-\left|k\right|}\Delta \phi [\frac{2\left|k\right|}{2-\left|k\right|};\left|\alpha \right|,\left|\epsilon \right|]\end{array}$

where the change of variable is $r={u}^{2/(2-|k\left|\right)}$ . The relationship between exponents is curiously the same as for pairs of power potentials in the Schrödinger equation, when allowance is made for the explicit separation of sign. In that presentation it was pointed out that a similar relationship applies to the classical radial action: here is another classical analog.

The method will now be to solve numerically for the endpoints at discrete sampling points, then convert the results of numerical integration at these discrete points into an interpolating spline function for the angular change. Numerical round-off error can produce results that have imaginary parts very close to zero, but since the angular change must be a real value such miniscule imaginary parts are ignored.

And without further ado, numerical evaluation of the angular change for positive exponents produces the following graph, which takes some time to evaluate:

Numerical evaluation of the angular change for negative exponents can be done by direct integration:

Given that there is a relationship between pairs of function values, one for positive and one for negative exponent, this interactive graphic could have been generated with the code for positive exponents along with the relationship. This is what will be used in combining the two graphs to reduce integration time:

The exact behavior is qualitatively the same as the inverse cubic polynomial above, but shows significantly more structure. In particular, the function does not appear to be differentiable at the origin, which will now be exhibited graphically.

Using a custom method of the spline function, the derivative of the angular change for positive exponents is

and the derivative of the angular change for negative exponents is

Manipulation of the two graphics indicates that the two derivatives do not meet at the origin, so that the function is not differentiable there. This rules out attempting to fit the function for both positive and negative exponents with a single differentiable function, such as an inverse fifth-order or higher polynomial.

The lack of differentiablity at the origin can also be shown by differentiating the relationship between positive and negative exponents with respect to the absolute value of the exponent,

$\begin{array}{c}\frac{\partial}{\partial \left|k\right|}[\Delta \phi (-|k\left|\right)=\frac{2}{2-\left|k\right|}\Delta \phi \left(\frac{2\left|k\right|}{2-\left|k\right|}\right)]\\ \Delta {\phi}^{\prime}(-|k\left|\right)=-\frac{2}{(2-|k|{)}^{2}}\Delta \phi \left(\frac{2\left|k\right|}{2-\left|k\right|}\right)-\frac{8}{(2-|k|{)}^{3}}\Delta {\phi}^{\prime}\left(\frac{2\left|k\right|}{2-\left|k\right|}\right)\end{array}$

so that as an exact result one has

$\Delta {\phi}^{\prime}(-0)=-\frac{1}{2}\Delta \phi \left(0\right)-\Delta {\phi}^{\prime}\left(0\right)$

and nondifferentiability at the origin is clear.

One thing remains to complete the presentation: a way of calculating the exponent for given values of reduced energy and coupling constant that will produce a specific angular change. This is done by subtracting the desired constant from the angular change function and using find_root. The function only needs to be recalculated upon modifications to the reduced constants.

Since the angular change function does not have any extrema for negative exponents, one inverse calculator will suffice for this part of the function inversion. The maximum in the function around *k* ≈ 1

Inverse calculator for *k* = −2 to *k* = 0 :

Inverse calculator for *k* = 0 to *k* ≈ 1 :

Inverse calculator for *k* ≈ 1 to *k* = 5 :

Specifying a rational decimal value for the angular change produces the exponent that would result in that angular change. The range of possible angular change is limited to between Δφ ≈ 0.4 and Δφ ≈ 0.6 for small positive exponents, and so is not so interesting. The possibilities for negative exponents are much more interesting, ranging from Δφ = 0.5 up to infinity. In particular, there is a negative exponent for which the orbiting body circles the center of attraction completely in moving from closest approach to farthest separation, which is pretty cool.

And how are these calculators used to determine when closed orbits are possible? A radial circuit was described at the start as the movement from closest approach to farthest separation back to closest approach. For a closed orbit, the total angular change after an integral number of radial circuits must be a multiple of 2π, so that the radial and angular variables return to their initial values at the same time. The condition for a closed orbit is thus

$2\Delta \phi \times {n}_{\mathrm{circuits}}=0\phantom{\rule{.5em}{0ex}}(mod2\pi )$

It is trivial to check by this condition that the simple harmonic oscillator takes two circuits for a closed orbit and the Kepler potential only one. This latter is true of any negative exponent leading to an angular change of an integral unit of π since

$2\left(n\Delta \phi \right)\times 1=0\phantom{\rule{.5em}{0ex}}(mod2\pi )\phantom{\rule{.5em}{0ex}},\phantom{\rule{.5em}{0ex}}n\ne 0\phantom{\rule{.5em}{0ex}}\mathrm{and\; integral}$

The exact value of the exponent will of course depend upon the values of reduced energy and coupling constant, and can be determined by appropriate settings of the calculator for negative exponents. For values of the system constants outside the limits set in the calculators, the complete code embedded in this document is easily modified.

More fun is finding negative exponents that lead to orbits requiring two or more circuits for a closed orbit. If the angular change is half-integral with respect to π, the orbit is closed when

$\Delta \phi =\frac{n}{2}\pi \phantom{\rule{.5em}{0ex}},\phantom{\rule{.5em}{0ex}}n\phantom{\rule{.5em}{0ex}}\mathrm{odd}\phantom{\rule{1em}{0ex}}\Rightarrow \phantom{\rule{1em}{0ex}}{n}_{\mathrm{circuits}}=2$

More generally, when the angular change is an irreducible fraction of π,

$\Delta \phi =\frac{n}{d}\pi \phantom{\rule{.5em}{0ex}},\phantom{\rule{.5em}{0ex}}n,d\phantom{\rule{.5em}{0ex}}\mathrm{coprime}\phantom{\rule{1em}{0ex}}\Rightarrow \phantom{\rule{1em}{0ex}}{n}_{\mathrm{circuits}}=d$

This results satisfies the modular math whether or not the denominator is even. To find conditions for an orbit closed in *d* circuits, set the value of the angular change to any irreducible fraction with this denominator. Again, the exact value of the exponent producing this change will depend upon the values of reduced energy and coupling constant.

*Uploaded 2014.02.02 — Updated 2017.02.05*
analyticphysics.com