Density functional theory approach for the description of PES

For the local density approximation of DFT, the total energy of the finite nuclei can be calculated using the space integral of the Hamiltonian density $\mathcal{H}(r)$, which consists of the kinetic energy τ, potential energy χt, and pairing energy ${\stackrel{⌣}{\chi }}_t$ densities: $\mathcal{H}(r)=\frac{{\hslash }^2}{2m}\tau (r)+{\displaystyle \sum _{t=\mathrm{0,1}}{\chi }_t}(r)+{\displaystyle \sum _{t=\mathrm{0,1}}{\stackrel{⌣}{\chi }}_t}(r)\mathrm{.}$ (1) Here, $\tau (r)$ is the density of the kinetic energy, and the symbol $t=\mathrm{0,1}$ indicates the isoscalar or isovector, respectively [41].

The mean-field potential energy in the Skyrme DFT usually has the following form: $\begin{array}{c}{\chi }_t(r)=C_t^{\rho \rho }{\rho }_t^2+C_t^{\rho \tau }{\rho }_t{\tau }_t+C_t^{J^2}{\mathbb{J}}_t^2\\ +C_t^{\rho \text{Δ}\rho }{\rho }_t\text{Δ}{\rho }_t+C_t^{\rho \nabla J}{\rho }_t\nabla \cdot J_t.\end{array}$ (2) Here, the particle density ρt, kinetic density τt, and spin current vector densities $J_t(t=\mathrm{0,1})$ can be obtained from the density matrix ${\rho }_t\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)$, depending on the spatial ($r$) and spin (σ) coordinates. In the aforementioned formula, $C_t^{\rho \rho }$, $C_t^{\rho \tau }$, and etc. are the coupling constants in the Hamiltonian density $\mathcal{H}(r)$, corresponding to the different types of densities, most of which are real numbers; $C_t^{\rho \rho }=C_{t0}^{\rho \rho }+C_{tD}^{\rho \rho }{\rho }_0^{\gamma }$ is an exception, which is the traditional density-dependence term. Expressions relating the coupling constants to the standard Skyrme parameters can be found in Ref. [42]. Specifically, the spin-orbit interaction in the Skyrme force corresponds to the term $C_t^{\rho \nabla J}{\rho }_t\nabla \cdot J_t$.

