The motion of the sun and the earth is an example of a “two-body problem”. This is a relatively simple problem that can be solved analytically (it interesting to notice here, however, that adding a new object to the problem, the moon for instance, make it completely intractable. This is the famous “three body problem”). We can assume that, to a good approximation, the sun is stationary and is a convenient origin of our coordinate system. This is equivalent to changing to a “center of mass” coordinate system, where most of the mass is concentrated in the sun. The problem can be reduced to an equivalent one body problem involving an object of reduced mass $\mu$ given by $$\mu=\frac{mM}{m+M}$$ Since the mass of the earth is $m=5.99\times 10^{24}$ kg and the mass of the sun is $M=1.99\times 10^{30}$ kg we find that for most practical purposes, the reduced mass of the earth-sun system is that of the earth. Hence, in the following we are going to consider the problem of a single particle of mass $m$ moving about a fixed center of force, which we take as the origin of the coordinate system. The gravitational force on the particle m is given by $${\mathbf F}=-\frac{GMm}{r^3}{\mathbf r},$$ where the vector ${\mathbf r}$ is directed from $M$ to $m$, and $G$ is the gravitation constant $$G=6.67\times 10^{-11} \frac{m^3}{kg.s^2}$$ The negative sign implies that the gravitational force is attractive, and decreases with the separation $r$. The gravitational force is a “central force”: its magnitude depends on the separation between the particles and its direction is along the line that connects them. The assumption is that the motion is confined to the $xy$ plane. The angular momentum ${\mathbf L}$ lies on the third direction $z$ and is a constant of motion, *i.e.* it is conserved: $$L_z=({\mathbf r}\times m{\mathbf v})_z=m(xv_y-yv_x)=\mathrm{const.}$$ An additional constant of motion ins the total energy $E$ given by $$E=\frac{1}{2}mv^2-\frac{GmM}{r}$$ If we fix the coordinate system in the sun, the equation of motion is $$m\frac{d^2{\mathbf r}}{dt^2}=-\frac{mMG}{r^3}{\mathbf r}$$ For computational purposes it is convenient to write it down in cartesian components: $$\begin{aligned} && F_x=-\frac{GMm}{r^2}\cos{\theta}=-\frac{GMm}{r^3}x, \\ && F_y=-\frac{GMm}{r^2}\sin{\theta}=-\frac{GMm}{r^3}y.\end{aligned}$$ Hence, the equations of motions in cartesian coordinates are: $$\begin{aligned} && \frac{d^2x}{dt}=-\frac{GM}{r^3}x, \\ && \frac{d^2y}{dt}=-\frac{GM}{r^3}y, \end{aligned}$$ where $r^2=x^2+y^2$. These are coupled differential equations, since each differential equation contains both $x$ and $y$.
Since many planetary orbits are nearly circular, it is useful to obtain the condition for a circular orbit. In this case, the magnitude of the acceleration ${\mathbf a}$ is related to the radius by $$a=\frac{v^2}{r}$$ where $v$ is the speed of the object. The acceleration is always directed toward the center. Hence $$\frac{mv^2}{r}=\frac{mMG}{r^2}$$ or $$v=(\frac{MG}{r})^{1/2}. $$ This is a general condition for the circular orbit. We can also find the dependence of the period $T$ on the radius of a circular orbit. Using the relation $$T=\frac{2\pi r}{v},$$ we obtain $$T^2=\frac{4\pi^2r^3}{GM}$$
An ellipse has two foci $F_1$ and $F_2$, and has the property that for any point the distance $F_1P+F_2P$ is a constant. It also has a horizontal semi-axis $a$ and a vertical $b$. It is common in astronomy to characterize an orbit by its “eccentricity” $e$, given by the ratio of the distance between the foci, and the length of the major axis $2a$. Since $F_1P+F_2P=2a$, it is easy to show that (consider a point $P$ at $x=0$,$y=b$) $$e=\sqrt{1-\frac{b^2}{a^2}},$$ with $0<e<1$. A special case is $a=b$ for which the ellipse reduces to a circle and $e=0$. The earth orbit has eccentricity $e=0.0167$.
It is useful to choose a system of units where the product $GM$ is of the order of unity. To describe the earth’s motion, the convention is to choose the earth’s semi-major axis as the unit of length, called “astronomical unit” (AU) and is $$1AU=1.496 \times 10^{11}m.$$ The unit of time is taken to be “one year”, or $3.15 \times 10^7$s. In these units, $T=1$yr, $a=1AU$, and we can write $$GM=\frac{4\pi ^2a^3}{T^2}=4\pi ^2 AU^3/yr^2.$$
Write a program to simulate motion in a central force field. Verify the case of circular orbit using (in astronomical units) ($x_0=1$, $y_0=0$ and $v_x(t=0)=0$. Use the condition ([circular] to calculate $v_y(t=0)$ for a circular orbit. Choose a value of $\Delta t$ such that to a good approximation the total energy $E$ is conserved. Is your value of $\Delta t$ small enough to reproduce the orbit over several periods?
Run the program for different sets of initial conditions $x_0$ and $v_y(t=0)$ consistent with the condition for a circular orbit. Set $y_0=0$ and $v_x(t=0)=0$. For each orbit, measure the radius and the period to verify Kepler’s third law ($T^2/a^3=\mathrm{const.}$). Think of a simple condition which allows you to find the numerical value of the period.
Show that Euler’s method does not shield stable orbits for the same choice of $\Delta t$ used in the previous items. Is it sufficient to simply choose a smaller $\Delta t$ or Euler’s method is no stable for this dynamical system? Use the average velocity $1/2(v_n+v_{n+1})$ to obtain $x_{n+1}$. Are the results any better?
Set $y_0=0$ and $v_x(t=0)=0$. By trial and error find several choices of $x_0$ and $v_y(t=0)$ which yield coonvinient elliptical orbits. Determine total energy, angular momentum, semi-major and semi-nimor axes, eccentricity, and period for each orbit.
You probably noticed that Euler’s algorithm with a fixed $\Delta t$ breaks down if you get to close to the sun. How are you able to visually confirm this? What is the cause of the failure of the method? Think of a simple modification of your program that can improve your results.
In [1]:
class particle2(object):
def __init__(self, mass=1., x=0., y=0., vx=0., vy=0.):
self.mass = mass
self.x = x
self.y = y
self.vx = vx
self.vy = vy
def euler(self, fx, fy, dt):
self.vx = self.vx + fx/self.mass*dt
self.vy = self.vy + fy/self.mass*dt
self.x = self.x + self.vx*dt
self.y = self.y + self.vy*dt
In [111]:
%matplotlib inline
import numpy as np
from matplotlib import pyplot
from matplotlib.colors import ColorConverter as cc
import math
GM = 4*math.pi**2
r = 1. # radius of the orbit
v0 = math.sqrt(GM/r) # This is the condition for circular orbits
x0 = r # initial position
y0 = 0. # we asume we start from the x axis
v0x = 0. # and the initial velocity
v0y = v0 # is perpendicular to the vector position.
dt = 0.01 # time step
tmax = 4.
nsteps = int(tmax/dt)
x = np.zeros(nsteps)
y = np.zeros(nsteps)
vx = np.zeros(nsteps)
vy = np.zeros(nsteps)
energy = np.zeros(nsteps)
x[0] = x0
y[0] = y0
vx[0] = v0x
vy[0] = v0y
energy[0] = 0.5*v0y**2 - GM/r
p = particle2(1., x0, y0, v0x, v0y)
for i in range(1,nsteps):
r = math.sqrt(p.x*p.x+p.y*p.y);
r3 = r * r * r;
fx = -GM*p.x/r3;
fy = -GM*p.y/r3;
p.euler(fx, fy, dt)
x[i] = p.x
y[i] = p.y
vx[i] = p.vx
vy[i] = p.vy
energy[i] = 0.5*(p.vx**2+p.vy**2) - GM/r
t = np.linspace(0.,tmax,nsteps)
pyplot.plot(t, energy, color='green', ls='-', lw=3)
pyplot.xlabel('time')
pyplot.ylabel('Energy');
In [112]:
pyplot.plot(t, x, color='red', ls='-', lw=3)
pyplot.xlabel('time')
pyplot.ylabel('position x');
In [43]:
pyplot.plot(x, y, color='blue', ls='-', lw=3)
pyplot.xlabel('position x')
pyplot.ylabel('position y');
The presence of other planets implies that the total force on a planet is no longer a central force. Furthermore, since the orbits are not exactly on the same plane, the analysis must be extended to 3D. However, for simplicity, we are going to consider a two-dimensional solar system, with two planets in orbit around the sun.
The equations of motion of the two planets of mass $m_1$ and $m_2$ can be written in vector form as $$\begin{aligned} && m_1\frac{d^2 {\mathbf r}_1}{dt^2}=-\frac{m_1MG}{r_1^3}{\mathbf r}_1+\frac{m_1m_2G}{r_{21}^3}{\mathbf r}_{21}, \\ && m_2\frac{d^2 {\mathbf r}_2}{dt^2}=-\frac{m_2MG}{r_2^3}{\mathbf r}_2+\frac{m_1m_2G}{r_{21}^3}{\mathbf r}_{21},\end{aligned}$$ where ${\mathbf r}_1$ and ${\mathbf r}_2$ are directed form the sun to the planets, and ${\mathbf r}_{21}={\mathbf r}_2-{\mathbf r}_1$ is the vector from planet 1 to planet 2. This is a problem with no analytical solution, but its numerical solution can be obtained extending our previous analysis for the two-body problem.
Let us consider astronomical units, and values for the masses $m_1/M=0.001$ and $m_2/M=0.01$. Consider initial positions $r_1=1$ and $r_2=4/3$ and velocities ${\mathbf v}_{1,2}=(0,\sqrt{GM/r_{1,2}})$.
Write a program to calculate the trajectories of the two planets, and plot them.
What would the shape and the periods of the orbits be if the don’t interact? What is the qualitative effect of the interaction. Why is one planet affected more by the interaction that the other? Are the angular momentum and energy of planet 1 conserved? Is the total momentum an energy of the two planets conserved?
In [108]:
NPLANETS = 2
m1 = 0.1
r1 = 1. # radius of the orbit
v1 = math.sqrt(GM/r1) # This is the condition for circular orbits
m2 = 0.01
r2 = 4./3. # radius of the orbit
v2 = math.sqrt(GM/r2) # This is the condition for circular orbits
dt = 0.001 # time step
tmax = 2.
nsteps = int(tmax/dt)
x = np.zeros(shape=(nsteps,NPLANETS))
y = np.zeros(shape=(nsteps,NPLANETS))
vx = np.zeros(shape=(nsteps,NPLANETS))
vy = np.zeros(shape=(nsteps,NPLANETS))
energy = np.zeros(shape=(nsteps,NPLANETS))
x[0,0] = r1
y[0,0] = 0.
vx[0,0] = 0.
vy[0,0] = v1
energy[0][0] = 0.5*v1**2 - GM/r1
x[0,1] = r2
y[0,1] = 0.
vx[0,1] = 0.
vy[0,1] = v2
energy[0][1] = 0.5*v2**2 - GM/r2
planets = [particle2]
for i in range(1,NPLANETS): # create a list of NPLANETS particle2's
planets.append([particle2])
planets[0] = particle2(1., x[0,0], y[0,0], vx[0,0], vy[0,0])
planets[1] = particle2(1., x[0,1], y[0,1], vx[0,1], vy[0,1])
for i in range(1,nsteps):
for n in range(0,NPLANETS):
r = math.sqrt(planets[n].x*planets[n].x+planets[n].y*planets[n].y);
r3 = r * r * r;
fx = -GM*planets[n].x/r3;
fy = -GM*planets[n].y/r3;
planets[n].euler(fx, fy, dt)
x[i,n] = planets[n].x
y[i,n] = planets[n].y
vx[i,n] = planets[n].vx
vy[i,n] = planets[n].vy
energy[i,n] = 0.5*(planets[n].vx**2+planets[n].vy**2) - GM/r
t = np.linspace(0.,tmax,nsteps)
In [109]:
pyplot.plot(x[:,0], y[:,0], color='blue', ls='-', lw=3)
pyplot.plot(x[:,1], y[:,1], color='red', ls='-', lw=3)
pyplot.xlabel('position x')
pyplot.ylabel('position y');
In [110]:
pyplot.plot(t, energy[:,0], color='blue', ls='-', lw=3, label='Planet 1')
pyplot.plot(t, energy[:,1], color='red', ls='-', lw=3, label = 'Planet 2')
pyplot.plot(t, energy[:,0]+energy[:,1], color='green', ls='-', lw=3, label = 'Both planets')
pyplot.xlabel('Energy')
pyplot.ylabel('time');
In this particular case, since both planets do not interact, the individual energies are conserved.
In this section, we will investigate the classical scattering of a particle of mass $m$ by a central potential. In a scattering event, the particle, with initial kinetic energy $E$ and impact parameter $b$ approaches the potential from a large distance. It is deflected during its passage near the force center and eventually emerges with the same energy, but moving at an angle $\Theta$ with respect to the original direction. This problem is very similar in many aspects to the orbital motion, but in this case the potential is repulsive, and it is not necessarily a function of the inverse square of the distance. The energy and momentum are conserved, and the trajectory lies in the plane.
Our basic interest is on the deflection function $\Theta (b)$, giving the final scattering angle $\Theta$ as a function of the impact parameter. This function also depends upon the incident energy. The differential cross section for scattering at an angle $\Theta$, $d\sigma / d\Omega$ is an experimental observable that is related to the deflection function by
$$\frac{d\sigma}{d\Omega}=\frac{b}{\sin{\Theta}}|\frac{db}{d\Theta}|. $$Thus, if $d\Theta /db$ can be computed, the cross section is known.
Expressions for the deflection function can be found analytically only for a few potentials, so that numerical methods usually must be employed. There are some simplification that can me made using the fact that the angular momentum is conserved, which connects the angular and the radial motion, making the problem one-dimensional. However, in this section we are going to use the tools learned in the previous sections, and solve the four first-order differential equations for the two coordinates and their velocities in the $xy$ plane.
In the following we are going to consider a Lennard-Jones potential: $$V(r)=4V_0[(\frac{a}{r})^{12}-(\frac{a}{r})^6], $$ The potential is attractive for long distances, and strongly repulsive approaching the core (see Fig. [lennard-jones]), with a minimum occurring at $r_{min}=2^{(1/6)}a$ with a depth $V_0$.
Before beginning any numerical computation, it is important to have some idea of what the results should look like. Sketch what you think the deflection function should look like at relatively low energies, $E \leq V_0$, where the the peripheral collisions at large $b\leq r_{max}$ will take place in a predominantly attractive potential and the more central collisions will “bounce” against the repulsive core. What happens at much higher energies $E\gg V_0$, where the attractive pocket in $V$ can be neglected? Note that the values of $b$ where the deflection function has a maximum or a minimum, Eq. ([cross]) shows that the cross section should be infinite,as occurs in the rainbow formed when light scatters from water drops.
Write a program that calculates, for a given kinetic energy $E$, the deflection function solving the equations of motion a a number of equally spaced $b$ values between 0 and $r_{max}$.
Use your program to calculate the deflection function for scattering from a Lennard-Jones potential at selected values of $E$ ranging from $0.1V_0$ to $100V_0$. Reconcile your answers in step 1) with the results obtained. Calculate the differential cross sections a function of $\Theta$ at these energies.
If your program is working correctly you should observe for energies $E\leq V_0$ a singularity in the deflection function where $\Theta$ appear to approach $-\infty$ at some critical value of $b$, $b_{crit}$, that depends on $E$. This singularity, which disappears when $E$ becomes larger that about $V_0$ is characteristic of “orbiting”, and the scattering angle becomes logarithmically infinite. What happens is that the particle spends a very long time spiralling around the center. Calculate some trajectories around this point and convince yourself that this is precisely what’s happening. Determine the maximum energy for which the Lennard-Jones potential exhibits orbiting by solving the correct set of equations involving $V$ and its derivatives.
In [ ]: