Hover over the menu bar to pick a physics animation.

Canvas 0: 2D Parabolic Potential Color Plot

Extra Width: 30

Maximum Potential: 0.5

`V_("electric")=(Z_(alpha)Z_N e^2)/(4pi epsilon_0 r)`

For this animation I will model the combined potential as a truncated parabola where the parabola is centered at the center of the original nucleus as shown in Canvas 0.

I also made an estimate (see spreadsheed 0) of the starting separation of the alpha particle from the center of the nucleus to provide the usual measured alpha particle kinetic energy which is in the range of 3-5 mega-electron-volts (MeV). The spreadsheet assumes that the alpha particle starts at the center of the potential with some kinetic energy. Note that we must distinguish between the atomic number (`Z_N`), which is the number of protons in the nucleus, and the mass number (`M_N`), which is the sum of the mass of the protons and neutrons in the nucleus. In a large nucleus the electric potential is proportional to `Z_N`.

The strong force decreases with distance as a negative exponential called the Yukawa potential. Thus, if the atomic number is larger than that of lead (`Z_N=82`) the nuclear radius is large enough that the force of the electric charges can be greater than the strong force and the nucleus can sometimes break up (fission). In fact the nuclear binding energy per nucleon peaks at a mass number of 56 (iron) and drops by about 10 % by a mass number of 200.

One might ask why some heavy nucleii are stable and others undergo fission. These are isotopes like lead 208 which have closed outer nuclear shells similar to the closed outer electron shells of the noble metals. These isotope numbers are called "magic numbers".

Spreadsheet 0:Alpha Emission Energy Calculation

Canvas 1: 2D Parabolic Potential Color Plot

with Wave Packet Superimposed

with Wave Packet Superimposed

Initial Kinetic Energy: 0.1

Initial Wave Vector Angle: 0

Initial Potential Energy: 0.1

Initial Particle Angle: 90

Plot Wave Only

Plot Wave over Potential

Plot Real

Plot Imaginary

Plot Magnitude

Iteration Time Increment: 0.2

`1/(2m) nabla cdot nabla Psi+V Psi=i h/(2pi) (d Psi)/dt (1)`

where m is the particle mass, chosen to be 1 here, `nabla *` is the divergence operator, ` nabla` is the gradient operator, `i=sqrt(-1)` V is the potential function, V(x,y), and `h` is Planck's constant, also chosen to be 1. The spatial units used here will be screen pixels and the time units will be animation update units, usually 1/60 of a second so that we see 60 frames per second. Equation 1 will be solved by iterating the values of `Psi(x,y,t)` in both space and time. The space iteration is done by using the finite difference method with higher order terms. The usual finite difference (FD) method for a 1D function vector `f_i` multiplies its values by a simple matix:`1/2 nabla cdot nabla f=((-1,1/2,0,0,0),(1/2,-1,1/2,0,0),(0,1/2,-1,1/2,0),(0,0,1/2,-1,1/2),(0,0,0,1/2,-1))((f_1),(f_2),(f_3),(f_4),(f_5)) (2)`

What will be used here is a matrix that has higher order terms than the simple `(1/2,-1,1/2)` surrounded by zeros that are shown in example 2. These give a better approximation to the second derivative needed for equation 1. For, 3rd, 5th and 7th order approximations, the coefficients are:`1/2 c_3=[1/2,-1,1/2]`

`1/2 c_5=[-1/24,4/6,-5/4,4/6,-1/24]`

`1/2 c_7=[1/180,-3/40,3/4,-49/36,3/4,-3/40,1/180]`

Typical `psi` vector arrays are 300 elements long. The FD matrix would be fairly large but it is zero except around its diagonal (called "sparse") so it still iterates reasonably fast on a desktop computer but slower on a mobile device.For the time derivative, it is sufficient to use a simple iteration scheme where the time increment, `dt` is chosen so that the integration remains stable. To improve on this, one could use a 4th order Runge-Kutta method.

The 2D "vector" is actually a matrix where the row variation is the vertical part of the 2D grid and the colums are the horizontal part of the grid. Example 3 shows how this works for the (x,y) coordinate system and a 7x7 grid:

