Microscopic Theory for Spontaneous Fission

Nuclear fission is a fascinating field of research that involves large-amplitude collective dynamics of a microscopic many-body system. Specifically, the process of spontaneous-fission decay can only be explained within the quantum tunneling phenomenon. The present review discusses recent advancements in the theoretical understanding of spontaneous fission. These concern precise prediction of the spontaneous fission observables like fission lifetime and distribution of fragment yields. The theoretical developments presented here are based on a coherent coupling between the adiabatic collective dynamics and the static inputs obtained from the nuclear energy density functional formalism.


INTRODUCTION
Nuclear spontaneous fission (SF) is a unique decay mechanism that has crucial applications in both basic and applied sciences [1][2][3]. Particularly, the stability of very heavy and superheavy nuclei strongly depends on the SF probability [4][5][6]. Although superheavy nuclei predominantly decay via α-emission at the beginning of a decay chain, SF leads to terminate the chain. This type of decay sequences are experimentally observed for isotopes of Fl and Ts [7,8]. Moreover, in comparison to α-emission, SF is predicted to be the preferred decay mode for neutron-rich superheavy nuclei [5,9,10]. In the case of nuclear astrophysics, SF strongly impacts the abundances of heavy elements in stars by participating in the r-process recycling mechanism [11][12][13]. Specifically, distributions of fission-fragment yields from different fission modes (SF, beta-delayed fission, and neutroninduced fission) are essential components of the r-process abundances [9,[14][15][16][17][18][19] and, therefore, very accurate prediction of these yields is required to improve the r-process network calculations. Further, as suggested in a recent study [20], the precise estimation of fission yields is indispensable for a better understanding of the chemical evolution of r-process elements produced in binary neutron-star mergers. In the application frontier, SF data are important to calibrate the nuclear material counting techniques relevant to power generation and international safeguards [21,22]. However, measurements in actinide nuclei are restricted due to safety issues. Therefore, for both basic science and applications, predictive modeling of SF observables is of utmost interest. The present scenario and the prospects of fission theory are described in a recent review [23].
In the SF process, a fissioning nucleus undergoes quantum tunneling through a single or multiple potential barriers generated by the coherent motion of strongly interacting nucleons. Ideally, the time-dependent density functional theory (TDDFT) provides the most realistic microscopic framework to deal with such large amplitude collective dynamics [24][25][26][27]. Specifically, in the characterization of fission yields, nuclear dissipation plays a crucial role near the scission configuration and would be best accounted for by TDDFT. Albeit very promising, TDDFT poses several limitations in its application to SF. Primarily, the quantum tunneling is energetically forbidden within this semiclassical approach. Besides, current implementations of TDDFT simulate only a single fission path for a given initial condition; reconstruction of a full yield distribution requires large-scale Monte Carlo sampling, which is beyond the current computational capabilities. The exascale computing platforms may open up the avenues to overcome such restrictions [23].
In general, the collective dynamics of a nucleus is believed to be a much slower process than the random motion of the constituent nucleons. Consequently, majority of the fission models are implemented within the adiabatic approximation that segregates the collective degrees of freedom from the intrinsic coordinates. The same approach is adopted in a static nuclear density functional theory (DFT) based model [28], where the collective motion is simulated by incorporating the DFT inputs within an appropriate equation of motion. Specifically, in case of SF, the dynamics of a fissioning system can be divided into two successive steps as depicted in Figure 1 [29]. In the first step, the system tunnels through a multidimensional space of collective coordinates. This process is mainly governed by the potential energy profile and the collective inertia, which is often calculated within the adiabatic time-dependent Hartree-Fock-Bogoliubov (ATDHFB) formalism [30,31]. The region beyond the outer turning point ("out" in Figure 1) is energetically accessible and the time-evolution in this space can be followed with a simple classical prescription, e.g., the Langevin dynamical model. Finally, the system breaks into two fragments at the scission configuration. The dynamics in the second phase involves the collective inertia and the dissipative forces as well. The SF halflife is primarily decided by the tunneling phase as, for most of the relevant nuclei, it is significantly slower than the subsequent propagation outside the barrier. The second dynamical phase controls the fission fragment properties like yield and total kinetic energy (TKE) distributions. FIGURE 1 | Variation in potential energy calculated along the most-probable fission path of 240 Pu. The region marked as WKB is classically forbidden as the fission excitation energy E 0 is less than the potential energy. Shape evolution: the collective action is minimized between the turning points "in" and "out" by using the WKB method and, from "out" to scission, dynamics is governed by Langevin dynamics. The figure is adapted from Sadhukhan et al. [29].
Apart from the standard inputs discussed in the previous paragraph, paring correlations play a critical role in controlling the SF lifetime and the connected fission pathway. The individual nucleonic motion leads to shell structures that guide both the nuclear shape and the potential energy along a fission path. Moreover, crossings of single-particle levels can modify the collective inertia through associated changes in the single-particle configurations [32][33][34]. The residual interaction among these crossing configurations is strongly influenced by the pairing force. Precisely, a larger pairing gap helps the collective motion to be more adiabatic [35][36][37][38][39][40]. The enhancement of pairing fluctuations along the fission path was first postulated in Moretto and Babinet [41] by using a simple parabolic potential. In fact, the collective inertia and potential energy show opposite trends as changes. The potential energy increases as deviates from the static value s , obtained from the self-consistent energy minimization procedure. In contrast, the collective inertia varies as −2 [35,[42][43][44][45] and, therefore, the dynamic , corresponding to the minimum collective action, may differ from s . This suggests that the parameter should be implemented as an independent dynamical variable to determine the least action trajectory. Indeed, a reduction in the collective action due to pairing fluctuations is observed in many macroscopicmicroscopic studies [46][47][48][49]. In addition, the pairing fluctuations are recently treated as dynamical variables in microscopic models based on the DFT formalism [50,51].

Fission Half-Life
The SF mechanism involves a very wide range of timescales depending on the choice of the fissioning nucleus. For example, the observed SF half-life for actinides varies from a few ms to 10 20 yrs. Therefore, it is impractical to develop a SF model based on the real-time quantum dynamics. The most common approach for the calculation of SF half-life is rooted in the one-dimensional WKB approximation to the quantum tunneling process. The corresponding half-life can be expressed as [52,53], T 1/2 = ln 2/(nP), where n is the rate of collision on the fission barrier and P is the barrier penetration probability given by In the above equation, S(L) is the action integral calculated along a predefined fission path L(s) in the multidimensional collective space. The expression for S(L) is given by, where V(s) and corresponds to the most probable fission path [35,54]. The M eff can be expressed in terms of the multidimensional collective inertia tensor M ij (q 1 , q 2 , ...) as [52,53,55]: Generalization of the WKB method to several dimensions is recently recommended as a future goal [23]. The potential energy V is obtained by subtracting the vibrational zero-point energy E ZPE [56] from the Hartree-Fock-Bogoliubov (HFB) energy E HFB (= Ĥ HFB ). In the DFT formalism, E HFB can be computed self-consistently by solving the constrained HFB equations for the Routhian: whereĤ HFB ,Q ij , andN τ represent the HFB hamiltonian, the mass multipole moment operators, and particle-number operators for neutrons (τ = n) and protons (τ = p), respectively. The Lagrange multipliers λ 2τ can be used to control the particlenumber fluctuation terms: Pairing correlations in nucleons are interconnected to the particlenumber fluctuations [57,58] and, hence, expected to be stronger for λ 2τ > 0, compared to its static value obtained with λ 2τ = 0. Further, the overall magnitude of pairing correlations is linked to the average pairing gap [31,59]. Therefore, λ 2n and λ 2p can be utilized as dynamical coordinates [50] to scan the configuration space over a wide range of . It is indeed more physical to select the isoscalar (λ 2n + λ 2p ) and the isovector (λ 2n − λ 2p ) variations as independent coordinates. Consequently, a one-dimensional path L(s) can be identified with the collective variables {q i } ≡ {Q 20 , Q 22 , Q 30 , ..., λ 2n + λ 2p , λ 2n − λ 2p } as functions of path's length s. In a recent study [50], the role of dynamical isovector pairing is found to be negligible and, therefore, the associated coordinate λ 2n − λ 2p can be set to zero. I denote λ 2n + λ 2p as λ 2 in the subsequent discussions. Also, I should mention that an appropriate normalization scheme for all the coordinates must be adopted to make ds dimensionless in Equation (2) [50].
Minimum-action paths can be simulated by using two different techniques named as the dynamic-programming method (DPM) [52,56] and the Ritz method (RM) [53,56]. In DPM, the dynamical space is first discretized into a multidimensional mesh. Then, at any intermediate step, minimum action paths are calculated for all the mesh points on a hypersurface perpendicular to the elongation coordinate (Q 20 ). Calculation is propagated along the Q 20 direction and finally the hypersurface of the desired outer-turning point is reached. In this method, path lengths are further divided into smaller segments whenever the distance between two adjacent mesh points is large. This is essential to ensure numerical accuracy as the collective inertia may vary quite sharply in certain regions of the collective space. In the case of RM, trial paths are defined in terms of Fourier series of dynamical coordinates and the coefficients of different Fourier components are obtained by minimizing the action. Although RM is easier to implement numerically, efficiency of this method depends on the number of Fourier components needed to reproduce the actual fission path. On the other hand, DPM is free of such limitations.
The HFB energy E HFB can be calculated by employing either covariant or non-relativistic energy density functionals (EDFs). In the case of non-relativistic EDFs, the SkM* parametrization [60] of the zero range Skyrme functionals is commonly used in fission studies together with the density-dependent mixed pairing interaction [61], where the pairing strengths are calculated locally by reproducing the odd-even mass differences around 240 Pu [62]. The parameters of this parametrization are benchmarked for large deformations relevant to fission. Recent optimizations of the Skyrme functionals are performed within the UNEDF project [63] and one of its variants, UNEDF1 HFB , Schunck et al. [64] is successfully applied to fission works [65][66][67]. Apart from these, other microscopic interactions like the finite-range Gogny interaction [68][69][70] and the Barcelona-Catania-Paris-Madrid EDFs [9,69,71,72] are widely used in the SF calculations. Despite different groups of EDFs exist, appropriate benchmarking with the experimental data are performed to ensure the consistency of model predictions. For example, it is recently shown that the SF yields of the superheavy 294 Og nucleus are robust against different choices of EDFs [66]. Covariant EDFs [51,[73][74][75][76][77][78] based frameworks for large-amplitude collective dynamics are becoming more accessible with the increasing computational resources. Interestingly, SF pathways, calculated within covariant EDFs [51,74], are found to be in close agreement with the non-relativistic results.

Calculation of Fission-Fragment Yields
It is desirable to use the TDDFT framework to study the evolution of a fissioning system in the collective space beyond the outer-turning point (Figure 1: region in-between "out" and "scission"). However, as pointed out in a recent study [79], the dynamics near scission is strongly dissipative and existing versions of TDDFT are not adequate since they are lacking fluctuations in collective coordinates. Also, an advanced TDDFT framework is recently proposed that incorporates fluctuations to generate TKE and yield distributions of fission fragments. In this approach, density fluctuations are assumed to prepare an ensemble of different configurations outside the barrier region and the subsequent propagation is followed with the standard TDDFT [80]. However, the fluctuations are considered in a somewhat restricted configuration space and a more exhaustive calculation may require huge computations. A feasible alternative to TDDFT is the time-dependent generator coordinate method (TDGCM) [77,[81][82][83][84]. However, in this approach, the Gaussian overlap approximation [81,85] is additionally assumed to derive simple expressions for the parameters of the collective Hamiltonian. As a result, structural details like large fluctuations in the collective inertia are diluted [86]. It is also pointed out in the recent proposal [23] by Bender et al. that a stochastic mean-field approach with large fluctuations is more suitable for calculating fragment yields. Furthermore, the requirement of strongly dissipative dynamics for yield distributions of excited nuclei is well-established [87].
The stochastic Langevin dynamical model is a plausible option to avoid all the above-mentioned difficulties. It is quite straightforward to implement this model even in a complicated multidimensional collective space. In Langevin dynamics, the intrinsic motion of the nucleons is assumed to form a heat bath. The collective coordinates interact with the heat bath through random and dissipative forces. This decoupling of the collective coordinates from the internal degrees of freedom is performed under the adiabatic approximation. Fluctuations (random forces) introduce stochasticity in the collective dynamics and dissipation hinders the motion by transferring collective energy to the heat bath. Also, the collective motion experiences the standard conservative force exerted by potential energy. First, a family of SF probabilities P(s out ) is obtained on the hypersurface of outer turning points s out . The hypersurface should contain mass octupole moment Q 30 as one of the coordinates since this variable defines different realizations of the fragment mass and charge. Subsequently, fission paths are computed for all the s out s by solving the Langevin equations [88,89]: where p i is the momentum conjugate to q i . η ij and g ij represent the dissipation tensor and the strength of the random force, respectively, and these two quantities are connected through Einstein's fluctuation-dissipation theorem: k g ik g jk = η ij k B T.
Here, T is the temperature of the nucleus. It is calculated at each instant of the Langevin evolution by assuming the nucleus as a non-interacting Fermi gas and the resulting formula is T = √ E * /a (T in MeV); a being the level density parameter and In the case of SF studies, a can be approximated as a shape independent constant given by a = A/10 MeV −1 [29,66]. The stochastic variable Ŵ j (t) signifies the Markovian nature of the random force with the time correlation property: Ŵ k (t)Ŵ l (t ′ ) = 2δ kl δ(t − t ′ ). The excitation energy E * increase as the system slides down to lower potential resulting stronger effects from fluctuations. The scission configuration is defined with the condition that the number of neck-particles (N q ) in the fissioning system is less than a critical value [28]. Each point on the scission hypersurface uniquely identifies the particle numbers of two fission fragments and these numbers can be calculated by integrating the nucleonic density distributions [62]. Owing to the random force, an ensemble of Langevin events with the same initial configuration (i.e., same s out ) yields different fission pathways. Finally, the charge and mass distributions of the yields can be extracted by counting the number of events terminating at a given fragmentation. The numbers are weighted with P(s out ) to account for the tunneling phase. Further, to incorporate the uncertainties in N q , Langevin yields are convoluted with Gaussian functions [90].
Although the Langevin model does not explicitly simulate the time evolution of nucleonic degrees of freedom, it incorporates all the essential microscopic effects through the input quantities like PES and collective inertia which are obtained from effective nucleon-nucleon interactions. A special characteristic of the Langevin formalism is the presence of fluctuating and dissipative forces. In the case of induced fission, the importance of fluctuations and dissipation is well-established [87,91] as the compound system is produced at least with a reasonable amount of excitation energy. In SF, the ground-state zero-point vibration is the initial source of excitation energy and it is inadequate to trigger any noticeable randomness in the collective motion. However, for all the relevant nuclei, potential energy drops rapidly below its ground state value as the deformation grows beyond the tunneling region. Consequently, nuclei acquire sufficient excitation energy that enhances fluctuations to a considerable level. In the following section, I have demonstrated the impact of the random force in Equation (5) on SF pathways.

Nucleonic Localization Function and Pre-fragments
For a better understanding of the structural evolution in a fissioning nucleus, nucleonic localization functions (NLF) are calculated [25,92,93]. NLF measures the probability of finding two nucleons with the same spin σ and isospin q at the same spatial localization. It is computed as described in references [93,94]: where ρ qσ , τ qσ , j qσ and ∇ρ qσ are the particle density, kinetic energy density, current density, and density gradient, respectively. The Thomas-Fermi kinetic energy density τ TF qσ = 3 5 (6π 2 ) 2/3 ρ

5/3
qσ is introduced as a normalization parameter. A value of C ∼ 1 indicates a large nucleon's localization, i.e., a low probability of finding two nucleons with the same quantum numbers at the same spatial location. On the other hand, C = 1/2 corresponds to the limit of a homogeneous Fermi gas. The concept of NLF was originally applied to characterize chemical bonds in electronic systems. In nuclear physics, it is first used to visualize the cluster structures in light nuclei [94]. As illustrated in references [93,94], the clustering of nucleons inside a nucleus can be predicted more precisely by NLFs in comparison to the scalar density distributions given by ρ qσ . This is because the spatial distributions of NLFs exhibit pattern of concentric rings that reflect the underlying shell structure, but such patterns are averaged out in the density distributions. Recent studies [66,67,92,93] suggest that NLFs can be utilized to identify prefragments in a fissioning nucleus. The method is described in Sadhukhan et al. [92] for a typical case of elongated 240 Pu. The corresponding NLFs are shown in Figure 2. Evidently, the parts of NLFs for z ≥ z L and z ≤ z H contain ring-like patterns delineating the presence of localized nucleons inside the deformed 240 Pu. This finding is further extended to define prefragments by integrating the densities of the localized portions and doubling the result to account for reflection symmetry. As replicated in Figure 2, the compound nuclear NLFs were found to be in remarkable agreement with the ground-state NLFs of the predicted prefragments ( 128 Sn and 80 Ge) [92]. Similarly, prefragments are predicted successfully by comparing the scalar densities [95,96]. Here, I should mention that the notion of prefragments is a purely theoretical concept as it cannot be measured experimentally. Moreover, other definitions of pre-fragment exist [25,97]. Therefore, the validity of a prefragment based description depends on its ability to reproduce the experimental observables.

RESULTS ON DIFFERENT ASPECTS OF SPONTANEOUS FISSION
Calculation of SF observables is a very active field of research as the theoretical capability is increasing. In parallel, appropriate sets of fission observables are required to benchmark theoretical models [98]. In the rest of this review, I will discuss selected results from the recent theoretical achievements pertinent to both of these aspects.

Effect of Collective Inertia in Fission Pathway
In recent studies [56,74], it is demonstrated that the SF pathways are strongly influenced by the characteristics of the collective inertia. The microscopic collective inertia is usually calculated within the ATDHFB formalism [30,31] and it is commonly known as the cranking inertia M C [86]. An exact calculation of M C requires derivatives of the particle and pairing densities with respect to the dynamical coordinates. These can be achieved by employing the three-point or a higher-order Lagrange formula [99,100]. On the other hand, in case of the commonly used perturbative cranking inertia M C p [86], these derivatives are reduced to matrix elements of mass multipole moments (Q ij ). Figure 3 demonstrates the variations of squareroot-determinants of both M C and M C p calculated for 264 Fm in a two-dimensional collective space of (Q 20 , Q 22 ). As discussed in 1, M C 1/2 shows large fluctuations as an outcome of crossings in single-particle levels at the Fermi energy [39]. To affirm this, single-particle energies for 264 Fm are displayed in Figure 4 along two straight lines defined by Q 22 = 0 and Q 20 = 61 b. Multiple level crossings near the Fermi energy are clearly visible at deformations where M C changes sharply. Similar features of the inertia tensor are observed within the covariant EDF formalism [74].
The dynamical minimum-action paths (or equivalently most probable paths), obtained with M C and M C p , are drawn in Figure 5. The same figure also shows the static path that traverses the minimized collective potential [56]. Due to strong dynamical hindrance by the perturbative inertia the corresponding minimum-action path avoids large triaxial shapes. M C p varies rather smoothly along both the deformation coordinates and, therefore, the minimum in the action integral in Equation (2) is achieved by minimizing the path-length. This weaker dependency on the triaxial shapes, imposed by collective inertia, is also observed in older fission studies [101][102][103][104][105]. On the other hand, due to localized large variations in M C , the non-perturbative path passes through the triaxial shapes that are fairly close to the static pathway. Apparently, both the non-perturbative and static trajectories always adhere to a configuration that tries to minimize the density of singleparticle levels on the Fermi energy by avoiding level crossings. In short, the underestimation of structural details in M C p results in an artificial restoration of axial symmetry, which is broken in both static and non-perturbative approaches. This conclusion is also verified within the relativistic mean-field formalism [74]. Although the inertia strongly influences the topology of the minimum action path near the first fission barrier in actinides, it is rather simple in the space outside the fission isomer. Here, the SF path usually follows the minimum distance from a mass-symmetric configuration to the nearest outer turning point [29] (shown in the following figure, Figure 8).  Most importantly, apart from modifying the fission pathways, collective inertia strongly impacts the fission lifetime. The SF half-life changes by orders of magnitude depending on the choice of collective inertia, even for the same fission trajectory [56]. For example, values of T 1/2 and S(L) corresponding to different selections of the fission path and inertia are given in Table 1.

Role of Pairing Correlations
In the previous subsection, I demonstrated that a fissioning system tries to always follow single-particle configurations with comparatively lower level density. This can be fulfilled by avoiding the regions of level-crossing. In contrast, pairing correlations increases with the single-particle level density and, as I have discussed in 1, it affects the potential and collective inertia  in the opposite way. The potential energy use to increase with pairing fluctuations, while the collective inertia diminishes as the pairing correlations become stronger than self-consistent values. The least-action path is determined dynamically by the interplay between these two inverse effects. Typical nature of a PES and M C along the pairing coordinate is shown in Figure 6 [50]. It portrays clear evidence of the opposite tendencies discussed above. Minimum action paths for two different nuclei are calculated in Sadhukhan et al. [50] by including the pairing degrees of freedom. Corresponding projections onto the (Q 20 , Q 22 ) and (Q 20 , λ 2 ) planes are shown in Figure 7. Also, two-dimensional (2D) fission paths calculated without pairing fluctuations are compared in this figure. In the case of 264 Fm, the three-dimensional (3D) pathway, calculated with pairing fluctuations, closely follows triaxial configurations of the corresponding 2D path. However, this scenario changes for 240 Pu, where the difference between the axial and triaxial barrier-heights is small and, as a result, pairing correlations could enforce the axial symmetry of the path in the region between the ground state configuration and superdeformed fission isomer. Nevertheless, irrespective of this system dependence, the pairing fluctuations are substantially enhanced in both the cases. Therefore, the dynamic coupling between pairing and deformation coordinates can produce dramatic changes in the SF process. Relativistic mean-field calculation [51] for Fm isotopes also shows a similar behavior of the fission pathway under the influence of dynamic paring fluctuations. Moreover, owing to the associated reduction in the action integral, the calculated SF half-life decreases by as much as three decades. These strong dynamical effects, predicted for SF, are however expected to disappear at higher excitation energies of the compound system. The subsection is concluded by demonstrating the average pairing gaps ( n and p ), in Figure 8, along the SF path of 240 Pu [29]. Although the dynamic (3D) path overlaps with the static pathway in the (Q 20 , Q 30 ) space, average pairing gaps are always higher in case of the 3D path except near the outer turning line where pairing correlations are quenched.

Fission-Fragment Yield Distributions
Stochastic Langevin dynamics is widely used to study fission fragment yield distributions of excited compound systems [89,91,106,107]. Only recently, it is successfully applied to calculate SF yields [29]. The detailed formalism is described in 1. Since the dynamics is stochastic in nature, it is difficult to understand the time evolution from a single fission event. Therefore, the concept of effective fission path (EFP) is devised in Sadhukhan et al. [92] for a better realization of the post-tunneling dynamics. First, for a given initial configuration, the local density of Langevin trajectories [108] is computed by counting the number of tracks in a small volume element of the collective space. Such distributions for two initial configurations are presented in Figure 9. Evidently, these two distributions are quite distinct in nature. The spreading of distribution is mainly governed by the interplay between the conservative and fluctuating forces. As I explained in 1, fluctuations become dominant near the scission and it leads to broader trajectory distributions. Next, the EFP is calculated by tracing the maxima in a trajectory-density distribution. Effectively, an EFP guides to the most probable fragmentation for the associated initial configuration. Eleven distinct EFPs are calculated in Sadhukhan et al. [92] for the representative system 240 Pu, and these are further shown in Figure 9. Also, the partial contribution  of each EFP to the total mass distribution is extracted by weighing the corresponding distribution with the appropriate tunneling probability. All the eleven partial mass distributions are plotted in Figure 10 along with the overall distribution. Due to higher tunneling probability, the peak region of the cumulative distribution is mostly contributed by EFPs close to the most probable path (EFP 5). On the other hand, contributions are negligible from those EFPs which originate far away from the most-probable outer turning point. For example, EFP 1 and EFP 11 hardly alters the total mass distribution. Most interestingly, certain EFPs (e.g., EFP 3 and EFP 4), associated with high tunneling probability, end up at large mass asymmetries and these constitute the tail part of the yield distribution. Such fission trajectories can only appear due to the presence of the random force in the Langevin dynamics.
In addition to the isolated yield distributions corresponding to either mass or charge of the fission fragments, DFT inputs enable the Langevin model to predict the correlation between mass and charge numbers of the fragments. In a recent study [66], it is calculated for the heaviest discovered element 294 118 Og. Three different EDFs are used for this purpose and, as shown in Figure 11, all of them predict a strongly asymmetric fission, or cluster emission, to be the dominant decay mode for this nucleus.

Uncertainties in Yield Distributions
It is necessary to estimate the uncertainties due to different model parameters [112]. In case of SF yield distributions, predicted within the hybrid WKB + Langevin method, uncertainties are primarily associated to the three input quantities: ground state zero-point energy E 0 , dissipation tensor η ij in Equation (5), and scission configuration which is defined with the neck-particle  [109,110] are shown by circles that include mirror points (open circles). Only the heavy-fragment parts are plotted in both (A,B). The figure is partially adapted from Sadhukhan et al. [92]. number N q . Moreover, the use of a particular EDF may induce additional bias in the results. In practice, the SF half-life is reproduced by tuning the free parameter E 0 [5,29] that effectively shifts the location of the inner and outer turning-points. As a result, both P(s out ) and fission paths are modified. Secondly, the fragment properties are expected to strongly depend on the scission configuration. However, no method exists that defines the scission configuration uniquely within the static adiabatic description of fission. Usually, different values of N q are used to identify the scission hypersurface in a multidimensional space. In the case of η ij , a microscopic theory is still missing and it is considered as an adjustable parameter in the SF model [29,66]. Considering all these limitations, a sensitivity analysis of the yield distributions with respect to all the model parameters is essential. Calculations are performed in Sadhukhan et al. [29] to illustrate the uncertainties in the yield distributions produced by E 0 , N q , and η ij . As demonstrated in Figure 12, the mass and charge distributions of 240 Pu are found to be robust against wide variations of all these input quantities. A similar response to the dissipation tensor is observed for the yield distributions of 294 118 Og [66]. Further, as plotted in Figure 13, both mass and charge distributions show weak dependency on the choice of the EDF and also on the dimensionality of the configuration space.   Figure 10. The shaded regions describe uncertainties in the yield distributions corresponding to different values of the model parameters: E 0 (narrow red band), dissipation tensor (wider cyan band), and scission configuration (linear hatch pattern). See Sadhukhan et al. [29] for details. The figure is adapted from Sadhukhan et al. [29]. Therefore, the hybrid model, with reasonable values of the input parameters, can be used for a reliable prediction of the SF yields. Of course, theoretical progresses toward a more precise calculation of input parameters is required [23].

Prediction of Fragments From
Localized Pre-fragment I have already argued how a pre-fragment can be defined from NLFs (see section 1). Properties of such prefragments are needed to be scrutinized very carefully to validate their applicability in predicting the fission fragments. To this end, the particle numbers of the prefragments are extracted along the EFPs shown in Figure 9 [92]. The results are reiterated here in Figure 14 and it displays that the pre-fragment particle-numbers remain remarkably stable as the deformation increases toward scission. Moreover, variations in these numbers, indicated by the bands in Figure 14, become fairly narrow (< ±2 particles) at large deformations. This suggests that the prefragments formed in an initial configuration of the fissioning nucleus hardly change. This early development of the prefragments is a manifestation of the freeze-out of single-particle energies along the fission pathway [25,39,56,113], since the system tries to retain its microscopic configuration by escaping the level crossings. The concept is further extended to predict the particle numbers of the fission fragments. This is accomplished by distributing the neck nucleons to each pre-fragments following a statistical prescription. The predicted yield distributions are found to agree well with the experimental data [67].
This fast and efficient method of fission-fragment identification will be very useful in the r-process network calculations that predict the astrophysical abundances of more than half of the elements heavier than iron. In a very neutron-rich environment, such at those exists in the ejecta of neutron-star mergers, neutron-induced and β-delayed fission are highly probable. It is speculated that nuclear fission terminates the r-process paths near A ∼ 300. The location of the r-process endpoint can significantly influence the final yield distribution [114]. Since these superheavy nuclei can not be studied with a standard laboratory procedure, reliable theoretical predictions are of utmost interest. Hence, the density functional formalism, which is deeply rooted in the underlying nucleon-nucleon interactions, provides an ideal platform. Moreover, there can be multiple fission cycles along the r-process path and these cycles involve a wide variety of fissioning nuclei. Therefore, the theoretical method needs to be very time-efficient. The model prescribed in Sadhukhan et al. [67] fulfills both the requirements.
In parallel, precise measurements of fission fragment yield distributions of neutron rich actinides and heavier elements may help to minimize the theoretical uncertainties associated with model parameters.

CONCLUSION
In this review, I have elaborated on a successful theoretical model for the spontaneous fission yields and lifetime. It is developed in a hybrid manner by combining the WKB approximation for quantum tunneling with the stochastic Langevin dynamics. Several recent advancements to enhance the predictive capabilities of this hybrid model are presented. In particular, I elucidated the intricate role of collective inertia and pairing correlations in guiding the dynamics during the tunneling phase. The inevitable presence of fluctuation and dissipation in the final stage of the fission dynamics is explained in connection with the calculation of fragment yields. This approach could be a prospective candidate for large-scale applications to a wide range of fissioning nuclei. In parallel, further improvements in different aspects of the model are in queue [23].
Although an accurate prediction of the fission observables is the foremost priority for a fission model, global calculations of fragment properties related to stellar nucleosynthesis processes additionally demand a faster and more reliable technique compared to existing models. This is because such calculations involve a large variety of fissioning nuclei most of which are outside the valley of nuclear stability. For this purpose, a quicker method [67] is recently proposed that utilizes the idea of shell-stabilized/localized prefragments. This model enables the identification of fragments by performing selfconsistent calculations within a very localized domain of the configuration space. Clearly, theoretical investigation of the fission process is a diverse field of research with a broad perspective.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this review work and has approved it for publication.