Microscopic self-consistent description of induced fission: dynamical pairing degree of freedom

Nuclear density functional theory (DFT) is the only microscopic framework that can be used over the entire table of nuclides in a self-consistent description of phenomena ranging from ground-state properties and collective excitations, to large-amplitude nucleonic motion, fission and low-energy collisions. In the case of nuclear fission Bender et al. (2020), in particular, a fully quantum mechanical many-body model can be constructed starting from the time-dependent generator coordinate method (TDGCM) Berger et al. (1991). In this approach the nuclear wave function is represented by a linear superposition of many-body generator states that are functions of collective coordinates. In most cases these coordinates parameterize the shape of the nuclear density. The Hill-Wheeler equation of motion determines the time evolution of the wave function in the restricted space of generator states Verriere and Regnier (2020). By employing the Gaussian overlap approximation (GOA), the GCM Hill-Wheeler equation reduces to a local, time-dependent Schrödinger equation in the space of collective coordinates. The microscopic input for the collective Schrödinger equation, that is, the nuclear potential and collective inertia, are determined by self-consistent mean-field calculations for a choice of the energy density functional (EDF) or effective interaction. The TDGCM+GOA method can be applied to the dynamics of induced fission, starting from the ground state and following the time evolution of collective degrees of freedom all the way to scission and the emergence of fission fragments. This framework has been very successfully implemented in a number of fission studies based on nonrelativistic Skyrme and Gogny functionals Berger et al. (1991); Verriere and Regnier (2020); Goutte et al. (2005); Regnier et al. (2016a, 2018, b, 2017); Verriere et al. (2017); Regnier et al. (2019); Zdeb et al. (2017). These studies have investigated the dependence of the predicted fission dynamics on the choice of the EDF, initial conditions, form of the collective inertia, and the definition of scission configurations.

Relativistic energy density functionals Vretenar et al. (2005); Meng et al. (2006); Meng (2016) have also been employed in the description of both spontaneous Zhao et al. (2015, 2016) and induced nuclear fission Tao et al. (2017); Zhao et al. (2019a, b, 2020). The microscopic input for these studies is generated using either the multidimensionally constrained relativistic mean-field (MDC-RMF) Lu et al. (2014) or the relativistic Hartree-Bogoliubov model (MDC-RHB) Zhao et al. (2017). By employing the TDGCM+GOA collective model, several interesting topics have been explored in this framework, such as the influence of static pairing correlations on fission yields, different approximations for the collective inertia tensor, and finite temperature effects.

Most applications of the TDGCM to fission dynamics have considered a two-dimensional space of collective coordinates such as, for instance, quadrupole and octupole shape degrees of freedom. The recently developed computer code FELIX Regnier et al. (2016a, 2018) offers the possibility of solving the TDGCM+GOA equation for an arbitrary number of collective variables. In Ref. Regnier et al. (2017) a preliminary calculation of induced fission dynamics of ²⁴⁰Pu isotope in the three-dimensional space of shape variables (quadrupole, octupole and hexadecupole intrinsic deformations) has been reported. In particular, this model can also be used for a quantitative analysis of the critical role of dynamical pairing correlations in the process of induced fission.

The importance of pairing correlations for both spontaneous and induced nuclear fission has been emphasized in a number of studies Negele et al. (1978); Łojewski and Staszczak (1999); Sadhukhan et al. (2014); Staszczak et al. (1985); Zhao et al. (2016); Bernard et al. (2019); Lazarev (1987); Sadhukhan et al. (2016); Bulgac et al. (2016); Qiang et al. (2021). For spontaneous fission it has been shown that the coupling between shape and pairing degrees of freedom has a pronounced effect on the calculated fission lifetimes Sadhukhan et al. (2014); Zhao et al. (2016). In particular, when the gap parameter is considered as a dynamical variable, pairing correlations are generally enhanced thus reducing the effective inertia and the action integral along the fission path. This effect can significantly reduce the estimated spontaneous fission lifetimes, and it has also been noted that pairing fluctuations can restore axial symmetry in the fissioning system. A study of induced fission of ²⁴⁰Pu, using the microscopic time-dependent superfluid local density approximation (TD SLDA) Bulgac et al. (2016), has shown that both shape and pairing modes determine the dynamics of the final stage of the fission process, from configurations close to the outer fission barrier to full scission.

