Solving the Hamilton-Jacobi Equation

In this section we're going to use the Hamilton-Jacobi equation in order to derive analytical formulae for a motion of a material point in the central Newtonian M/r potential.

The derived formulae can then be used to compare an approximate numerical solution against an analytical, i.e., exact solution to the problem.

The problem is easiest to describe in spherical coordinates, r, $\theta$ and $\phi$. In these coordinates the Hamiltonian assumes the following form:

$$H = \frac{p_r^2}{2} + \frac{p_\theta^2}{2 r^2} + \frac{p_\phi^2}{2 r^2 \sin^2{\theta}} - \frac{M}{r}$$ (4.41)

where we have assumed that the gravitational constant G = 1, and the mass of the material point m = 1 too.

The Hamilton-Jacobi equation that corresponds to that Hamiltonian is

$$\frac{1}{2}\left(\frac{\partial S}{\partial r}\right)^2 + \f... ...}\right)^2 - \frac{M}{r} + \frac{\partial S}{\partial t} = 0$$ (4.42)

The method that is commonly use here is called the separation of variables. The Hamiltonian does not depend explicitly on

• time, t
• angle $\phi$

Therefore a solution is sought in the following form:

$$S = -E t + p_\phi \phi + S_r(r) + S_\theta(\theta) + C,$$ (4.43)

where C is a constant.

Substituting this into the Hamilton-Jacobi equation yields:

$$\frac{1}{2}\left(\frac{\textrm{d} S_r}{\textrm{d} r}\right)^2... ... + \frac{p_\phi^2}{2 r^2\sin^2{\theta}} - \frac{M}{r} - E = 0$$ (4.44)

Multiply this equation by 2 r². Now we can rewrite this equation placing all terms that depend on $\theta$ on the left hand side and all terms that depend on r on the right hand side:

$$\left(\frac{\textrm{d} S_\theta}{\textrm{d} \theta}\right)^2 ... ...r^2 - r^2 \left(\frac{\textrm{d} S_r}{\textrm{d} r}\right)^2$$ (4.45)

Because expressions on both sides of this equation depend on different variables, the equality can hold only if they are equal to the same constant, L²:

$$\left(\frac{\textrm{d} S_\theta}{\textrm{d} \theta}\right)^2 ... ...r^2 - r^2 \left(\frac{\textrm{d} S_r}{\textrm{d} r}\right)^2$$ (4.46)

And this implies that:

$$\frac{\textrm{d} S_r}{\textrm{d} r} \sqrt{2\left(\frac{M}{r} + E - \frac{L^2}{2r^2}\right)}$$ = (4.47)

$$\frac{\textrm{d} S_\theta}{\textrm{d} \theta} \sqrt{L^2 - \frac{p_\phi^2}{\sin^2{\theta}}}$$ = (4.48)

These, in turn, are first order ordinary differential equations, which can be readily integrated.

But recall that there are additional conditions that function S must satisfy. In this case

$$S = S(r, \theta, \phi, P_r, P_\theta, P_\phi, t)$$ (4.49)

where P_r, $P_\theta$ and $P_\phi$ are constants and

$$\frac{\partial S}{\partial P_r}$$ = Q_r (4.50)

$$\frac{\partial S}{\partial P_\theta} Q_\theta$$ = (4.51)

$$\frac{\partial S}{\partial P_\phi} Q_\phi$$ = (4.52)

where Q_r, $Q_\phi$ and $Q_\theta$ are also constants. In our case we have $P_r \sim E$, $P_\phi \sim p_\phi$, and $P_\theta \sim L$, and we can set Q_r, $Q_\phi$ and $Q_\theta$ to zero so that:

   

$$\frac{\partial S}{\partial E}$$ = 0 (4.53)

$$\frac{\partial S}{\partial L}$$ = 0 (4.54)

$$\frac{\partial S}{\partial p_\phi}$$ = 0 (4.55)

The last equation (4.59) yields:

$$0 = \frac{\partial S}{\partial p_\phi} = \phi - \int_{\theta... ...}\theta} {\sin{\theta} \sqrt{\sin^2{\theta} - p^2_\phi/L^2}}$$ (4.56)

It can be proven that this relation is satisfied by a flat motion, i.e., that the material point moves in a plane with vector L perpendicular to that plane. We can therefore change our system of coordinates so that, say, $p_\phi = 0$. Then

$$S = -E t + L\theta + \int_{r_0}^{r_1} \sqrt{2\left(E + \frac{M}{r} - \frac{L^2}{2 r^2}\right)} + C$$ (4.57)

Equation (4.58) yields the shape of the orbit:

$$0 = \theta - \int_{r_0}^{r_1} \frac{L\,\textrm{d}r / r^2} {\sqrt{2\left(E + \frac{M}{r} - \frac{L^2}{2r}\right)}}$$ (4.58)

with the following solution:

$$r = \frac{L^2/M}{1 + e\cos{\theta}}$$ (4.59)

where e is the eccentricity of the orbit:

$$e = \sqrt{1 + \frac{2EL^2}{M^2}}$$ (4.60)

Other parameters pertaining to the orbit are:

$$a = - \frac{M}{2E}$$ (4.61)

Remember that for a trapped particle the energy is negative. a is the semimajor axis of the orbit (when elliptic).

$$b = \frac{L}{\sqrt{-2E}}$$ (4.62)

b is the semiminor axis of the orbit (when elliptic)

$$r_{\textrm{{\scriptsize min}}} = \frac{L^2/M} {1 + \sqrt{1 + \frac{2EL^2}{M^2}}}$$ (4.63)

$r_{\textrm{{\scriptsize min}}}$ is the distance of the closest approach.

Equation (4.57) yields the time dependence of r versus t:

$$0 = -t + \int_{r_0}^{r_1} \frac{\textrm{d}r} {\sqrt{2\left(E + \frac{M}{r} - \frac{L^2}{2 r^2}\right)}}$$ (4.64)

The solution to this equation is quite complicated and can be given in terms of Bessel functions and harmonic motion with the mean circular frequency of

$$\frac{2\pi}{T} = \omega = \sqrt{\frac{M}{a^3}}$$ (4.65)