The aim of this lesson is to deduce the orbit of a mass say a planet moving under the influence of the inverse square law gravitational attraction of another body say the sun. The solution will imply a deep physical fact namely that a particle in an attractive central force field obeying the inverse square law moves in a path which is a Keplerian conic section. This leads to further questions as to why nature uses inverse square law forces or what the orbit would look like if she instead used an inverse cube law force. Before getting into the solution we first state nine remarkable properties which central forces and the resulting orbits have:

(i) Torque is always zero

Given

$$\mathbf{F}=F(r)\hat{\mathbf{r}}$$

and

$$\mathbf{r}=r\hat{\mathbf{r}}$$

we calculate the torque as the following cross product

$$\mathbf{\tau}=\mathbf{F}\times\mathbf{r}=F(r)\hat{\mathbf{r}}\times r\hat{\mathbf{r}}=F(r)r(\hat{\mathbf{r}}\times\hat{\mathbf{r}})=0$$

remembering the useful pair of relations $\mathbf{a}\cdot\mathbf{a}=a^2$ and $\mathbf{a}\times\mathbf{a}=0$

(ii) Angular Momentum is Conserved

Torque is defined as the rate of change of angular momentum i.e.

$$\mathbf{\tau}=\frac{d\mathbf{L}}{dt}$$

since the torque is zero by (i) we have 

$$\frac{d\mathbf{L}}{dt}=0$$

which has the solution $\mathbf{L}=$ constant. Thus a particle in a central force field moves so that its angular momentum is constant. In other words angular momentum is conserved.

(iii) The motion is confined to a plane perpendicular to the angular momentum vector

Angular momentum is defined as the moment of a particle’s momentum which we calculate as the following cross product

$$\mathbf{L}=\mathbf{r}\times\mathbf{p}=\mathbf{r}\times m\mathbf{v}=m(\mathbf{r}\times\mathbf{v})=m(\mathbf{r}\times\dot{\mathbf{r}})$$

Now since $\mathbf{a}$, $\mathbf{b}$ and $\mathbf{a}\times\mathbf{b}$ are mutually perpendicular and form a right-handed triad we see that $\mathbf{r}$ is perpendicular to $\mathbf{L}$ which is fixed by (ii). This means the radius vector $\mathbf{r}$ and hence the motion of the particle is confined to a plane perpendicular to $\mathbf{L}$

(iv) Kepler’s second law of areas hold

It can be shown that if the force is central then the radius vector sweeps out equal areas in equal time intervals or in other words the areal velocity $\frac{dA}{dt}$ is constant and in fact is equal to $\frac{|\mathbf{L}|}{2m}$

We have four properties so far but can get five more from the fact the central forces that depend only on the distance of the particle from the center are conservative. If you remember your vector calculus you would know that a conservative force satisfies five equivalent conditions (from any one the other four follows). Thus it follows that central forces satisfy these five properties which are stated below.

(v) Central Forces are curl-free

This is straight-forward to show from the definition of the curl of a vector field as follows

$$\mathbf{F}=F(r)\hat{\mathbf{r}}=\frac{F(r)}{r}\mathbf{r}=\frac{F(r)}{r}(x\mathbf{i}+y\mathbf{j}+z\mathbf{k})$$

$$\nabla\times\mathbf{F}=\frac{F(r)}{r}\begin{vmatrix}\mathbf{i}&\mathbf{j}&\mathbf{k}\\ \frac{\partial}{\partial x}&\frac{\partial}{\partial y}&\frac{\partial}{\partial z}\\x&y&z\end{vmatrix}=\frac{F(r)}{r}\left[\left(\frac{\partial z}{\partial y}\mathbin{-}\frac{\partial y}{\partial z}\right)\mathbf{i}+\left(\frac{\partial x}{\partial z}\mathbin{-}\frac{\partial z}{\partial x}\right)\mathbf{j}+\left(\frac{\partial y}{\partial x}\mathbin{-}\frac{\partial x}{\partial y}\right)\mathbf{k}\right]=0$$

Where we have used the fact that $x$, $y$ and $z$ are independent meaning if we vary say $y$ then $x$ does not change so $\frac{\partial x}{\partial y}=0$ etc.

(vi) The work done by $\mathbf{F}$ getting the particle from point A to point B in the field is independent of the path taken.

(vii) The work done by $\mathbf{F}$ when the particle traverses a closed loop is zero