The influence of ground-state (static) pairing correlations on charge yields and total kinetic energy of fission fragments for the case of induced fission of ²²⁶Th isotope was analyzed in Ref. Tao et al. (2017) using the TDGCM+GOA framework. It has been shown that an increase of the strength of the pairing interaction, beyond the range determined by empirical pairing gaps obtained from the experimental masses of neighboring nuclei, reduces the asymmetric peaks and enhances the symmetric peak in charge yields distribution. This is a very interesting result, and thus it is important to explore dynamical pairing correlations in induced fission. In this work we explicitly include the isoscalar pairing degree of freedom in the space of TDGCM+GOA collective coordinates, and perform the first realistic three-dimensional calculation of induced fission of ²²⁸Th. The theoretical framework and methods are reviewed in Sec. II. The details of the calculation and principal results are discussed in Sec. III. Section IV contains a short summary and outlook for future studies.

In the TDGCM+GOA framework induced fission is described as a slow adiabatic process determined by a small number of collective degrees of freedom. The initial step in modeling the fission of a heavy nucleus is a self-consistent mean-field (SCMF) calculation of the corresponding deformation energy surface as a function of few selected collective coordinates. Such a calculation provides the microscopic input, that is, the single-quasiparticle states, energies, and occupation factors, that determine the parameters of a local equation of motion for the collective wave function.

As in our previous studies, here we use the point-coupling relativistic energy density functional DD-PC1 Nikšić et al. (2008) in the particle-hole channel, while pairing correlations are taken into account in the Bardeen-Cooper-Schrieffer (BCS) approximation by a separable pairing force of finite range Tian et al. (2009):

where $\mathbf{R}=(\mathbf{r}_{1}+\mathbf{r}_{2})/2$ and $\mathbf{r}=\mathbf{r}_{1}-\mathbf{r}_{2}$ denote the center-of-mass and the relative coordinates, respectively. $P(\mathbf{r})$ reads

The parameters of the interaction were originally adjusted to reproduce the density dependence of the pairing gap in nuclear matter at the Fermi surface computed with the D1S parameterization of the Gogny force Berger et al. (1991). To reproduce the empirical pairing gaps in the mass region considered in the present study, the strength parameters of the pairing force have been increased with respect to the original values by the following factors: $G_{n}/G_{0}=1.12$ and $G_{p}/G_{0}=1.08$ for neutrons and protons, respectively.

In the present analysis the self-consistent deformation energy surfaces are calculated using the multidimensionally constrained relativistic mean-field (MDC-RMF) model Lu et al. (2012, 2014); Zhou (2016); Zhao et al. (2016) with constraints on mass multipole moments $Q_{\lambda\mu}=r^{\lambda}Y_{\lambda\mu}$, and the particle-number dispersion operator $\Delta\hat{N}^{2}=\hat{N}^{2}-\langle\hat{N}\rangle^{2}$. The Routhian is therefore defined as

where $E_{\rm RMF}$ denotes the total RMF energy that includes static BCS pairing correlations. The amount of dynamic pairing correlations can be controlled by the Lagrange multipliers $\lambda_{2\tau}$ ($\tau=n,p$), Vaquero et al. (2011, 2013); Sadhukhan et al. (2014); Zhao et al. (2016). To reduce the number of collective degrees of freedom and, therefore, the considerable computational task, here we only consider isoscalar dynamical pairing; $\lambda_{2n}=\lambda_{2p}\equiv\lambda_{2}$ is employed as the collective coordinate. The nuclear shape is parameterized by the deformation parameters

with $R=1.2~{}A^{1/3}$ fm. The shape is assumed to be invariant under the exchange of the $x$ and $y$ axes, and all deformations $\beta_{\lambda\mu}$ with even $\mu$ can be included simultaneously. The constrained RMF+BCS equations are solved by an expansion in the axially deformed harmonic oscillator (ADHO) basis Gambhir et al. (1990). In the present study calculations have been performed in an ADHO basis truncated to $N_{f}=20$ oscillator shells.

The dynamics of the fission process is thus governed by a local, time-dependent Schrödinger-like equation in the space of collective coordinates $\bm{q}$:

where $g(\bm{q},t)$ is the complex wave function of the collective variables $\bm{q}$ and time $t$. The collective Hamiltonian $\hat{H}_{\rm coll}(\bm{q})$

