A Desktop-computer Simulation for Exploring the Fission Barrier

A model of a fissioning nucleus that splits symmetrically both axially and equatorially is used to show how one can predict the presence of a fission barrier of several tens of MeV for nu-clides of mass number A ~ 90 and of ~ 10 MeV for elements such as uranium. While the present model sacrifices some physical realism for the sake of analytic and programming simplicity, it does reproduce the general behavior of the run of fission barrier energy as a function of mass number as revealed by much more sophisticated models. Its intuitive appeal and tractability make it appropriate for presentation in a student level " Modern Physics " class.


INTRODUCTION
In discussions of the theory of nuclear fission in undergraduate texts, one learns that there are two key aspects of the energetics of that phenomenon: 1) That there exists a limit (Z 2 /A) lim ~ 50 beyond which nuclei are unstable against spontaneous fission (Z is the atomic number and A the mass number), and 2) That whether or not a nucleus will fission upon absorbing a neutron depends on whether or not the binding energy so released exceeds a certain fission barrier.This is also known as the activation energy or the threshold energy.In older texts, such discussions were often accompanied by schematic plots of energy versus fission-product separation which showed the barrier as a rise of a few MeV over which the compound nucleus had to 'leap' before fissioning, with the fission products repelling each other and then 'sliding down' the Coulomb potential to a final system energy some 200 MeV lower than the initial spherical configuration [1].The maximum system energy in these illustrations is usually indicated as occurring right at the moment of fission.The fission barrier is a profoundly important business.In the case of uranium, its value (about 6 MeV) amounts to about 1 MeV less than the binding energy released when a nucleus of 235 U absorbs a neutron, but is about 1.5 MeV greater than that in the case of 238 U.This difference is sufficient to render the lighter isotope fissile by neutrons of any kinetic energy whereas 238 U can only be disrupted by those of energy ~ 1.5 MeV or greater.
These fundamental concepts date from a landmark paper published by Niels Bohr and John Wheeler (B&W) in 1939 wherein they modeled nuclei as liquid drops [2].Expressing the shape of a deformed nucleus as a sum of Legendre polynomials configured to conserve volume, they found that the essential aspects of fission could be explained by considering the total energy of a nucleus to be the sum of a surface energy U S proportional its surface area plus an electrostatic contribution U C arising from the Coulombic inter-potentials of the protons.B&W were able to derive a general expression for (Z 2 /A) lim of the form (Z 2 /A) lim = 2(a S /a C ), where a S and a C are surface and Coulomb energy parameters whose values are about 18 MeV and 0.72 MeV, respectively.B&W's algebra involved complicated expansions and integrals of Legendre polynomials and so their work is usually considered to be far too complex for presentation to undergraduates.Many texts skirt the detailed calculations by presenting simplified partial derivations of the Z 2 /A limit by modeling nuclei as ellipsoids and invoking an expression for the self-potential of such a geometry that is itself not trivially derived [3].The full algebra of the (Z 2 /A) lim calculation can be found in a paper recently published by this author [4].Through yet more involved calculations, Bohr and Wheeler also developed an expression for the barrier height, but this applied only for cases where the nuclear deformation represented a small departure form sphericity.In principle, to determine a barrier energy one needs to examine all possible 'deformation trajectories' between an initial spherical configuration and a final two-product configuration and then isolate the minimum excitation energy.
With the rapid development of electronic computers in the postwar era, numerical simulations naturally came to the fore, and models of fissioning nuclei involving up to tenth-order polynomials were being published as early as 1947 [5].The modeling of nuclear masses and barrier energies is now a very complex field of work, where shell effects, pairing corrections, energy levels, and deformation and mass asymmetries all play roles.For example, in a well-known paper, Myers and Swiatecki developed a seven-parameter model for nuclear masses [6].In addition to the usual Coulomb and surface energies, they included a volume-energy term dependent on A, a parity term, and three parameters for describing shell effects; only the Coulomb and surface terms are incorporated in the present simple model.
As a result of such complexities, many texts simply skip over the details of barrier energetics.It would thus appear that a model that could demonstrate the presence of the fission barrier while being tractable at an undergraduate level would be a useful contribution to physics pedagogy.The purpose of the present paper is to develop a straightforward numerical integration scheme for investigating the fission barrier.This scheme is based on a simplified model of fission where the nucleus divides symmetrically into two identical daughter products.In reality, symmetric fission is actually quite uncommon; the most likely mass ratio for the products is about 1.6.However, while a symmetric model does represent some sacrifice of physical reality, it does possess one very important redeeming feature: the shape of the nucleus can be modeled in such a way that its radius of curvature is never discontinuous.The significance of this can be appreciated on recalling that nuclear physicists model surface energies in a manner akin to describing surface tension: a discontinuous curvature would correspond to an infinite force.Despite the simplified approach adopted here, the model proves to give results in fairly respectable agreement with much more sophisticated analyses and so serves its educational purpose.
The model is developed in section 2 below.Section 3 describes the development and results of computer programs written to simulate the fission process.