This is easy to show using (vi) 

Let $W_{AB}$ represent the work done in getting from point A to point B by any path, then we can calculate the work done around a closed loop as

$$W=W_{AB}+W_{BA}=W_{AB}\mathbin{-}W_{AB}=0$$

(viii) There exists a scalar function $U$ called the potential defined up to the addition of a constant such that $\mathbf{F}=\mathbin{-}\nabla U$

(ix) The total mechanical energy $E=$ K.E. + P.E. is constant

To calculate the particle’s orbit we need to write the Lagrangian $L=T-U$ and then write and solve the resulting Euler-Lagrange equations of motion. The kinetic energy $T$ is given by

$$T=\frac{1}{2}mv^2=\frac{1}{2}m\dot{r}^2=\frac{1}{2}m(\dot{x}^2+\dot{y}^2+\dot{z}^2)$$

Now our equation for potential uses spherical coordinates which uses the variables $r, \theta, \phi$ whereas our equation for kinetic energy uses cartesian coordinates and involves $x, y, z$ so the next order of business is to write $T$ in spherical coordinates. The transformation equations are

$$x=r\sin\theta\cos\phi$$

$$y=r\sin\theta\sin\phi$$

$$z=r\cos\theta$$

We find the velocities $\dot{x}, \dot{y}, \dot{z}$ by differentiation using the product rule of elementary calculus

$$\dot{x}=\dot{r}\sin\theta\cos\phi+r\cos\theta\cos\phi\dot{\theta}\mathbin{-}r\sin\theta\sin\phi\dot{\phi}$$

$$\dot{y}=\dot{r}\sin\theta\sin\phi+r\cos\theta\sin\phi\dot{\theta}+r\sin\theta\cos\phi\dot{\phi}$$

$$\dot{z}=\dot{r}cos\theta-r\sin\theta\dot{\theta}$$

Now to find the total velocity squared $v^2=\dot{x}^2+\dot{y}^2+\dot{z}^2$ we need to square these three equations and add. When we square the first (or the second) equation we get six terms i.e. three product terms and three terms that are squares of the three terms in the original equation for $\dot{x}$ (or $\dot{y}$). When we square the third equation we get one product term and two terms that are squares of the two terms in the equation for $\dot{z}$. As an exercise we ask you to show that the seven product terms neatly cancel giving zero. The eight remaining terms involving squares can be simplified with repeated use of $\sin^2\theta+\cos^2\theta=1$ to give

$$v^2=\dot{x}^2+\dot{y}^2+\dot{z}^2=\dot{r}^2+r^2\dot{\theta}^2+r^2\sin^2\theta\dot{\phi}^2$$

Finally we can write the Lagrangian for a particle in the field of an attractive inverse square law force in all its glory

$$L=\frac{1}{2}m(\dot{r}^2+r^2\dot{\theta}^2+r^2\sin^2\theta\dot{\phi}^2)+\frac{a}{r}$$

We can use this Lagrangian to write the three Euler-Lagrange equations for the system namely the $r$ equation, the $\theta$ equation and the $\phi$ equation. For now we will just write the $\phi$ equation. Notice that $\phi$ does not appear in the Lagrangian which depends only on $r$, $\dot{r}$, $\theta$, $\dot{\theta}$, and $\dot{\phi}$. Thus $L$ is independent of $\phi$ which is called a cyclic coordinate and it follows that

$$\frac{\partial L}{\partial\phi}=0$$

It’s straight-forward to calculate that

$$\frac{\partial L}{\partial\dot{\phi}}=mr^2\sin^2\theta\dot{\phi}^2$$

The Euler-Lagrange Equation for the $\phi$ coordinate is 

$$\frac{d}{dt}\frac{\partial L}{\partial\dot{\phi}}=\frac{\partial L}{\partial\phi}$$

Thus we get

$$\frac{d}{dt}(mr^2\sin^2\theta\dot{\phi}^2)=0$$

this means the expression in parentheses is a constant which we will call $l$ and therefore write

$$mr^2\sin^2\theta\dot{\phi}=l$$

Now it’s for us to decide what this equation means. Recall it was stated above that the angular momentum vector $\mathbf{L}$ is constant. This means all three cartesian components $L_x$, $L_y$ and $L_z$ are constant. The equation $mr^2\sin^2\theta\dot{\phi}=l$ merely says that the $z$ component of angular momentum is constant because the left-hand side is exactly $L_z$ as shown below.