governs the time-evolution of the nuclear wave function from an initial state at equilibrium deformation, up to scission and the formation of fission fragments. $B_{ij}(\bm{q})$ and $V(\bm{q})$ denote the inertia tensor and collective potential, respectively, that are computed using the self-consistent solutions (total energy, single-quasiparticle states and occupation factors) for the RMF+BCS deformation energy surface. Here we assume axial symmetry with respect to the axis along which the two fragments eventually separate, and consider the three-dimensional (3D) collective space of quadrupole $\beta_{2}$ and octupole $\beta_{3}$ deformation parameters, and the dynamical pairing coordinate $\lambda_{2}$. The inertia tensor is the inverse of the mass tensor, that is, $B_{ij}(\bm{q})=(\mathcal{M}^{-1})_{ij}(\bm{q})$. The mass tensor is calculated using the adiabatic time-dependent Hartree-Fock-Bogoliubov (ATDHFB) method in the cranking approximation Baran et al. (2011):

where

$U$ and $V$ are the self-consistent Bogoliubov matrices, and $\rho$ and $\kappa$ are the corresponding particle and pairing density matrices, respectively. The derivatives of the densities are calculated using the Lagrange three-point formula for unequally spaced points Yuldashbaeva et al. (1999); Baran et al. (2011). The cranking expression Eq. (7) can be further simplified in the perturbative approach Brack et al. (1972); Nilsson et al. (1969); Girod and Grammaticos (1979); Bés and Szymański (1961); Sobiczewski et al. (1969), and this leads to the perturbative cranking mass tensor:

where

$|\mu\nu\rangle$ are the two-quasiparticle states with the corresponding quasiparticle energies $E_{\mu}$ and $E_{\nu}$. The details of the derivation of cranking formulas for the mass tensor can be found in Ref. Baran et al. (2011).

At this point we introduce an approximation that is not entirely consistent, but is necessary to reduce the computational task and stabilize the time-evolution of the collective state. Because we consider the particle-number dispersion operator $\Delta\hat{N}^{2}=\hat{N}^{2}-\langle\hat{N}\rangle^{2}$, the pairing part of the mass tensor must be calculated using the non-perturbative cranking expression Eq. (7). In the recent study of the differences between the perturbative and non-perturbative ATDHFB collective masses in the TDGCM+GOA description of induced fission Zhao et al. (2020), using the axial quadrupole and octupole intrinsic deformation as dynamical variables, we have shown that the structure of non-perturbative collective masses is much more complex due to changes in the intrinsic shell structure across the deformation energy surface, and is characterized by pronounced isolated peaks located at single-particle level crossings near the Fermi surface. It has been shown that the choice of non-perturbative cranking collective mass leads to a reduction of symmetric charge yields and, generally, to a better agreement with data. Even though in the analysis of Ref. Zhao et al. (2020) both non-perturbative and perturbative mass tensors were used in modeling induced fission dynamics, in the present 3D study the number of mesh points required to accurately calculate all the non-perturbative collective masses becomes prohibitively large for the available computational resources. The reason is the occurrence of pronounced peaks in the collective masses related to single-particle level crossings near the Fermi surface, and the corresponding abrupt changes of the occupation factors of single-particle configurations. This would require a major refinement of the grid, leading to very large number of mesh points and possible instabilities. Therefore, to be able to quantitatively analyze the effect of dynamical pairing on induced fission, here we use the perturbative cranking expression Eq. (9) to calculate the elements of the mass tensor that correspond to the quadrupole and octupole deformations while, as noted above, the pairing element of the mass tensor must be computed using the non-perturbative cranking formula. With the indices $1$, $2$ and $3$ corresponding to the $\beta_{20}$, $\beta_{30}$, and $\lambda_{2}$ collective coordinates, the following elements of the mass tensor are used: $\mathcal{M}^{Cp}_{11}$, $\mathcal{M}^{Cp}_{12}$, $\mathcal{M}^{Cp}_{22}$, $\mathcal{M}^{C}_{33}$, and we neglect the coupling terms $\mathcal{M}_{13}$ and $\mathcal{M}_{23}$.

To model the fission dynamics we follow the time-evolution of an initial wave packet $g(\bm{q},t=0)$ ($\bm{q}\equiv\{\beta_{2},\beta_{3},\lambda_{2}\}$), built as a Gaussian superposition of the quasi-bound states $g_{k}$,