In the DFT, the pairing correlation is usually incorporated by the Hartree-Fock-Bogoliubov (HFB) method [37]. For the Skyrme energy density functional, a commonly used pairing force is the density-dependent, zero range potential, which can be expressed as follows: Refs. [8, 43]: ${\widehat{V}}_{\text{pair}}\left(r\mathrm{,}{r}^{\text{'}}\right)=V_0^{\left(n\mathrm{,}p\right)}\left[1-\eta {\left(\frac{\rho (r)}{{\rho }_0}\right)}^{\gamma }\right]\delta \left(r-r^{⊵}\right)\mathrm{,}$ (3) Here, $V_0^{\left(n\mathrm{,}p\right)}$ is the pairing strength for neutrons (n) and protons (p), the exponent γ of the density-dependence affects the appearance of the neutron skins and halos [44], for which γ =1 is widely used [45, 46]. ρ₀ is the average density inside the nucleus (often considered as the saturation density of nuclear matter, and set as 0.16 fm^-3), and $\rho (r)$ is the total density. Different types of pairing forces can be obtained by choosing different values of η. The pairing force is ”volume-like”, which indicates that if η = 0, there is no explicit density dependence. The pairing force acts equivalently inside the nuclear volume. In contrast, the pairing force will be “surface-like” for η = 1, which has a significant effect around the nuclear surface and a small impact in the center area of the nucleus. The choice of $\eta =\frac{1}{2}$ is often called a mixed-pairing force, which is the average of these two types of pairing forces. To test the sensitivities on the fission-related properties based on the form of the pairing force, we considered η = 0, 0.25, 0.5, 0.75, and 1. η = 0.25 and 0.75 were used for test purposes, and the other choices of the η values have been frequently used for structural studies.

Time-dependent generator coordinate method for fission dynamics

Fission is a large-amplitude collective motion and can be described as a slow adiabatic process driven by a few collective degrees of freedom within the framework of TDGCM. In this approach, the many-body state wave function of the fissioning system can be expressed by the following generic form: $|\Psi (t)\rangle ={\displaystyle {\int }_q}f(q\mathrm{,}t)|\text{Φ}(q)\rangle \text{d}q\mathrm{.}$ (4) where $|\text{Φ}(q)\rangle$ consists of a set of known many-body wave functions parameterized by a vector of continuous variables $q$. Each of these $q$ is a collective variable and is chosen based on the specific physics problem. The quadrupole moment ${\widehat{Q}}_{20}$ and octupole moment ${\widehat{Q}}_{30}$ are usually chosen as the collective variables for the fission study. $f(q\mathrm{,}t)$ is the weighted function, which is solved by using the time-dependent Schrö dinger-like equation in the space of the coordinates q. Under the GOA, this equation can be expressed as follows: $i\hslash \frac{\partial g(q\mathrm{,}t)}{\partial t}={\widehat{H}}_{\text{coll}}\left(q\right)g(q\mathrm{,}t)\mathrm{.}$ (5) The collective Hamiltonian ${\widehat{H}}_{\text{coll}}\left(q\right)$ is as follows: ${\widehat{H}}_{\text{coll}}\left(q\right)=-\frac{{\hslash }^2}{2}{\displaystyle \sum _{ij}\frac{\partial }{\partial q_i}}{B}_{ij}\left(q\right)\frac{\partial }{\partial q_j}+V\left(q\right)\mathrm{,}$ (6) where V(q) is the collective potential, and the inertia tensor $B_{ij}\left(q\right)={\mathcal{M}}^{-1}\left(q\right)$ is the inverse of the mass tensor $\mathcal{M}$. The potential and mass tensor are determined by the Skyrme DFT in this study. g(q, t) is the complex collective wave function of the collective variables q and contains all the information regarding the dynamic of the fission system.

For the description of fission, the collective space is divided into an inner region with a single nuclear density distribution, and an external region that contains the two fission fragments. The scission contour defines the hyper-surface that separates the two regions. The flux of the probability current through this hyper-surface provides a measure of the probability of observing a given pair of fragments at time $t$. The integrated flux $F(\xi \mathrm{,}t)$ for the surface element ξ on the scission hyper-surface is calculated as follows: $F(\xi \mathrm{,}t)={\displaystyle {\int }_{t=0}^t}\text{d}t{\displaystyle {\int }_{qϵ\xi }}J(q\mathrm{,}t)\cdot \text{dS}\mathrm{,}$ (7) as in Ref. [40], where J(q, t) is the current $J(q\mathrm{,}t)=\frac{\hslash }{2i}B(q)\left[g^{*}(q\mathrm{,}t)\nabla g(q\mathrm{,}t)-g(q\mathrm{,}t)\nabla g^{*}(q\mathrm{,}t)\right]\mathrm{.}$ (8) The yield of the fission fragment with mass number A can be calculated as follows: $Y(A)=C{\displaystyle \sum _{\xi ϵ\mathcal{A}}\underset{t\to +\propto }{{\displaystyle \lim }}}F(\xi \mathrm{,}t),$ (9) where $\mathcal{A}$ indicates a set of all the surface elements ξ on the scission hyper-surface with a fragment mass of A, where C is the normalization constant to ensure that the total yield is normalized to 200 as usual. The yield of the fission fragment with the charge number Z can also be obtained from the integrated flux $F(\xi \mathrm{,}t)$. The mass number A can be replaced with the charge number Z in Eq. 9; the summation is then over the set of all the ξ values on the scission frontier with the fragment charge number Z. In this study, the FELIX(version 2.0) [47] computer code was used for modeling the time evolution of the fissioning nucleus in the framework of TDGCM+GOA.

Potential energy surface

As indicated in Sec.2.B, considering the adiabatic approximation approach for fission dynamics, obtaining the precise multidimensional PESs is the first step toward the dynamical description of fission. In this study, we chose the quadruple moment (q₂₀) and octuple moment (q₃₀) as the collective parameters, which are the most important collective degrees of freedom for the nuclear fission study; they describe the elongation of the nucleus and mass asymmetry, respectively. Figure 2 presents the PESs of ²⁴⁰Pu calculated by using the HFB method with five different types of pairing forces (η = 0, 0.25, 0.5, 0.75, and 1.0) in the collective space of (q₂₀, q₃₀). The collective variables ranged from 0 to 600 b for q₂₀, and from 0 to 60 b^3/2 for q₃₀ with the step of $\text{Δ}{q}_{20}=2$ b and $\text{Δ}{q}_{30}$=2 b^3/2.

Fig. 2Full-screen View

(color online) The PESs of ²⁴⁰Pu in the (q₂₀, q₂₀) plane calculated with the SkM*-DFT and the five different types of pairing forces in the HFB approximation. The value of η for the different types of pairing forces is provided in each panel. The energies are provided relative to the ground-state energy. The contours join the points on the surface with the same energy, and the interval of the contour is 1.0 MeV. The least-energy fission pathway is indicated by a solid red curve

As shown in Fig. 2, there is no notable difference in the topological properties of PESs with different types of pairing forces. Double-humped fission barriers were predicted for all the cases. An inner symmetric fission barrier followed by an outer asymmetric barrier was clearly distinguished. At q₂₀>200 b, symmetric valleys with large elongations were found. The symmetric and asymmetric fission valleys were well-separated by a ridge from (q₂₀, q₃₀) ≈ (150 b, 0 b^3/2) to (350 b, 20 b^3/2), and the height of the ridge gradually decreased as the η value increased. Therefore, the density-dependent surface pairing force led to the reduction of the ridge height. In addition, the asymmetric fission channel was favored for all the least-energy fission pathways, as indicated by the red lines in Fig. 2.

Energies of the symmetric and asymmetric fission paths as a function of the quadruple moment (q₂₀) are provided in Figs. 3(a) and (b), respectively. The value of η varied from 0 to 1, indicating that the ”volume-like” pairing force transitioned into a surface pairing force. Fig. 3 demonstrates that the fission barrier heights and isomeric-state energy decrease as η increases. Specifically, when η>0.5, the fission barriers explicitly decrease. For the least-energy fission pathway shown in Fig. 3(b), a smaller quadruple moment is needed for the occurrence of the scission for a larger η.

Fig. 3Full-screen View

(color online) Energies along the symmetric and least-energy fission pathways in ²⁴⁰Pu, with different types of pairing correlations, are shown in panels (a) and (b), respectively. All these energies are relative to their ground-state values

Table 1 lists the energies of the ground state, isomeric states, and fission barrier heights for the different types of pairing forces, along with the corresponding quadruple and octuple moments for each state. The energies of the isomeric state and heights of the fission barrier decrease as the η value increases. Owing to the lack of triaxial deformation in the collective space in our calculation, the heights of the fission barriers would be higher than those experimentally obtained, especially for the inner fission barrier [52, 53]. As indicated in the table, for η=1, the inner fission barrier height was near that demonstrated by the data, and the outer fission barrier was lower than that in the data, leaving no room for the triaxial degree of freedom. Thus, η=1, that is, the surface pairing force, may not be a good choice for the fission study. Based on Table 1, the deformations of these states, including the ground state, isomeric state, and inner and outer barriers, are generally not influenced by the type of pairing force. These deformations are mainly determined by the shell structure given by the mean-field potential. In our previous study [22], we also found that these deformations were relatively stable against the variations of the pairing strength.

Pairing energies

Figure 4 presents the pairing energies at different deformations for various types of pairing forces. The pairing energies at the ground and isomeric states are smaller than those at the fission barriers. At the same state, the pairing energy increases as the value of η increases, especially for η=1, which is the surface-type pairing force. The pairing gaps at different deformations are provided in Fig. 5. Once again, the pairing gap has a minimum at the ground state and a second minimum at the isomeric states. The pairing gaps are large around the fission barriers. For the adjustment of the strength of the different pairing forces used in this study, the same values of the pairing gaps in ²⁴⁰Pu were used. The figure demonstrates that at the smaller deformation region, the pairing gaps from the different types of pairing forces were relatively similar. However, when the deformation was large, explicit discrepancies appeared (q₂₀>150 b). For the pairing force with a smaller η value, the neutron and proton pairing gaps were generally smaller.

Fig. 4Full-screen View

(color online) Pairing energy of the ground state, fission isomer state, and inner and outer fission barriers in ²⁴⁰Pu as a function of the η parameter in the pairing force

Fig. 5Full-screen View

(color online) Pairing gap along the least-energy fission path for the neutron (a) and proton (b) with different types of pairing forces

Mass tensor

The mass tensor $\mathcal{M}$ reflects the response of the fissioning system to the collective coordinate changes. In this study, mass tensors were obtained from the static calculations by using the Skyrme DFT with the perturbative cranking approximation. As shown in Fig. 6, the elements of the mass tensor M₂₂, M₃₃, and M₂₃ were plotted as functions of q₂₀ along the lowest-energy fission path for the different choices of η. The mass tensor given by ATDHFB was larger than that by GCM for all types of pairing forces. As indicated in Ref. [8], this was caused by the missing correlations in the GCM method. Generally, as the quadruple moments increase, M₂₂ and M₃₃ gradually decrease. M₂₃ is negative, and its absolute value increases and has explicit fluctuations when the deformation is large. For different types of pairing forces tuned by η, the mass tensor apparently decreases and demonstrates reduced fluctuation against deformations when η is large. As indicated in Ref. [28], considering the mixed pairing force of η=0.5, the mass tensor decreases and fluctuates less when the pairing strength is decreased. The systematic behaviour of the η>0.5 pairing force resembles that of increasing the strength of the pairing force.

Fig. 6Full-screen View

(Color online)M₂₂ (top panel), M₂₃ (middle panel), and M₃₃ (bottom panel) of the mass tensors are shown as the elongated deformations along the static fission path for different η values. The resulst in the left column are obtained by using GCM method, and those in the right column are derived by the ATDHFB method

Scission lines

Determining the scission frontier is critical for describing the fission dynamics. In DFT, the operator $q_{\text{N}}=\langle {\widehat{Q}}_{\text{N}}\rangle =\langle e^{-{\left(\frac{z-z_{\text{N}}}{a_{\text{N}}}\right)}^2}\rangle$ is often used to evaluate the neck size of the fissioning nuclei. a_N= 1 fm is often chosen. z_N is the neck position, which has the lowest density between the two fragments. Generally, the neck size smoothly decreases as the fissioning nucleus elongates, and decreases to nearly zero after the scission, where the two fragments are sufficiently separated. In this study, we chose q_N= 4 as the critical value for determining the scission line in ²⁴⁰Pu, which has been used for ²⁴⁰Pu in Refs. [40, 28]. This value was chosen at the edge of the sudden decrease in the neck size, maintaining most of the pre-fission configurations for further fission dynamic calculations. The scission lines in the PES of the (q₂₀, q₃₀) collective spaces obtained by DFT using the different types of pairing forces are shown in Fig. 7. Generally, the scission contours of the different η values display similar patterns. For different η values, in or around the region of the symmetric fission, the increase of η leads to a smaller quadrupole moment at the scission point. Considering η= 1.0, symmetric fission occurred at the quadrupole moment q₂₀ ~ 480 b for η =1.0, and at approximately 550 b for the other η values. The shortest elongation occurred at q₂₀ ~ 300 b when the asymmetry increased to the octupole moment q₃₀ ~ 30 b^3/2. Subsequently, the scission lines turned toward the upper-right direction until significant asymmetry appeared. The pre-fission region for η=1 is explicitly smaller than the other cases.

Fig. 7Full-screen View

(color online) The scission lines using the criterion of q_N= 4 are shown for the different η values

Total kinetic energy

An important quantity in induced fission is the total kinetic energy (TKE), which is obtained by the fission fragments. In this study, the total kinetic energy of the two separated fragments at the scission point can be estimated as the Coulomb repulsive interaction $E_{\text{TKE}}=\frac{e^2{Z}_H{Z}_L}{d_{\text{ch}}}$, where e indicates the proton charge, Z_H and Z_L denote the charge number of the heavy and light fragments, respectively, and d_ch is the relative distance between the centers of charge of the two fragments at the scission point. This approximation for TKE has been frequently used for simplification, as demonstrated in Refs. [15, 20, 22, 57-59]. However, it neglects the dissipation and shell effects, among others, which can lead to an overestimated TKE compared to the experimental data [15, 20, 22, 57-59]. The dissipation effect has been recently considered [60], thus the calculated TKE can better agree with the data. The TKE values of the ²⁴⁰Pu fission fragments with different types of pairing forces were plotted as functions of the heavy fragment mass, as shown in Fig. 8. The open circles represent the calculated results for different η values, whereas the solid circles indicate the experimental data from the thermal neutron-induced ²³⁹Pu fission experiments [61-63]. Considering the general trend, a qualitative dip is reproduced for all types of pairing forces at A_H = 120, as well as a peak at A_H = 134.Near the dip or peak, the TKE is larger for larger η values. For A_H > 144, the discrepancies are fairly small for different pairing forces.

Fig. 8Full-screen View

(color online) Total kinetic energies of the nascent fission fragments as functions of the heavy fragment mass for ²⁴⁰Pu based on the different η values. Open symbols represent the calculated results, and solid symbols indicate the experimental data [61-63] for the thermal neutron

Fission yields

Figure 9 presents the mass and charge yields obtained with the code FELIX (version 2) [47] based on the TDGCM+GOA framework using the different types of pairing forces, which are compared with the experimental data. As a critical microscopic input of fission dynamic calculations, the mass tensor is calculated by the GCM or ATDHFB methods. In this calculation, qN= 4 was used to determine the scission line. Generally, the discrepancies between the calculated pre-neutron mass distribution and charge distributions obtained by using the mass tensors by the GCM and ATDHFB methods are small. Furthermore, the mass and charge yields calculated by using the mixed-pairing force with η =0.5 combined with the mass tensor by the ATDHFB method demonstrated the best agreement with the experimental data.

Fig. 9Full-screen View

(color online) Calculated pre-neutron mass yields [panels (a) and (b)] and charge yields [panels (c) and (d)] are compared with the experimental data. Only the heavy fragments are displayed. In panels (a) and (c), the mass tensor from the GCM method is used, and in panels (b) and (d), the mass tensor is calculated using the ATDHFB method. Data for the pre-neutron mass yields were obtained from Refs.[64, 65], and those of the charge yield were obtained from Ref.[66]

The impact of the different types of pairing forces on the mass and charge distributions is apparent. For the calculated results obtained by the ATDHFB mass tensor, the position of the peak was nearly constant for η= 0.0, 0.25, and 0.5, and moved toward the heavy fragment as η= 0.75 and 1.0. For the results with GCM mass tensor, the mass and charge distributions of the fission fragment shifted toward the more heavy fragment as η increased (panels (a) and (c)). Furthermore, the theoretical calculations obtained by the TDGCM with ATDHFB mass tensors (panels (b) and (d)) demonstrated that the yields from the symmetric fission channel increased as η increased, which was related to the decrease in the height of the ridge as η increased, as shown in Fig. 2.