$$\mathbf{L}=\mathbf{r}\times\mathbf{p}=\mathbf{r}\times m\mathbf{v}=m(\mathbf{r}\times\dot{\mathbf{r}})=m\begin{vmatrix}\mathbf{i}&\mathbf{j}&\mathbf{k}\\x&y&z\\ \dot{x}&\dot{y}&\dot{z}\end{vmatrix}$$

So 

$$L_z=m(x\dot{y}\mathbin{-}y\dot{x})=mr^2sin^2\theta\dot{\phi}$$

We ask the reader to verify the last step using the equations for $x$, $y$, $\dot{x}$ and $\dot{y}$ given above.

Our next step is to orient our $z$ axis so that it is parallel to the constant angular momentum vector $\mathbf{L}$. When this is done the $x$ and $y$ components of $\mathbf{L}$ become zero

$$L_x=L_y=0$$

and the $z$ component is just

$$L_z=|\mathbf{L}|\mathbf{k}$$

so we may write the angular momentum vector $\mathbf{L}$ as

$$\mathbf{L}=L_x\mathbf{i}+L_y\mathbf{j}+L_z\mathbf{k}=0\mathbf{i}+0\mathbf{j}+|\mathbf{L}|\mathbf{k}=|\mathbf{L}|\mathbf{k}$$

we next take the dot product of the angualar momentum vector $\mathbf{L}$ with the position vector $\mathbf{r}$. This dot product must equal zero since the vectors are perpendicular owing to the fact that the plane of the motion is perpendicular to $\mathbf{L}$

$$\mathbf{L}\cdot\mathbf{r}=|\mathbf{L}|\mathbf{k}\cdot(x\mathbf{i}+y\mathbf{j}+z\mathbf{k})=|\mathbf{L}|z=0$$

Now since $|\mathbf{L}|$ is not zero we have

$$z=0$$

or in spherical coordinates

$$r\cos\theta=0$$

i.e. 

$$\theta=\frac{\pi}{2}$$

$$\dot{\theta}=0$$

Thus it turns out that when we choose the $z$ axis to be parallel to the angular momentum vector we have $z=0$ or equivalently $\theta=\frac{\pi}{2}$. This just means the motion is confined to the $xy$ plane. We may use the values of $\theta$ and $\dot{\theta}$ to perform some crucial simplifications.

$$mr^2\sin^2\theta\dot{\phi}=l$$ 

becomes

$$mr^2\dot{\phi}=l$$

The Lagrangian

$$L=\frac{1}{2}m(\dot{r}^2+r^2\dot{\theta}^2+r^2\sin^2\theta\dot{\phi})+\frac{a}{r}$$

simplifies to

$$L=\frac{1}{2}m(\dot{r}^2+r^2\dot{\phi}^2)+\frac{a}{r}$$

We now write the $r$ or radial equation for the simplified Lagrangian. First we calculate

$$\frac{\partial L}{\partial r}=mr\dot{\phi}^2\mathbin{-}\frac{a}{r^2}$$

$$\frac{\partial L}{\partial\dot{r}}=m\dot{r}$$

Since the Euler-Lagrange equation for $r$ is

$$\frac{d}{dt}\frac{\partial L}{\partial\dot{r}}=\frac{\partial L}{\partial r}$$

we get

$$m\ddot{r}=mr\dot{\phi}^2\mathbin{-}\frac{a}{r}$$

or dividing throughout by $m$ and rearranging

$$\ddot{r}\mathbin{-}r\dot{\phi}^2=\mathbin{-}\frac{a}{mr^2}$$

Our task from now on is to solve this second order ordinary differential equation. First we make the following change of variable

$$u=\frac{1}{r}$$

We rewrite the equation $mr^2\dot{\phi}=l$ we found above as

$$\dot{\phi}=\frac{lu^2}{m}$$

We will refer to this as the $\dot{\phi}$ equation for brevity. Over the next three lines we deduce an expression for $\ddot{r}$ in terms of $u$ and $\phi$

$$\frac{dr}{d\phi}=\frac{dr}{du}\frac{du}{d\phi}=\frac{d}{du}\frac{1}{u}\frac{du}{d\phi}=\mathbin{-}\frac{1}{u^2}\frac{du}{d\phi}$$

