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{·}{\phi }=L\to \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{Ldr}{r\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}\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\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 the energy and coupling constant are both negative. Explicitly separating signs, the integral and its evaluation are

$\int \frac{dr}{r\sqrt{\left|\alpha \right|r-\left|\epsilon \right|r^2-1}}={tan}^{-1}\left[\frac{\left|\alpha \right|r-2}{2\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\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 the energy and coupling constant are both positive. With a change of variable $u=r^2$ , the integral to be evaluated is

$\int \frac{dr}{r\sqrt{\epsilon r^2-\alpha r^4-1}}=\frac{1}{2}\int \frac{du}{u\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}\underset{u_{-}}{\overset{u_{+}}{\int }}\frac{du}{u\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 the angular change can be expressed in terms of an inverse Weierstrass elliptic function with integration endpoints determined from an algebraically soluble cubic denominator polynomial. Since the integration endpoints and angular change can always be evaluated numerically to exhibit overall functional behavior, these additional exactly soluble cases do not add enough information for their attendant formulaic overhead.

There are, however, three limiting cases that are very useful. As k → ∞ the two positive roots of the denominator polynomial can be estimated reliably. For a root with value less than unity the second term of the polynomial approaches zero, so that the limiting lower root is determined by

$\epsilon r_{-}^2-1\approx 0\to 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\to 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\sqrt{\epsilon r^2-1}}=\underset{1}{\overset{\sqrt{\epsilon }}{\int }}\frac{du}{u\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 through positive values, the potential disappears and the motion of the orbiting body becomes unconstrained. This can be treated analytically by allowing the upper endpoint to approach infinity. The change in angular position for this limit is

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

Since this is the same value as at k = 2 there must be an extremum in Δφ(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 through negative values

$\Delta \phi (k\to 0^{-})\approx \underset{\frac{1}{\sqrt{\left|\alpha \right|-\left|\epsilon \right|}}}{\overset{\infty }{\int }}\frac{dr}{r\sqrt{\left(\right|\alpha |-|\epsilon \left|\right)r^2-1}}=\underset{1}{\overset{\infty }{\int }}\frac{du}{u\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 the potential is subsumed into the angular momentum term and the motion is again unconstrained. The appropriate endpoints of integration are found by first rationalizing the denominator polynomial,

$\begin{array}{c}\left(\right|\alpha |-1)r^2-\left|\epsilon \right|r^4=0\\ r_{-}=0r_{+}=\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\sqrt{\left|\alpha \right|-1-\left|\epsilon \right|r^2}}\\ \phantom{\Delta \phi (k\to -2)}=\frac{1}{\sqrt{\left|\alpha \right|-1}}\underset{0}{\overset{1}{\int }}\frac{du}{u\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, where the angular change is independent of the other constants of the system. The exponent of the power potential clearly has priority over the other constants, so that depicting the angular change as a function of the exponent will be more informative than as a function of the other constants. The angular change will depend on these other constants between the fixed points as well as the asymptotic value for infinite exponent.

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 is given by a real pole there, and the horizontal asymptote at k = ∞ is a simple additive constant. The remaining behavior can be captured by a term with no real poles and three constants that goes to zero for large 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}\to \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:

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 the discriminant is ${\epsilon }^2-4\alpha$ , while for the Kepler potential k = −1 the discriminant is $|\alpha {|}^2-4|\epsilon |$ , and keeping these discriminants positive corresponds to the choices of reduced constants indicated.

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. Rotating the surface with and without the plane visible also indicates that this separatrix is not exactly the same as the plane, but for purposes of this presentation it is pretty accurate. The same can be see for contours of negative exponent:

This approximate separatrix has the following application: setting r = 1 in the expressions for the contours and requiring that the denominators in the integral for angular change remain real gives the linear inequalities

$\begin{array}{llll}\hfill \epsilon -\alpha -1>0& \to & \alpha <\epsilon -1\hfill & ,k>0\\ \left|\alpha \right|-\left|\epsilon \right|-1>0& \to & \left|\epsilon \right|<\left|\alpha \right|-1& ,k<0\end{array}$

It should be stressed that these linear inequalities on the constant combinations are not exact, but are useful in ensuring that interactive graphics will 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\sqrt{\left|\alpha \right|r^{2-\left|k\right|}-\left|\epsilon \right|r^2-1}}\\ =\frac{2}{2-\left|k\right|}\underset{u_{-}}{\overset{u_{+}}{\int }}\frac{du}{u\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. It is in the nature of cubic splines to exhibit some irregularity in areas where the data points differ markedly, and this may be visible in the interactive graphics.

Numerical evaluation of the angular change for positive exponents produces the following graph:

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.

Returning the derivative of the spline function rather than its functional value, 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 finding the root of the combination numerically.

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 means that two separate calculators will be needed for positive exponents, one specifying a search from the origin to the maximum and the other from the maximum outward.

Uploaded 2014.02.02 — Updated 2017.09.05 analyticphysics.com