where the value of the parameter $\sigma$ is set to 0.5 MeV. The collective states $\{g_{k}(\bm{q})\}$ are solutions of the stationary eigenvalue equation in which the original collective potential $V(\bm{q})$ is replaced by a new potential $V^{\prime}(\bm{q})$ that is obtained by extrapolating the inner potential barrier with a quadratic form. The mean energy $\bar{E}$ in Eq. (11) is then adjusted iteratively in such a way that $\langle g(t=0)|\hat{H}_{\rm coll}|g(t=0)\rangle=E_{\rm coll}^{*}$, and this average energy $E_{\rm coll}^{*}$ is set to be 1 MeV above the fission barrier. The TDGCM+GOA Hamiltonian of Eq. (6), with the original collective potential $V(\bm{q})$, propagates the initial wave packet in time.

The time propagation is modeled using the TDGCM+GOA computer code FELIX (version 2.0) Regnier et al. (2018). The time step is $\delta t=5\times 10^{-4}$ zs (1 zs $=10^{-21}$ s), and the charge and mass distributions are calculated after $4\times 10^{4}$ time steps, which correspond to 20 zs. As in our recent calculations of Refs. Tao et al. (2017); Zhao et al. (2019a, b, 2020), the parameters of the additional imaginary absorption potential that takes into account the escape of the collective wave packet in the domain outside the region of calculation Regnier et al. (2018) are: the absorption rate $r=20\times 10^{22}$ s^-1 and the width of the absorption band $w=1.0$.

The collective space is divided into an inner region with a single nuclear density distribution, and an external region that contains two separated fission fragments. The scission hyper-surface that divides the inner and external regions is determined by the expectation values of the Gaussian neck operator $\displaystyle\hat{Q}_{N}=\exp[-(z-z_{N})^{2}/a_{N}^{2}]$, where $a_{N}=1$ fm and $z_{N}$ is the position of the neck Younes and Gogny (2009). We define the pre-scission domain by $\langle\hat{Q}_{N}\rangle>3$, and consider the frontier of this domain as the scission surface. 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$. Each infinitesimal surface element is associated with a given pair of fragments $(A_{L},A_{H})$, where $A_{L}$ and $A_{H}$ denote the lighter and heavier fragments, respectively. The integrated flux $F(\xi,t)$ for a given surface element $\xi$ is defined as Regnier et al. (2018)

where $\bm{J}(\bm{q},t)$ is the current

The yield for the fission fragment with mass $A$ is defined by

The set $\mathcal{A}(\xi)$ contains all elements belonging to the scission hyper-surface such that one of the fragments has mass number $A$. The charge yields are obtained by convoluting the raw flux with a Gaussian function of the number of particles Regnier et al. (2016b); Zhao et al. (2019a), with a width of 1.6 units.

As an illustrative example, the fission of ²²⁸Th is considered. For this nucleus the charge distribution of fission fragments exhibits a coexistence of symmetric and asymmetric peaks Schmidt et al. (2000). In the first step a large scale MDC-RMF calculation is performed to generate the potential energy surface, single-nucleon wave functions and occupation factors in the $(\beta_{2},\beta_{3},\lambda_{2})$ collective space. The intervals for the values of the collective variables are: $-1\leq\beta_{2}\leq 7$ with a step $\Delta\beta_{2}=0.04$; $0\leq\beta_{3}\leq 3.5$ with a step $\Delta\beta_{3}=0.05$; and $-0.1\leq\lambda_{2}\leq 2.0$ with a step $\Delta\lambda_{2}=0.1$. The relativistic energy density functional DD-PC1 Nikšić et al. (2008) is used in the particle-hole channel, and particle-particle correlations are described by the separable pairing force (1) in the BCS approximation.

Figure 1 displays the 2D projections of the 3D deformation energy surface of ²²⁸Th on the collective plane $(\beta_{2},\beta_{3})$, for several values of the collective coordinate $\lambda_{2}$. Note that the value $\lambda_{2}=0$ corresponds to static BCS pairing, while positive $\lambda_{2}$ equates to enhanced pairing correlations. Only configurations with $Q_{N}\geq 3$ are plotted, and the frontier of this domain determines the scission contour. The red curves correspond to static fission paths of minimum total energy. The deformation surfaces for $\lambda_{2}=0$ and $\lambda_{2}=0.3$ are very similar, with a pronounced ridge separating the asymmetric and symmetric fission valleys. This ridge decreases with increasing values of the pairing coordinate $\lambda_{2}$. The scission contour for $\lambda_{2}=0$ starts from an elongated symmetric point at $\beta_{2}\approx 6$, and evolves to a minimal elongation $\beta_{2}\approx 3$ as reflection asymmetry increases. For larger values of $\lambda_{2}$ the scission contour is not modified significantly, although the starting point of the scission contour on the quadrupole axis shifts to smaller elongations at $\beta_{2}\approx 5$.