SYMMETRIC FISSION MODEL
Suppose that our nucleus begins as a sphere of radius R O .As a consequence of some disturbance, it begins to distort in a manner that is at all times both axially and equatorially symmetric, as illustrated in Figure 1.It is assumed that, at any moment, the ends of the distorted nucleus can be modeled as spheres of radius R whose centers are separated by distance 2d and which are connected by an equatorial 'neck'.The outer edge of the neck is taken to be defined by a circular arc of radius R which is part of a circle that just tangentially touches the spheres which constitute the product nuclei; imagine rotating the figure around the central vertical axis.As described above, modeling the fission process in this way allows us to avoid introducing any curvature discontinuities into the shape of the distorted surface; this would be very difficult to do if the process were not both axially and equatorially symmetric.As the spheres move apart, R must decrease to conserve nuclear volume as nuclei are considered to be incompressible fluids.The radius of the neck will decrease and eventually reach zero, at which point the teardrop-shaped products will break contact and fly apart due to mutual repulsion.A convenient coordinate for parameterizing the process is the angle θ, which is measured counterclockwise from the polar axis to a line joining the centers of the upper sphere and the imaginary one that defines the equatorial neck.θ decreases from π/2 at the start of the process to π/6 at the moment of fission.

Volume and Surface Areas; Volume Conservation
The volume of the partially-fissioned nucleus is given by [7]   where Demanding conservation of volume, Eqs.1 and 2 give the radius R in terms of the initial radius R O and θ as The center-to-center separation of the spheres at any moment is The surface area S of the distorted nucleus is given by where At the moment when the equatorial neck has shrunk to zero radius (θ= π/6), the radii of the spherical caps is The separation of the centers at this time will be 2 2 3 2.737 , and the full height of the distorted nucleus will be 2R fiss + 2d fiss ~ 4.317R O .

Surface Energy
The surface energy of the deformed nucleus is assumed to be directly proportional to its surface area S at any moment.To formulate this in terms of atomic and mass numbers, it is helpful to invoke the empirical result that nuclear radii behave as R ~ a o A 1/3 , with a o ~ 1.2 fm.Invoking this result along with Eqs.3, 5, and 6 gives the surface area as where S  is purely a function of the emergence angle θ.If we introduce  as the value of the conversion factor from surface area to energy in MeV, the product 2 4π o a  is the surface energy parameter a S ~ 18 MeV discussed in the Introduction above.We can then write the surface energy as

Coulomb Energy
The Coulomb energy of the deformed nucleus of Figure 1 is determined by direct numerical integration: The fissioning nucleus is imagined to be situated within a lattice of volume elements, and the self-energy is com-puted by adding up the Coulomb energies of all pairs of elements that find themselves within the nucleus.While this procedure is conceptually straightforward, there are various practical issues to deal with, such as ensuring that the number of cells is sufficiently large to obtain reasonable accuracy without rendering the computation prohibitively time-consuming, taking advantage of the axial symmetry of the situation to minimize the number of computations, and adopting a convenient system of units.
It was decided to imagine surrounding the configuration of Figure 1 with a cylindrical lattice of N cells comprising N Z vertical layers, each of which consists of N R radial rings and N  angular wedges: The lattice remains of constant volume throughout the integration (as do the individual cells), and is set to have a radius equal to that of the initial undeformed nucleus (R O ) and a height just great enough to accommodate the just-fissioned system, h = βR O , where β is the factor 4.317 … discussed just after Eq.7.The bottom of the lower sphere is taken at every step in the calculation to be sitting at (x, y) = (0, 0).R O is used as the unit of distance for the numerical integration.The radial rings are set up to become thinner with increasing radius in order that all cells are of the same volume; the radius of the i'th ring (in units of R O ) is given by i R r i N  .Suppose that the nucleus contains Z protons.Its charge density will be , and, from the radius, height, and number of cells in the cylinder as described above, the charge contained in a cell that lies within the nucleus will be Q = 3eZβ/4N.The Coulomb potential between any two 'occupied' cells (i, j) separated by distance r ij will be 4π and summing over all pairs of cells, the total Coulomb potential energy emerges as where a C is the Bohr-Wheeler Coulomb energy parameter and where where δ i = 1 if cell i lies within the nucleus, and zero if it does not.Like S  , C  is purely a function of the emergence angle θ.Combining Eqs.10-13, the total (surface + Coulomb) configuration energy of the nucleus becomes With a S = 18 MeV and a o = 1.2 fm, 15a C /16a S ~ 0.0375.These values are assumed in what follows.

