nuclear_binary_system_EJS.xml - by Javier E. Hasbun and Benjamin E. Hogan (2012)

For a binary system, such as the one shown below, we can use the following equations to determine the motion of m1 and m2:

r1 = rcm - (m2/M) r
r
2 = rcm + (m1/M) r
r
cm = (m1r1+m2r2)/M
μr'' = f r/r
μ = m1m2/m1+m2

Since we know that

μr'' = f r/r

we can solve for r'' giving

r'' = (f/μ) r/r

Seperating this into x- and y-components gives

d2x/dt2 = (f/μ)*x/√(x2+y2)

d2y/dt2 = (f/μ)*y/√(x2+y2)

Next, we need to find the force. Using the WOODS-SAXON spherically symmetric potential which is given by

V(r) = -V0*ff(r,R,a)

where ff(r,R,a) is a Fermi-function form factor given by

ff(r,R,a) = [1+exp(r-R/a)]-1,

and the parameters are

Vo = 52.06 MeV

a = 0.662 fm

R = R0*A^(1/3), Ro = 1.26 fm.

Since we know that a potential is defined as V(r) = ∫ f(r)∙dr,

we can find f(r) by taking the derivative of V(r) w.r.t. r. We get

f(r) = -V0*exp(r-R/a)/{a*[1+exp(r-R/a)]}.

This is then inserted into the equations of motion giving

d2x/dt2 = -V0*exp((r-R)/a)*x/{a*μ*[1+exp((r-R)/a)]*√(x2+y2)}

d2y/dt2 = -V0*exp((r-R)/a)*y/{a*μ*[1+exp((r-R)/a)]*√(x2+y2)}

Solving these equations give the x- and y-coordinates of the center of mass which in turn give the x- and y-coodinates of m1 and m2, respectively.

To ensure unit compatibility, let

x = X*a0

y = Y*a0

a = A*a0

R = R*a0

t = T*tau

V0 = V0*Eb

m1 = M1*m0

m2 = M2*m0

μ = μ*m0

The equations of motion then become

d2X/dT2 = -V0*exp((sqrt(X^2+Y^2)-R)/A)*X/{A*μ*[1+exp((sqrt(X^2+Y^2)-R)/A)]^2*√(X2+Y2)}

d2Y/dT2 = -V0*exp((sqrt(X^2+Y^2)-R)/A)*Y/{A*μ*[1+exp((sqrt(X^2+Y^2)-R)/A)]^2*√(X2+Y2)}

after arbitrarily setting (Eb*tau^2 / m0*a0^2) =1.

Solving the above expression for tau, we see that the one second is equivalent to 1.02E-22 seconds of simulation time.

To determine the initial conditions, we set t=0, rCM=0, and drCM/dt=vCM=0.

This gives

r1,0 = -(m2/m1)*r2,0

and

v1,0 = -(m2/m1)*v2,0 ,

or in component form:

x1,0 = -(m2/m1)*x2,0

y1,0 = -(m2/m1)*y2,0

vx1,0 = -(m2/m1)*vx2,0

vy1,0 = -(m2/m1)*vy2,0


Setting either component of x=0 makes the y-component a maximum. For y we choose 1.

Similarly this would make the y-component of velocity zero, the the x-component would be a maximum.

For vx we chose 0.085.