To illustrate the dynamical effect on the pairing correlations, in Fig. 2 we display the values of the neutron (upper panel) and proton (lower panel) pairing gap along the static fission paths, as functions of the quadrupole coordinate $\beta_{2}$, for several values of the isoscalar pairing collective coordinate $\lambda_{2}$ (cf. Fig. 1). In the interval of values of $\lambda_{2}$ considered here, the values of the pairing gaps increase by a factor $\approx 2-3$ and, characteristically, for $\lambda_{2}\geq 0.6$ all traces of the underlying shell effects along the static fission path vanish.

In Fig. 3 we plot the values of the perturbative cranking collective masses $M^{Cp}_{11}$, $M^{Cp}_{22}$, and the non-perturbative cranking mass $M^{C}_{33}$, along the static fission path as functions of the quadrupole deformation, for different pairing collective coordinates $\lambda_{2}$. As noted above, the indices $1$, $2$ and $3$ refer to the $\beta_{2}$, $\beta_{3}$ and $\lambda_{2}$ coordinates, respectively. The collective masses $M^{Cp}_{11}$ and $M^{Cp}_{22}$ exhibit a more complex structure for smaller values of deformations $\beta_{2}<1$ and decrease for larger deformations. One notices that $M^{Cp}_{22}$ exhibits a sharper decrease for $\beta_{2}>1$, compared to $M^{Cp}_{11}$. The behaviour of the $M^{C}_{33}$ collective mass is much more complex, due to the fact that it has been calculated using the non-perturbative cranking formula Eq. (7). This leads to the appearance of prominent peaks at the locations of single-particle level crossing near the Fermi surface. However, in general also $M^{C}_{33}$ decreases for large deformations and the peaks are less pronounced. As the value of the dynamical pairing coordinate $\lambda_{2}$ increases, all three collective masses decrease and the underlying shell effects gradually vanish. This effect is consistent with the results shown in Fig. 6 of Ref. Tao et al. (2017).

The scission contours in the ($\beta_{2},\beta_{3}$) plane are shown in Fig. 4 for several collective pairing coordinates $\lambda_{2}$. The contours are generally not very different, especially for asymmetric fission. In particular, scission points that are close to the static fission path are not sensitive to dynamical pairing. For larger values of $\lambda_{2}$, however, the scission contour is shifted towards smaller quadrupole deformations $\beta_{2}$ values for nearly symmetric fission.

In Fig. 5 we compare the theoretical predictions for the charge yields with the data for photo-induced fission of ²²⁸Th. The 3D calculation employs the collective space built from the deformation $\beta_{2}$, $\beta_{3}$ and pairing $\lambda_{2}$ coordinates, while the 2D calculation includes only the shape degrees of freedom $\beta_{2}$ and $\beta_{3}$, and static pairing correlations. As in our study of Ref. Tao et al. (2017), the calculation in the 2D collective space corresponds to normal and enhanced static pairing, that is, to 100% and 110% of the normal pairing strength, determined by the empirical pairing gaps. The theoretical predictions follow the general trend of the data, except that our model obviously cannot reproduce the odd-even staggering of the experimental charge yields. The calculation that includes only the 2D collective space with a static pairing strength adjusted to empirical ground-state pairing gaps in this mass region (100%), predicts yields that are entirely dominated by asymmetric fission with peaks at $Z=35$ and $Z=55$. By increasing static pairing (110%), the asymmetric peaks are reduced and a contribution of symmetric fission develops, but not strong enough to reproduce the data. It is interesting to notice that a very similar distribution of charge yields is predicted by the 3D model calculation that includes dynamical pairing. On a quantitative level, even the 3D calculation does not completely reproduce the experimental yields. The model predicts tails of the asymmetric peaks that are not seen in experiment, and thus fails to quantitatively match the symmetric contribution. It has to be noted, however, that in the present study the collective potential and inertia tensor have been calculated at zero temperature. In our recent study of finite temperature effects in TDGCM+GOA Zhao et al. (2019a), a calculation of induced fission of ²²⁶Th has shown that, although the model can qualitatively reproduce the empirical triple-humped structure of the fission charge and mass distributions already at zero temperature, the position of the asymmetric peaks and the symmetric-fission yield can be described much better when the potential and collective mass are determined at a temperature that approximately corresponds to the internal excitation energy of the fissioning nucleus.

