Mass distribution of fission fragments within the Born-Oppenheimer approximation

The fission fragments mass-yield of U is obtained by an approximate solution of the eigenvalue problem of the collective Hamiltonian that describes the dynamics of the fission process whose degrees of freedom are: the fission (elongation), the neck and mass-asymmetry modes. The macroscopic-microscopic method is used to evaluate the potential energy surface. The macroscopic energy part is calculated using the liquid drop model and the microscopic corrections are obtained using a Woods-Saxon single-particle levels. The four-dimensional modified Cassini ovals shape parametrization is used to describe the shape of the fissioning nucleus. The mass tensor is taken within a cranking-type approximation. The final fragment mass distribution is obtained by weighting the adiabatic density distribution in the collective space with the neck-dependent fission probability. The neck degree of freedom is found to play a significant role in determining the final fragment mass distribution.


Introduction
Within the last few decades a number of theoretical approaches for the description of fission process were developed. They can be divided into two groups: timedependent and quasi-static approaches. Among the timedependent approaches one could mention the time-dependent density functional theory [1], the time-dependent generator coordinate method [2], the Brownian shapemotion approach for nuclear fission [3], the models based on the solution of Langevin equations for the shape degrees of freedom [4][5][6][7]. The quantum mechanical timedependent approaches, like time-dependent density functional theory (TDDFT) [8], are unfortunately very time consuming. In order to reduce the computation time to a manageable level, various simplifications are introduced. By now the results of the time-dependent approaches are available for a very small number of fissioning nuclei. The quasi-static methods like the scission point model [9][10][11] are simpler and allow the calculations for many fissioning systems. However in the scission point model, it is difficult to define in a formal way (free of the fitting parameters) the scission point configuration. So, any new approach to the description of fission process would be quite welcome.
A very stringent test of any theoretical model which describes the nuclear fission process should be a proper a e-mail: pomorski@kft.umcs.lublin.pl reproduction of the fission fragments mass distribution. The goal of the present paper is to obtain such a distribution by an approximate solution of the eigenvalue problem of the 3-dimensional collective Hamiltonian with degrees of freedom corresponding to elongation, neck formation and mass asymmetry of the nuclear shape. The present model is similar to the 2-dimensional one of refs. [12,13] but the non-adiabatic and dissipative effects are not taken here into account since their effect is rather small for lowtemperature fission.
The potential energy surface (PES) is obtained in the present work using the macroscopic-microscopic method with the liquid drop model for the macroscopic part of the energy while the microscopic shell and pairing corrections are calculated using the Woods-Saxon (WS) singleparticle levels [14,15]. The shape of the fissioning nucleus is described by the four-dimensional modified Cassini ovals (MCO) [14,16]. It was shown in ref. [17] that the MCO describe very well the so-called optimal nuclear shapes obtained through a variational description [18], even those close to the scission configuration.
The mass tensor is taken within the cranking-type approximation (cf. e.g. sect. 5.1.1 of ref. [19]). The Born-Oppenheimer approximation (BOA) is used to describe the coupling of the fission mode with the neck and mass asymmetry degrees of freedom. It will be shown that, in order to obtain a fission fragment mass distribution in agreement with the experimental data, the fission probability should depend on the neck size. The paper is organized in the following way. First we shortly present the details of our theoretical model, then we show the collective potential energy surface evaluated in the macroscopic-microscopic model for 236 U and the components of the mass tensor. The calculated fission fragments mass distribution is compared with the experimental data in the next section. Conclusions and possible extensions and applications of our model are presented in the summary.

