Alpha Particle Emission from a Large Nucleus

Alpha particle emission is a very frequent mode of disintegration (fission) of very large unstable nucleii. For example, uranium-238 may decay to thorium-234 which yields a nucleus (alpha particle) of mass 4 which is the same as the helium 4 nucleus. Electric charge on the alpha particle is not necessarily the same as
Canvas 0: 2D Parabolic Potential Color Plot
Extra Width: 30
Maximum Potential: 0.5
that of the helium nucleus (-2e). The decay is by the way of quantum tunneling through the combined strong force attractive potential, V, of the large nucleus and the much weaker repelling potential of due to the combined electric potentials of the alpha particle (`Z_(alpha)=2`) and the residual charge (`Z_N=N-2`) where N is the original charge (atomic number) of the large nucleus.

`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
We will compute the trajectory of a classical particle which has the same initial position and momentum as the quantum particle and show that both tend to move as expected in the 2D parabolic potential. The difference will be that the quantum particle will have some fraction of its probability density outside the range of the classical particle. This is called tunneling and is the same as occurs in optics when the total internal reflection (TIR) is frustrated by having a second surface of the same index of
Canvas 1: 2D Parabolic Potential Color Plot
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
refraction within a wavelength of the TIR surface.

Numerically Propagating the Wave Packet

The quantum wave function that we will propagate is really a linear combination of many intrinsic eigenfunctions (eigenvectors) of the parabolic potential. Since it is not a single eigenvector, we will henceforth call it a wave packet. The quantum particle wave packet `psi` is computed by time and space iteration of the Schrodinger equation:

`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. A mobile device is adequate for social use but is NOT an acceptable device on which to animate physics in a detailed manner.
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)`

where this is the complete second partial derivative at coordinates `(x_4,y_3)`. Note from equation 1 that `Psi` has both real and imaginary parts. The time derivative on the right side of the equation feeds a component that is 90 degrees out of phase into the left side of the equation. The complex nature of the equation actually can be an advantage that improves the accuracy of the iteration. The present animation owes its speed and accuracy to an algorithm that iterates the imaginary part first and then uses the modified complex wave packet to iterate the real part. This was shown to reduce many errors in the following paper: P.B. Visscher, A Fast Explicit Algorithm for the Time-Dependent Schrodinger Equation, Computers In Physics, 596–598 (Nov/Dec 1991). The method results in significant reduction of error terms. Without this method the iteration of `Psi` would tend to diverge at much shorter iteration time increments. To implement this method, we do the following steps:

`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`

where dt is the time increment. Then using the newly updated `Psi_(imag)`

`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)`

As a starting wave packet setup I have chosen one that will make circular orbits about the center of the axisymmetric potential. This particular setup does not yield an emission of the alpha particle represented by the wave packet, although a circular one could if given a large enough k vector. However one might expect the wave packet size to remain constant as it propagates around its orbit. That is not the case here. The wave packet greatly diminishes in size at 90 degrees and 270 degrees. Even while it diminishes in size, the integral of `PsicdotPsi` remains constant as reported in the upper left of Canvas 1. The reason for that is that the wave packet is not chosen to be a simple sum of the wave eigenfunctions of the axisymmetric parabolic potential. It is still a legitimate wave packet, however, since the eigenvectors of the potential are a completed set and therefore any wave packet shape can be constructed from them. To implement a starting wave packet that will remain constant in size around its orbit, I may need to change the coordinates from Cartesian to polar. That will be a future animation under "Quantum Physics" menu item.

Results

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.
For this version I chose to initialize the wave packet as a curved gaussian where the curvature is that of a circular orbit with the given kinetic energy. This results in 4 fold periodicity for the circular orbit.
I also computed the classic orbit using polar coordinates. This is rather involved and I referred heavily to 1

Summary

The learner is welcome to press F12 (Windows) or Option-Command-I (Mac OS) and then select Sources to view the code that I used for this animation. Also the learner may not be familiar with all the terms I have used but she obviously has an internet connection so she can search for and read very good explanations of those terms.