PROGRAMS AND RESULTS
Two double-precision FORTRAN programs, ANGLE and FISSION, have been developed to carry out the calculations described above.Since the geometry of the process depends only on the emergence angle θ, it is convenient to have one program (ANGLE) generate a tabulation of R, d, and S  as functions of θ, and a second one (FISSION) carry out the numerical integration for C  .The programs need only be run once, after which the run of total energy for any (A, Z) as a function of θ can be obtained by scaling S  and C  according as Eq.14.ANGLE is written to step  in increments of −0.01 radians from π/2 at the start of the process down to π/6 at the moment of fission.FISSION has been written to take advantage of the axial symmetry of the model.The internal and inter-pair potentials of the N  angular wedges need not be computed for all individual and pairs of wedges since the internal potentials of all wedges will be identical and the potentials between all possible pairs of wedges can be determined by computing the potential between any one and half of all the others.The equatorial symmetry of the model was not built into the program in the anticipation that it might later be used to simulate non-symmetric fissions.After some experimentation, it was found that choosing (N Z , N R , N  ) = (200, 100, 100) provided for both reasonably expedient run times and sensible accuracy; a complete run requires about 12 hours on a Macintosh G5 and typically returns a starting-configuration energy accurate to better than 0.1% (see below).
Figure 2 shows results for (A, Z) = (90, 40) and (235, 92), that is, nuclides of zirconium and uranium, respectively.The ordinate is the change in total configuration energy in the sense (deformed nucleus -original nucleus).Zirconium is used as an example because for real nuclei, much more sophisticated models, such as that of Myers and Swiatecki [6], reveal that fission barriers peak at a value of ~ 55 MeV for nuclei with A ~ 90.While there is some noise in the simulation due to the finite number of lattice cells, the computed maximum of ~ 47 MeV for zirconium is in good agreement with this peak value.For uranium, the computation gives a barrier of ~ 13 MeV, somewhat high compared to the true value of ~ 6 MeV but still respectably accurate considering the simplicity of the model.At the moment of fission, ∆E in the case of uranium has dropped to about −45 MeV.One feature to note in Figure 1 is that in the case of uranium the maximum distortion energy is reached at about the middle of the fission process, whereas for zirconium the maximum occurs at very near the end of the process.
The numerical precision achieved with the model can be judged by examining how well it reproduces the configuration energy for the initial undeformed nucleus.This can be computed analytically as For (A, Z) = (235, 92) and (a C , a S ) = (0.72 MeV, 18 MeV), this gives U start = 1673.0MeV; the simulation gives 1673.74MeV, only 0.04% high.
Figure 3 shows how the results of the present model compare to those Myers and Swiatecki's seven-parameter model.For the present model, this plot was formed by computing ∆E() for values of A = 10, 15, 20, •••, 300 and searching for the maximum value of ∆E for each A; the curve is interpolated.In computing these maximum ∆E values it was assumed that the atomic number Z for a given value of A could be modeled as Z ~ 0.627A 0.917 (1 < Z < 98), a purely empirical result based on the (Z, A) trend of 352 nuclides with half-lives > 100 years [8].The Myers and Swiatecki curve was generated by fitting the empirical formula  successfully reproduce the overall trend of ∆E(A).
The key conclusions here are that the simulation is physically appealing, convincingly demonstrates the existence of fission barriers, reproduces the trend of barrier heights with mass number, and gives plausible values for the barrier heights for interesting nuclides like uranium.

SUMMARY
A model for simulating symmetric fission with a desktop computer has been developed.This model convincingly demonstrates both the presence of and approxi-mately correct magnitudes for fission barriers of real nuclides.The only physical concepts involved are Coulomb and areal self-energies, and, as such, the model is presentable to undergraduates in order to help them understand this important phenomenon.A spreadsheet giving results for S  and C  as functions of θ which users can employ to run calculations for any desired (A, Z) values is freely available for download [9].
16) to Figure 18 of their paper by a least-squares calculation.The values of the parameters are (C, p, A o , n) = (5.11517, 0.592 21, 156.630 3, 2.537 96).This fit predicts a slightly high peak value for the curve (~ 58 MeV vs. a true value of ~ 55 MeV), but is convenient for comparing the models.The present model predicts somewhat high values of E at both low and high values of A and somewhat low ones for intermediate values, but it does

Figure 3 .
Figure 3.Comparison of results of present model (dashed line) with those of Myers and Swiatecki (solid line; Ref. [6]).See discussion of Eq. 16.