`Psi_(7x7)(x_i,y_i)= ((psi(x_1,y_1),psi(x_2,y_1),psi(x_3,y_1),psi(x_4,y_1),psi(x_5,y_1),psi(x_6,y_1),psi(x_7,y_1)), (psi(x_1,y_2),psi(x_2,y_2),psi(x_3,y_2),psi(x_4,y_2),psi(x_5,y_2),psi(x_6,y_2),psi(x_7,y_2)), (psi(x_1,y_3),psi(x_2,y_3),psi(x_3,y_3),psi(x_4,y_3),psi(x_5,y_3),psi(x_6,y_3),psi(x_7,y_3)), (psi(x_1,y_4),psi(x_2,y_4),psi(x_3,y_4),psi(x_4,y_4),psi(x_5,y_4),psi(x_6,y_4),psi(x_7,y_4)), (psi(x_1,y_5),psi(x_2,y_5),psi(x_3,y_5),psi(x_4,y_5),psi(x_5,y_5),psi(x_6,y_5),psi(x_7,y_5)), (psi(x_1,y_6),psi(x_2,y_6),psi(x_3,y_6),psi(x_4,y_6),psi(x_5,y_6),psi(x_6,y_6),psi(x_7,y_6)), (psi(x_1,y_7),psi(x_2,y_7),psi(x_3,y_7),psi(x_4,y_7),psi(x_5,y_7),psi(x_6,y_7),psi(x_7,y_7))) (3)`

To obtain the 5th order second derivative required for equation 1 for the third row, `y_3`, of example 3, we would successively multiply by `c=1/2 c_5` as shown below:`1/2 (del ^2 Psi(x,y_3))/(del x^2)=c_1 psi(x_2,y_3)+c_2 psi(x_3,y_3)+c_3 psi(x_4,y_3)+c_4 psi(x_5,y_3)+c_5 psi(x_6,y_3) (4)`

and there would be a similar multiplication for 7th order accuracy but that requires more matrix rows than are shown in example 3. For, say, the top `y_1` row of example 3, the multiplication by `c_i` in equation 4 is truncated on the left to obtain:`1/2 (del ^2 psi(x,y_1))/(del x^2)=c_3 psi(x_1,y_1)+c_4 psi(x_2,y_1)+c_5 psi(x_3,y_1) (5)`

Similarly, for the bottom `y_7` (7th) row, the multiplication is truncated on the right:`1/2 (del ^2 psi(x,y_7))/(del x^2)=c_1 psi(x_5,y_7)+c_2 psi(x_6,y_7)+c_3 psi(x_7,y_7) (6)`

Of course, for `(del ^2 Psi(x_4,y))/(del y^2)` we would apply the same coefficients to the 4th column and get an expression similar to equation 4. More explicitly the complete second derivative `Psi''(x_4,y_3)` is
`1/2 nabla cdot nabla(x_4,y_3)=c_1 psi(x_2,y_3)+c_2 psi(x_3,y_3)+c_3 psi(x_4,y_3)+c_4 psi(x_5,y_3)+c_5 psi(x_6,y_3)`+

`c_1 psi(x_4,y_2)+c_2 psi(x_4,y_3)+c_3 psi(x_4,y_4)+c_4 psi(x_4,y_5)+c_5 psi(x_4,y_6) (7)`

`h/(2pi) (d Psi_(imag))/dt=-[-1/(2m) nabla cdot nabla Psi_(real)+V Psi_(real)]`

`Psi_(imag)(t+dt)=Psi_(imag)(t)-2pi/h[-1/(2m) nabla cdot nabla Psi_(real)+V Psi_(real)]dt`

`h/(2pi) (d Psi_(real))/dt =-1/(2m) nabla cdot nabla Psi_(imag)+V Psi_(imag) `

` Psi_(real)(t+dt)=Psi_(real)(t)+2pi/h[-1/(2m) nabla cdot nabla Psi_(imag)+V Psi_(imag)]dt (9)`

If you choose an initial kinetic energy that is large enough and/or a maximum potential that is small enough,
you will see that the truncated parabolic potential no longer contains the wave packet and it breaks up. When that happens, it
runs into boundary of the containing square box. Since the boundary of the square box represents an infinite potential change, the wave packet is
reflected from the boundary and may re-enter the parabolic potential only to leave again on the other side.