Today we begin the study of Weber's potential in the isolated two-body system consisting of an electron and positron pair $e^-$ and $e^+$. We assume the particles $e^\pm$ have equal mass $m=m_{e^{\pm}}$. The reduced mass is concentrated at the centre-of-mass $\mu=m/2$.
Weber's force is attractive between the pair $e^\pm$ at all distances.
The particles $e^\pm$ do not indefinitely spiral inwards. Simulations indicate that the radial distance between $e^\pm$ stays strictly bounded between two upper and lower limits $$0 < r_{lower} \leq r \leq r_{upper} < + \infty .$$ This is rigorously proved in Weber-Clemente, 1990.pdf).
If the electron is indivisible particle, then the above two-body problem models a pair $e^-$ and $e^+$ of isolated electron and positron.
But do the particles $e^\pm$ ever 'collide' and annihalate?
In the standard physics textbooks, it seems well known that annihalation between $e^\pm$ occurs and two gamma rays are ejected in opposite directions when $e^\pm$. conserving momentum, etc., and converting all their mass into energy. Thus it's determined that two gamma rays of energy $0.511 keV$ are released, where Einstein's formula $E=m_ec^2$ is applied, where $m_e$ is the reduced mass. [ref] The annihalation of $e^\pm$ is apparently an experimental test of the validity of Einstein's "mass-energy" hypothesis.
But what does Weber's potential say about the annihalation of $e^+$ and $e^-$ ?
If we know the centre of mass has zero net force, then we can replace the positions $r_1$, $r_2$ of the particles by their relative distance $r_{12}$ from the centre of mass. This yields $$r_1=R +\frac{m_2}{m_1+m_2} r_{12}$$ and $$r_2=R -\frac{m_1}{m_1+m_2} r_{12}.$$
Applying Newton's Second Law that $F_{21}=-F_{12}$ yields the following equation for $r_{12}''$: $$ \mu . r_{12}'' = F_{21},$$ where $\mu$ is the reduced mass of the system, namely $\mu=\frac{m_1 m_2}{m_1+m_2}=0.5$.
In the following equations we use numpy.odeint to solve Weber equations of motion of the relative distance $r_{12}$. Therefore we have reduced the two-body problem to a one-body problem. This is a standard reduction.
Given the solution for $r_{12}$, how do we reconstruct the paths/positions of the particles $r_1$, $r_2$ ? Answer: via the relation $r_1=R+\frac{m_2}{m_1+m_2} r_{12}$ and $r_2=R-\frac{m_1}{m_1 + m_2}r_{12}$.
Now the relative distance $r_{12}$ is a type of radial distance, and if $r, \omega$ is spherical coordinates, then we have $$r'^2=|v|^2= x'^2+y'^2+z'^2=(r')^2+r^2 (\theta')^2. $$ The above formula is the usual $|v|^2=v_r^2+v_t^2$, and the tangent velocity $v_t$ satisfies $v_t=r\theta'$, where $\theta'$ is the angular velocity.
The conservation of angular momentum says that the angular moment $L=\mu r\times v$, where $v$ is the linear velocity of $r$, is constant along the motion. Moreover one has $$|L|=\mu r^2 \theta'.$$ Thus we find the formula $$\theta'=\frac{|L|}{\mu r^2}.$$ This implies $$T=\frac{\mu}{2}v^2= \frac{\mu}{2}[(r')^2+\frac{|r\times v|^2}{r^{2}}]$$ represents the kinetic energy of the system.
The conservation of energy says $T+U$ is constant along trajectories.
# Here we define basic functions.
def cross(v1, v2):
x1, y1, z1 = v1
x2, y2, z2 = v2
return [y1*z2 - z1*y2, -(x1*z2 - z1*x2), x1*y2 - y1*x2 ]
def rho(rel_position):
x,y,z = rel_position
return (x*x+y*y+z*z)**0.5
def dot(vector1, vector2):
x1, y1, z1 = vector1
x2, y2, z2 = vector2
return x1*x2+y1*y2+z1*z2
def rdot(position, vector):
return dot(position, vector)/rho(position)
def norm(rel_velocity):
return rho(rel_velocity)
mu=0.5 ## reduced mass of the system. We assume m1 and m2 are equal, hence mu=1/2.
c=1.0 ## speed of light constant in Weber's potential
# Define the angular momentum
def AngMom(rel_position, rel_velocity):
return cross(rel_position, rel_velocity)
def L(rel_position, rel_velocity):
return norm(cross(rel_position, rel_velocity))
# Linear Kinetic Energy
def T(rel_position, rel_velocity):
vt = norm(cross(rel_position, rel_velocity))
# next formula decomposes v^2=(vr)^2+(vt)^2, where vt=r*θ'=|L|/(mu*r)
return (mu/2)*(rdot(rel_position, rel_velocity)**2) + (mu/2)*(rho(rel_position)**-2)*(vt**2)
## Weber Potential Energy
## Negative sign given -1=q1*q2
def U(rel_position, rel_velocity):
x,y,z = rel_position
vx,vy,vz = rel_velocity
rdot=dot(rel_position, rel_velocity)/rho(rel_position)
return -(1/rho(rel_position))*(1-(rdot*rdot)/2)
## import the basic packages
import numpy as np
from scipy.integrate import odeint, solve_ivp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
## Integrating two-body isolated system of oppositely charged particles
## i.e. positron+electron pair.
## The product of the charges q1*q2 is factor in Weber's force law, and appears twice
## in the formula of Newton's F=ma.
def weber(t, rel_state):
x, y, z, vx, vy, vz = rel_state
r=(x*x + y*y + z*z)**0.5
rdot=(x*vx+y*vy+z*vz)/r
A=(-1)*r**-2 ## minus sign from q1*q2
B=1-(rdot*rdot)/2
C=(mu+((c*c*r)**-1))**-1 ## +plus instead of -minus.
dxdt = vx
dydt = vy
dzdt = vz
dvxdt = (x/r)*A*B*C
dvydt = (y/r)*A*B*C
dvzdt = (z/r)*A*B*C
return [dxdt, dydt, dzdt, dvxdt, dvydt, dvzdt]
t_span = (0.0, 100.0)
t = np.arange(0.0, 100.0, 0.1)
y1=[2.0,0,0,] # initial relative position
v1=[-0.4, 0.4, 0] # initial relative velocity
result = odeint(weber, y1+v1, t, tfirst=True) #here odeint solves the weber equations of motion relative y1+v1 for t.
Energy=T(y1,v1) + U(y1,v1)
print('The initial total energy T+U is equal to:', Energy)
print('The initial angular momentum is equal to', norm(AngMom(y1,v1)))
fig = plt.figure()
ax = fig.add_subplot(1, 2, 1, projection='3d')
ax.plot(result[:, 0],
result[:, 1],
result[:, 2])
ax.set_title("position")
ax = fig.add_subplot(1, 2, 2, projection='3d')
ax.plot(result[:, 3],
result[:, 4],
result[:, 5])
ax.set_title("velocity")
What does the above plot demonstrate?
It reveals a precession motion around the centre of mass. This is not predicted by Coulomb's force, which bounds the trajectories to elliptical orbits like Newton's Law of Gravitation.
The force is central, therefore we have conservation of angular momentum, and implies the system is constrained to a plane, namely orthogonal to the angular moment of the system.
Moreover the system satisfies a conservation of linear momentum, namely the sum $T+U$ is constant.
- Problem: Verify Assis-Clemente's 1990 formula for the lower and upper limits of the relative distance along the orbits.
import matplotlib.pyplot
import pylab
r_list=[]
for j in range(1000):
sample_position=[result[j,0], result[j,1], result[j,2]]
sample_velocity=[result[j,3], result[j,4], result[j,5]]
r_list.append(
(
int(j), rho(sample_position) )
)
prelistr = list(zip(*r_list))
pylab.scatter(list(prelistr[0]),list(prelistr[1]))
pylab.xlabel('time')
pylab.ylabel('rho')
pylab.title('Solutions have Upper and Lower Limits')
pylab.show()
#TU_list=[]
#for j in range(1000):
# sample_position=[result[j,0], result[j,1], result[j,2]]
# sample_velocity=[result[j,3], result[j,4], result[j,5]]
# TU_list.append(
# (rho(sample_position),T(sample_position,sample_velocity)+U(sample_position, sample_velocity)))
#prelist1 = list(zip(*TU_list))
#pylab.scatter(list(prelist1[0]),list(prelist1[1]))
#pylab.xlabel('distance r')
#pylab.ylabel('T+U')
#pylab.title('')
#pylab.show()
# The plot below demonstrates the conservation of angular momentum.
# Note that rdot is directly equal to the sample_velocity. I.e. there is
# no need to define rdot=v.hatr/r. This was error.
#A_list=[]
#for j in range(1000):
# sample_position=[result[j,0], result[j,1], result[j,2]]
# sample_velocity=[result[j,3], result[j,4], result[j,5]]
# A_list.append(
# (rho(sample_position), norm(cross(sample_position, sample_velocity)) ) )
#prelist2 = list(zip(*A_list))
#pylab.scatter(list(prelist2[0]),list(prelist2[1]))
#pylab.xlabel('rho')
#pylab.ylabel('Angular Momentum')
#pylab.title('Conservation of Angular Momentum')
#pylab.show()
#rho_list=[]
#for j in range(180):
# rho_list.append(
# (int(j), rho([result[j,0], result[j,1], result[j,2]]),
# )
# )
from sympy import *
t=symbols('t')
m=symbols('m')
c=symbols('c')
r=Function('r')(t)
P=Function('P')(r,t)
F=Function('F')(r,t)
U=-(r**-1)*(1-((r.diff(t))**2)*(2*c*c)**-1)
F=(-1)*(U.diff(t))*((r.diff(t))**-1)
pprint(simplify(U))
print()
pprint(simplify(F)) ## symbolic computation of the Force law.