Shape parameterization
We define the shape of fissioning nucleus by the parameterization developed in [14], where some cylindrical coordinates {ρ, z} related to the lemniscate coordinates system {R, x} by the equations are introduced. The relation between {ρ, z} and {R, x} depends on the parameter ε ≡ s/R 2 0 , where s is the squared distance between the focus of Cassinian ovals and the origin of coordinates. The surface in the lemniscate coordinates system is obtained by the establishing relation between R and x, R = R(x). The coordinate surfaces R(x) = R 0 are the Cassini ovals (shown in the bottom row of fig. 1).
The deviation of the nuclear surface from Cassini ovals is defined by an expansion of R(x) in series of Legendre polynomials P n (x), where R 0 is the radius of the spherical nucleus. The cylindrical coordinates {ρ, z} are related to {ρ, z} by where z cm is the z-coordinate of the center of mass of Cassini ovaloid (2) and the constant c is introduced in order to insure the volume conservation. Instead of ε, it turns out convenient, see [15], to introduce another parameter, α, where z L (z R ) is the minimum (maximum) value of the cylindrical coordinate z and r neck is the radius of the cross section of the nuclear surface in the plane z = 0. The explicit relation between α and ε can be derived from (4) using the connection of the {R, x} coordinates with the cylindrical ones, The parameter α is defined so that at α = 1 the neck radius turns into zero for any value of all other deformation parameters α n . The value α = 0 corresponds to the spherical shape, α = 1 to the shape with zero neck. In the case that all parameters α n are equal to zero, α coincides with ε.
The coefficients α and α n are considered as the deformation parameters. Examples of the shapes in the Cassini parametrization are shown in fig. 1.

Born-Oppenheimer approximation
Below we use the following collective coordinates to describe the fission dynamics: Here R 12 is the distance (in units of the radius R 0 of the corresponding spherical nucleus) between the mass center of left and right parts of the nucleus. Within the Cassini parameterization it is convenient to divide the shape into left and right parts by the position of the point x = 0 (x is one of the two lemniscate coordinates, see (1)). Such a definition makes sense both for pear-like shapes and for the shapes with a neck. In the last case the point x = 0 coincides with a good accuracy with the position of the neck. The a in (6) is the mass asymmetry coordinate, V 1 and V 2 are the volumes of left and right parts of the nucleus. The parameter α 4 describes the neck degree of freedom (see ref. [16]).
With these coordinates the classical energy of the system becomes where M ij and V ({q i }) denote the inertia tensor and the potential energy, respectively. The quantized form of this Hamiltonian is the following: The eigenproblem of this Hamiltonian will be solved here in the Born-Oppenheimer approximation (BOA) in which one assumes that the motion towards fission is much slower than the one in the two other collective coordinates. The effect of the nonadiabatic terms on the fission massyield is rather small at low-energy fission as it was found in ref. [13] in a 2D model and this is the reason why we shall omit them in the following. In the BOA the Hamiltonian could be written as follows: T fis in eq. (9) is the fission mode kinetic energy operator where M (R 12 ) is the average inertia related to the fission mode andĤ perp is the collective Hamiltonian related to the neck and mass asymmetry coordinates. The eigenfunction of Hamiltonian (9) can be written as a product, where ϕ n are the eigenfunctions ofĤ perp , and they are evaluated for each mesh-point value in the R 12 direction. Using the above relations one can rewrite the eigenequation of the Hamiltonian (9) in the following form: The average inertia M (R 12 ) in the kinetic energy operator T fis (eq. 10) is evaluated as follows: where ϕ 0 is the eigenfunction of the HamiltonianĤ perp corresponding to its lowest eigenenergy for a given value of R 12 . The integration takes place over the whole region of allowed values of a and α. The approximate solution of the above eigenvalue problem can be obtained using the WKB formalism. The energies e n (R 12 ) in eq. (13) define the fission potential for different channels, which is important when one describes the nonadiabatic fission process in the coupled channels approach [13]. In the following we shall take only the lowest energy channel, which corresponds to the adiabatic approximation. Within this approximation the wave function of the fissioning nucleus is written in the form of a product of the wave function u 0E (R 12 ) describing the motion towards fission and the function ϕ 0 (a, α 4 ; R 12 ) which corresponds to the lowest eigenenergy e 0 of the Hamiltonian (12). The probability of finding a nucleus, at a given value of R 12 , in the (a, α 4 ) point is equal to |ϕ 0 (a, α 4 ; R 12 )| 2 .
In refs. [12,13] it is shown how to go beyond the BOA and include nonadiabatic and dissipative effects, but in the following we are going to omit such effects, which are expected to be small at low temperatures, and limit our discussion to the effect of the neck degree of freedom on the fission fragment mass distribution.