$$\dot{r}=\frac{dr}{dt}=\frac{dr}{d\phi}\frac{d\phi}{dt}=\mathbin{-}\frac{1}{u^2}\frac{du}{d\phi}\dot{\phi}=\mathbin{-}\frac{1}{u^2}\frac{du}{d\phi}\frac{lu^2}{m}=\mathbin{-}\frac{l}{m}\frac{du}{d\phi}$$

$$\ddot{r}=\mathbin{-}\frac{l}{m}\frac{d}{dt}\frac{du}{d\phi}=\mathbin{-}\frac{l}{m}\frac{d}{d\phi}\frac{du}{d\phi}\frac{d\phi}{dt}=\mathbin{-}\frac{l}{m}\frac{d^2u}{d\phi^2}\dot{\phi}=\mathbin{-}\frac{l}{m}\frac{d^2u}{d\phi^2}\frac{lu^2}{m}=\mathbin{-}u^2\frac{l^2}{m^2}\frac{d^2u}{d\phi^2}$$

where in the penultimate step in lines two and three we made use of the $\dot{\phi}$ equation. Substituting the expression we just found for $\ddot{r}$ as well as $r=\frac{1}{u}$ and $\dot{\phi}=\frac{lu^2}{m}$ into the second order ordinary differential equation we are trying to solve now gives

$$\mathbin{-}u^2\frac{l^2}{m^2}\frac{d^2u}{d\phi^2}\mathbin{-}\frac{1}{u}\left(\frac{lu^2}{m}\right)^2=\mathbin{-}\frac{au^2}{m}$$

multiplying throughout by $\mathbin{-}\frac{m^2}{l^2u^2}$ we finally get a second order ordinary differential equation that is easy to solve:

$$\frac{d^2u}{d\phi^2}+u=\frac{am}{l^2}$$

The solution to this equation is found by patching together the general solution satisfying the homogeneous equation with right-hand side equal to zero and any particular solution that satisfies the original equation. Notice that $u=\frac{am}{l^2}$ is a particular solution because this means that $u$ is a constant so all orders of derivatives with respect to $\phi$ are zero. The general solution to 

$$\frac{d^2u}{d\phi^2}+u=0$$

is found by solving the characteristic equation

$$\lambda^2+1=0$$

which gives 

$$\lambda=\pm i$$

which means that the general solution can be written as

$$u=A\cos(\phi+\epsilon)$$

Tacking the particular solution unto this gives a final solution of

$$u=\frac{am}{l^2}+A\cos(\phi+\epsilon)=\frac{am}{l^2}(1+e(\cos\phi+\epsilon))$$

Where we have renamed $\frac{Al^2}{am}$ as $e$ for convenience. We now have

$$r=\frac{1}{u}=\frac{l^2}{am(1+e\cos(\phi+\epsilon)}=\frac{\rho}{1+e\cos(\phi+\epsilon)}$$

Where we have set $\rho=\frac{l^2}{am}$. The constant $\epsilon$ is not a new physical paramenter. It merely specifies the orientation of the orbit to whichever direction you choose as $\phi=0$. Since the central potential is rotationally symmetric you can redefine $\phi$ as $\phi=\phi+\epsilon$. This makes $\epsilon$ go away giving

$$r=\frac{\rho}{1+e\cos\phi}$$

We now have the polar equation of a conic section with one focus at the origin. The constants $\rho$ and $e$ depend on initial conditions namely the particle’s position and velocity at time $t=0$. $\rho$ is positive, has the dimensions of length and is known as the semi-latus rectum. $0\leq e\leq 1$ is the eccentricity of the orbit. Using conservation of Energy from the Lagrangian one can show that the particle’s energy is given by

$$E=\frac{am^2(e^2\mathbin{-}1)}{2l^2}$$

The following table summarizes the various orbits and energies for different values of $e$. Negative Energy just means the negative P.E. is greater in magnitude than the positive K.E.

\begin{array}{|c|c|c|}\hline\textit{e}&\text{Orbit}&\text{Energy}\\ \hline\text 0&\text{circle}&\text{negative}\\ \hline \text{(0, 1)}&\text{ellipse}&\text{negative}\\ \hline 1&\text{parabola}& 0\\ \hline (1, \infty)&\text{hyperbola}&\text{positive}\\ \hline\end{array}