Finally, to illustrate the effect of dynamical pairing on the flux of the probability current through the scission hyper-surface, in Fig. 6 we plot the time-integrated flux through the scission contour in the $(\beta_{2},\beta_{3})$ plane, for a given value of the pairing collective coordinate $\lambda_{2}$

The set $\mathcal{B}(\xi\equiv\beta_{2},\beta_{3})$ contains all elements of the scission contour with a given value $\lambda_{2}$. Even though it appears that dynamical pairing does not significantly modify the scission contour (cf. Fig. 4), nevertheless its effect on the collective flux and, therefore, on the occurrence of fission, is remarkable. For negative values of $\lambda_{2}$, that is, for correlations weaker than static pairing at $\lambda_{2}=0$, the flux rapidly decreases to zero. For positive values of $\lambda_{2}$ the flux exhibits a steep increase and a prominent peak at $\lambda_{2}\approx 0.3$. Note that this value corresponds to an increase of $\approx 20\%$ with respect to the static proton and neutron pairing gaps (cf. Fig. 2). The collective flux through the scission contour weakens with a further increase of pairing, and eventually vanishes for $\lambda_{2}>1$.

The influence of dynamical pairing degrees of freedom on induced fission has been investigated in a unified theoretical framework based on the generator coordinate method (GCM) with the Gaussian overlap approximation (GOA). In an illustrative calculation of fragment charge yields for induced fission of ²²⁸Th, the collective potential and inertia tensor have been computed using the self-consistent multidimensionally constrained relativistic mean field model, based on the energy density functional DD-PC1, and with pairing correlations treated in the BCS approximation with a separable pairing force of finite range. The fission fragment charge distributions are obtained by propagating the initial collective state in time with the time-dependent GCM+GOA. The flux of the probability current through the scission hyper-surface determines the probability of observing a given pair of fragments.

In this work, for the first time, the dynamics of induced fission has been consistently described in a three-dimensional space of collective coordinates that, in addition to the axial quadrupole and octupole intrinsic deformations of the nuclear density, also includes an isoscalar pairing degree of freedom. A number of studies has already demonstrated the importance of dynamical pairing for the calculation of spontaneous fission lifetimes, and static pairing correlations for modeling induced fission. As this work has also shown, a much more difficult problem is the inclusion of dynamical pairing degrees of freedom in a time-dependent description of induced fission. One expects, of course, that the model becomes more realistic as the Hilbert space of collective coordinates is expanded. However, as the present analysis has illustrated, it can be difficult to numerically stabilize the time-evolution of the fissioning system when shape and pairing collective coordinates are considered simultaneously in a three-dimensional calculation. Here it was necessary to compute the collective inertia using different cranking approximations for the shape and pairing degrees of freedom. Nevertheless, this approach enables a qualitative study of the effect of dynamical pairing on induced fission.

This initial study has clearly demonstrated the important effect that dynamical pairing correlations have on the induced-fission fragment distribution. In particular, the charge distribution of fission fragments of ²²⁸Th is characterized by symmetric and asymmetric peaks, but this structure cannot be reproduced in a two-dimensional calculation that only includes shape collective variables. In that case, and with a static pairing strength adjusted to empirical ground-state pairing gaps, the calculated yields correspond to a completely asymmetric fission. Only by artificially increasing the static pairing correlations or, more naturally, by including the dynamical pairing degree of freedom in the three-dimensional calculation, the asymmetric peaks get reduced and a contribution of symmetric fission develops in agreement with the empirical trend. It is also interesting to note that the time-integrated collective flux through the scission contour in the $(\beta_{2},\beta_{3})$ plane, exhibits a characteristic functional dependence on the pairing collective coordinate, with a prominent peak at a value that correspond to an increase of $\approx 20\%$ with respect to the static pairing gaps.

Future advances in computational capabilities will open the possibility of more quantitative applications of multi-dimensional TDGCM+GOA to fission dynamics. An immediate task will be to consider shape and pairing degrees of freedom on an equal footing, and consistently compute the corresponding collective inertia tensor and metric. A more challenging problem is to include dynamical pairing degrees of freedom in recently developed fission models that attempt to incorporate restoration of symmetries broken by the intrinsic densities in constrained mean-field calculations (rotational, reflection, and particle number symmetry) Marević and Schunck (2020); Verriere et al. (2019, 2021).