Numerical results
The potential energy surface is evaluated for 236 U at zero temperature within the macroscopic-microscopic model in which the macroscopic part of the energy E LD (q) is obtained using the liquid drop formula and the microscopic shell δE shell and pairing δE pair corrections are calculated using the Woods-Saxon single-particle potential. All parameters of the potential are described in ref. [16]. The potential energy was calculated in the 4-dimensional space of deformation parameters α, α 1 , α 4 , α 6 and minimized then with respect to α 6 keeping the parameter R 12 fixed. I.e., for each fixed value of α 1 , α 4 , α 6 the dependence of the potential energy on α (α is the main elongation parameter in the Cassini parametrization) was transformed by interpolation into a dependence on R 12 . Then V (R 12 , α 1 , α 4 , α 6 ) was minimized with respect to α 6 keeping R 12 , α 1 , α 4 fixed. In this way we obtain the potential energy as a function of three deformation parameters, defines uniquely the mass asymmetry. The mass parameter M ij (q) for the fission process is commonly calculated by the Inglis formula where |0 and |m denote the ground and an excited state of the system. In the case that the ground and the excited states of the system are described by the BCS approximation, M ij (q) is given by [20], where the term P ij stands for the contribution due to the change of occupation numbers, when the deformation varies. The quantities E μ , u μ and v μ in (17) are the quasi-particle energies and coefficients of the Bogolyubov-Valatin transformation, correspondingly, and Unfortunately, the expression (17) has the very unpleasant feature that it does not turn into the mass parameter of a system of independent particles when the pairing vanishes, Δ → 0. More precisely, the nondiagonal sum over single-particle states in (17) does turn into the mass parameters of the system of independent particles when Δ → 0, but the diagonal sum goes to infinity in that limit (it is proportional 1/Δ 2 , as demonstrated in [20]). Thus, at some points in the deformation space, where the density of single-particle states is very low, the mass parameter (17) becomes unreasonably large. The same happens in excited systems, when the temperature is close to its critical value T crit at which the pairing gap disappears.
One should also keep in mind that the diagonal contribution to the sum in (17) comes from the matrix elements between the ground state and the pair excited states that correspond to a particle number different from that of the ground state. In a particle number conserving theory such contribution could not appear.
In order to avoid the problems related to the diagonal contributions to (17) we have omitted in (17) the diagonal matrix elements μ|∂H/∂q i |μ and taken into account only the nondiagonal part of (17), The suppression of dangerous diagonal terms leads to a diminishing of the inertia tensor. The relative contribution of diagonal terms to (17) depends very much on the point in the space of deformation parameters, but, in principle, it is not so small. For a particular case a = 0.1, α 4 = 0, the diagonal part of M R12R12 is by one order of magnitude smaller than the nondiagonal one, but the diagonal part of M aa is comparable with the nondiagonal ones. From the calculations it follows that the diminishing of inertia tensor increases the width of mass distribution of fission fragments. Luckily, when the diagonal part to the mass tensor is omitted the width of mass distribution gets very close to the experimental data.
The mass parameter (18) was calculated with respect to three deformation parameters, q = {R 12 , α 1 , α 4 }. We have checked that the account of α 6 does not change much the value of the mass parameter (18). So, (18) was calculated first for the coordinates α, α 1 , α 4 , α 6 = 0, then the α-mass parameter was transformed into the R 12 -mass parameter using the relation between the derivatives ∂ ∂R 12 = ∂α ∂R 12 and, eventually, the value of R 12 , α 1 , α 4 -mass parameter was interpolated from the grid point in α, α 1 to the grid point in R 12 , a like it was done for the potential energy. The potential energy surface (PES) related to the spherical liquid drop energy (upper row) and the M R12R12 component of the inertia tensor (middle row) as well as the neck radius κ = r neck /R 0 (lower row) are plotted in fig. 2 on the (A f , α 4 ) plane for three different values of the relative distances of the center of the fragments R 12 . Here A f is the atomic mass of fission fragment, A f = A(1+a)/2.
Comparing the PES maps presented in the upper row of fig. 2 one observes a transition from the minimum at A f = 140 at R 12 /R 0 = 2.1 to A f = 132 at R 12 /R 0 = 2.3. The calculations with the smaller step in R 12 show that the mass of the heavy fragment changes continuously from A f = 140 to A f = 132. The minimum at A f = 132 appears to be due to the strong shell effects in the heavier fragment with the magic numbers N = 82 for neutrons and Z = 50 for protons. At large elongation the fragment gets more spherical and the shell structure of double magic 132 Sn comes into play.
The fission fragment mass yield is obtained by solving the quantum mechanical problem of the collective Hamiltonian which describes the fission process in the three dimensional space (3D) composed of the deformation parameters R 12 , a, α 4 defined above.
In our calculation all transport parameters: potential energy surface and six components of the inertia tensor, were tabulated in 33915 grid points in the effective 3D spaces spanned by the R 12 4)) spanned by the α, α 1 , α 4 , and α 6 deformation parameters. The eigenproblem of theĤ perp Hamiltonian (12) was solved by the diagonalization method in the basis of 2D harmonic oscillator wave functions from the major shells with the phonon number smaller than or equal to N max = 16.
First one has to prepare a set of the 2D mass distribution probabilities on the plane (a, α 4 ) by solving the eigenproblem of a corresponding 2D collective Hamiltonian for fixed elongations R 12 [21]. One has to bear in mind that it is very unlikely that fission occurs at some fixed R 12 or when the system reaches the scission line/surface. The problem is much more complicated and one has to take into consideration the size of the neck. Looking at the probability distributions integrated over α 4 for 236 U presented in the top part of fig. 3 w(a, R 12 ) = |Ψ (a, α 4 ; R 12 )| 2 dα 4 , one cannot see any qualitative change with respect to the results which we have published in ref. [13] for the calculation made in the (R 12 , a) plane. Both distributions, i.e. the one corresponding to A f (1) = 140 for smaller R 12 and the one for A f (1) = 132 at R 12 close to the scission line are only slightly broader. So, the problem to reproduce the experimental distribution of fission fragments, seen in the 2D space [13], remains also in the 3D space. The only solution is to assume that fission occurs with a certain probability before (or after) reaching the critical elongation R crit 12 at which the liquid drop splits into two fragments. Depending on the neck radius a fissioning nucleus has to make its choice to fission or not to fission. When it decides for fission it would leave the phase-space of collective coordinates. Of course this is not a Hamlet dilemma, where there is only the choice between yes or no. We are rather faced here with a statistical problem and the answer yes is given with a certain probability which one then will have to take into account in the distribution probability (20) in the phase space. The geometrical criteria introduced above for the appearance of fission are in line with the idea of the macroscopic-microscopic model used here to evaluate the potential energy surface, where the charged nuclear liquid drop plays a dominant role.
Following such an assumption a part of the events (read trajectories in the Langevin approach, or distributions in our quantum mechanical model) disappears from the phase-space and leads to a kind of weighting of the mass distribution corresponding to different elongations R 12 . To do this one has to evaluate the neck ra-dius in the whole 3D space. The neck radius parameter κ = r neck /R 0 is plotted in the lower row of fig. 2 on the (A f , α 4 ) plane for three different values of the elongation parameter R 12 . The slight wiggles in fig. 2 are caused by the approximate minimization with respect to the α 6 on the mesh with finite grid size. It is seen that on average the neck radius decreases with growing R 12 . For a constant R 12 the neck radius varies strongly with α 4 and A f . One commonly agrees that fission takes place when the neck radius becomes of the order of the size of a nucleon. This is the case for κ ≈ 0.2, which is realized at R 12 /R 0 = 2.0 for α 4 = −0.18; R 12 /R 0 = 2.25 for α 4 = 0, and R 12 /R 0 > 2.5 for α 4 = 0.18 and the asymmetry parameter a ≈ 0.2. From the optimal shape approach [18] one knows that the scission shape corresponds to r neck ≈ 0.3R 0 and R crit 12 ≈ 2.3R 0 . In the case of the Cassini parametrization used in the present paper the r neck can be somewhat smaller.
One could try to parametrize the neck-rupture probability P in the following form: where k is the momentum in the direction towards fission (or simply the velocity along the elongation coordinate R 12 ), while κ = κ(a, α 4 , R 12 ) is the deformationdependent relative neck size.  plays no essential role, and will disappear from the final expression of the mass distribution when one normalizes it. The geometry-dependent part of the neck breaking probability is taken in the form of a Fermi function The parameters κ 0 and d have to be fixed by comparing the theoretical fission fragment mass distributions with the experimental ones. Our goal is to fix these parameters in a kind of universal way, independent of the specific fission reaction that one wants to investigate. The present investigation has to be treated only as a first attempt in this direction. The momentum k which appears in the denominator of eq. (22) has to ensure that the probability depends on time in which one crosses the subsequent intervals in R 12 coordinates: is the velocity towards fission. The value of k depends on the difference E −V (R 12 ) and on the part of the collective energy which is converted into heat Q, In the quantum mechanical picture the heat Q can be replaced by the imaginary part of the collective potential [13]. In the Langevin picture the method should be almost the same but one has also to work at least in the 3D space. In our present calculations we have put Q = 0, i.e., we assumed a "complete acceleration" scenario: no dissipation takes place, which is reasonable since at low excitation energies the friction force is very week.
The M in (24) is the cranking inertia relative to the R 12 deformation parameter. In principle, we used the definition of cranking inertia, but for the reasons explained above the contributions from the diagonal matrix elements were removed.
The fission probability w at a given R 12 and a will be given by the integral The dependence of the fission probability (25) on the mass asymmetry is shown in the bottom part of fig. 3 for a few values of R 12 . One observes that due to the Fermi function in (23) the contribution of larger R 12 (smaller necks) is enhanced and contributions from smaller R 12 are suppressed.
From the maps of the potential energy surface shown in the top part of fig. 2 one observes that for R 12 ≤ 2.1 the neck parameter is α 4 = −0.05 while for R 12 ≥ 2.2 the minimum at the PES is at α 4 = −0. 15. This means that a large part of the distribution probability |Ψ (a, α 4 ; R 12 )| 2 will undergo fission also at smaller elongations and one has to subtract this part from the phase-space, i.e. to diminish the initial distribution by subtracting the events which have already fissioned. The final (measured) mass distribution of the fission fragments will be the sum of those subtracted events.
Such an approach means that the fission process should be spread over some region of R 12 and that for given R 12 at fixed mass asymmetry one has to take into account the probability to fission at previous R 12 points. I.e., one has to replace w(a, R 12 ) by w (a, R 12 ) = w(a, R 12 ) (26) The effect of the replacement (26) is demonstrated in fig. 4. It is seen there that it substantially reduces the magnitude of the fission probability at large R 12 as expected.
The mass yield will be the sum of all partial yields at different R 12 , As it is seen from (27) the scaling factor k 0 in the expression for P , eq. (22), has vanished and does not appear in the definition of the mass yield. Our model will thus only have two adjustable parameters, κ 0 and d, that appear in the neck-breaking probability (23).   A comparison of the measured [22] and here calculated (eq. (27)) fission fragment mass distributions is shown in fig. 5 for the thermal neutron-induced fission of 235 U. One can see that the calculated mass distribution is very close to the experimental values. The double-peak structure, the position and the relative magnitude of the peaks are reproduced rather well.

Conclusions
The extended Cassini ovals deformation parameters and the macroscopic-microscopic model (E LD plus the Woods-Saxon single-particle potential) gives yields for 236 U a PES with an asymmetric fission valley corresponding to A f ≈ 140 when the relative distance between the fragment mass centers is smaller than R 12 = 2.3R 0 . At larger elongations the bottom of the fission valley is shifted towards A f ≈ 132. This shift has an important consequence, it allows to describe the second (inner) peak in the measured mass distribution.
We have shown that the three-dimensional quantum mechanical model which couples the fission, neck and mass asymmetry modes is able to describe the main features of the fragment mass distribution when the neck-dependent fission probability is taken into account. The obtained mass distribution is slightly shifted, by approximately 2 mass units, towards symmetric fission as compared with the experimental mass yield, but reproduces nicely the structure of the distribution observed in the experiment. This shift could be partly due to a too large stiffness of the LD energy in the mass asymmetry degrees of freedom and/or to a lack of the nonadiabatic effects (beyond the Born-Oppenheimer approximation) [13] which makes the distributions slightly wider than the sole adiabatic ones. Also the energy dissipation, not taken into account in the present investigation, could modify somewhat the theoretical distribution.