Does Classical Physics Predict Nuclear Fusion?

Introduction

The occurrence of nuclear reactions in the Sun is invoked to account for the huge energy output of our home star. However, for the reactions to occur the nuclei in the Sun must have enough energy to overcome the electrostatic repulsion between the positively-charged nuclei (this repulsion is called the Coulomb barrier). If one treats the solar core as a hot, ideal gas then one can use the equipartition theorem to calculate the average energy of the nuclei in the gas. For typical stellar core temperatures of about two million kelvins one finds the kinetic energies of the nuclei in the gas are far below the Coulomb barrier. Quantum mechanical tunneling can explain the occur of reactions below this barrier, but classical physics offers an alternative explanation. For the constituents of an ideal gas the distribution of velocities in the gas can be described by

                                                  n(v) = 4π N((((m/(2πkT))))^(3/2))(v^2)(e^(-(mv^2)/2kT))

where v is the velocity of a gas particle,  m is the mass of the individual gas constituents, N is the number of particles, T is the temperature, and k is Boltzmann's constant. Even if the average kinetic energy of the particles is well below the Coulomb barrier there is still a finite probability of two nuclei colliding with sufficient energy to surmount the barrier. In this laboratory you will determine if this classical argument is a viable alternative.

The Maxwellian Velocity Distribution

The expression above is called the Maxwellian velocity distribution. Use it to plot the distribution of velocities for an ideal gas of protons at a temperature of two million kelvins. Assume that N=1 so the function behaves as a normal probability. How would you test that your expression is correct? What is the most probable velocity for a proton to have?

The Coulomb Barrier

The proton can be treated as a sphere with a radius of about 1 fm. The interaction between two protons is strictly via the Coulomb force until they are touching and then they fuse. What is the height in energy units of the Coulomb barrier between the two protons in the center of mass frame of reference?  The relationship between the energy in the center of mass frame and the laboratory frame of reference (where one of the particles is stationary) is

                                                                   (Subscript[E, cm])  = ((Subscript[m, tgt])/((Subscript[m, proj]) + (Subscript[m, tgt])))(Subscript[E, lab])

where (Subscript[m, proj]) is the mass of a moving particle (the projectile) and (Subscript[m, tgt])is the mass of a stationary particle (the target). Treat one of the protons as a projectile and the other as a stationary target. What is the velocity of a proton projectile with enough energy to overcome the Coulomb barrier in the center-of-mass frame? How does this velocity compare with the most probable energy at the stellar core you calculated above?

The Classical Fusion Probability

We want to determine the probability of a proton in our stellar core having enough energy to overcome the Coulomb barrier between two protons. To find this we have to integrate the probability distribution of the protons from a lower limit derived from the Coulomb barrier out to infinity. This lower limit is the minimum velocity a proton can have incident on a stationary proton and still pass the Coulomb barrier. 

A sample of the Mathematica commands you need to perform this task is shown below. Notice the Integrate command is 'wrapped' by the N[ ] command. This tells Mathematica to print the final numerical result instead of keeping it in symbolic form.

f[v_] := v*Exp[-v^2]
lowV = 10;
N[ Integrate[f[v], {v,lowV, Infinity}] ]

All computers have a lower limit to the smallest number they can calculate. If you get a result of zero for one of your calculations you can determine the computer limit with the command

$MinMachineNumber

If you obtained zero for a calculation it may mean your result was below this lower limit.

Depending on what version of Mathematica you are using, the transmission coefficient may be treated as a complex number even if the imaginary part is zero. This means you can’t plot it or manipulate it like you would a real number. To fix this use the Re function that will take only the real part of your value.

myT = 1.0 + I*2.0
Re[myT]

The Quantum Mechanical Result

We now have to compare our classical result with the quantum mechanical one. For now we will make a very crude estimate of the transmission coefficient for two colliding protons by treating the potential energy curve as a rectangular barrier. We will assume that if the proton 'tunnels' through the barrier fusion will occur. Make a plot of the Coulomb potential energy between two protons as a function of the distance r between their centers and choose the height and width of the rectangular barrier to approximate the Coulomb potential. Later on, we will make a more accurate approximation of the Coulomb potential energy.

To make a plot of the rectangular barrier you can use the Which function to construct a discontinuous function.

V0=10;
a0=1.0;
Vrec[x_] := Which[x<0,0,0≤x≤ 2*a0, V0, x>2*a0,0]
Plot[Vrec[x],{x,-a0,4*a0},PlotRange->{-0.1,1.5*V0},Frame->True,FrameLabel->{"x","V(x)"}, PlotStyle->RGBColor[1,0,0],Exclusions->None]

We found in class that the transmission coefficient is 

                                                                      T = (1/(|(Subscript[t, 11])((|)^2)))

where

                    (Subscript[t, 11]) = (1/4)[ (1 + ((Subscript[k, 2])/(Subscript[k, 1])))(1 + ((Subscript[k, 1])/(Subscript[k, 2])))(e^(-(Subscript[ik, 2])2a)) + (1 - ((Subscript[k, 2])/(Subscript[k, 1])))(1 - ((Subscript[k, 1])/(Subscript[k, 2])))(e^((Subscript[ik, 2])2a))]
                      
and

                                               (Subscript[k, 1]) = (Sqrt[((2mE)/(ℏ^2))])      and      (Subscript[k, 2]) = (Sqrt[((2m(E - (Subscript[V, 0])))/(ℏ^2))])

To calculate the transmission coefficient (with its complex conjugate) in Mathematica use the command

Trans[En_] := 1/(Conjugate[t11[En]]*t11[En])

Where En is the energy of the projectile and Trans[En] is the transmission coefficient. 

How can you check your calculation? Use that method to check your calculation. What is the probability of tunneling occurring at the classical most probable energy in your gas at the stellar core?

Conclusion

Which mechanism is a more likely candidate to yield a higher reaction rate. Explain. In either case is the reaction probability very high? Are their any flaws in either of these models? How would you rectify them?