8 Computational strategy

A number of computer codes are available for modelling various aspects of fission dynamics. The most computationally intensive ones are those used for the calculation of HFB configurations and related quantities (inertias, energy overlaps, time-dependent evolution, etc.). Two basic implementations of HFB/TDHFB solvers differ in the representation of the single-particle wave functions:

Space Discretisation

This family of codes employs a finite volume, either in coordinate space or momentum space, that is discretised using a suitable mesh of points. Many different choices can be found in the literature, from the Lagrange mesh based on orthogonal polynomials Baye (2015), to B-spline-based Blazkiewicz et al. (2005); Pei et al. (2008), and adaptive wavelet-based methods Pei et al. (2014). Among the codes relevant for fission we mention: Sky3D Maruhn et al. (2014); Afibuzzaman et al. (2018); Schuetrumpf et al. (2018), and EV8 Ryssens, Hellemans, Bender and Heenen (2015) (publicly available) as well as HFB-AX Pei et al. (2008), MADNESS-HFB Pei et al. (2014), MOCCa Ryssens (2016), HFB-2D-LATTICE Terán et al. (2003); Blazkiewicz et al. (2005), Skyax Reinhard et al. (2020), and the code of Jin et al. (2017).

Basis Expansion

This family of codes is based on an expansion of single-nucleon wave functions in a finite set of suitable basis functions. Most often this is a basis of the harmonic oscillator (HO) eigenfunctions. There have also been attempts to use the transformed HO basis states Stoitsov et al. (2003), or a two-centre HO basis for improving the description of elongated shapes Dubray et al. (2008). Some of the principal codes that use this representation for fission modelling are: HFODD Schunck et al. (2017) and HFBTHO Perez et al. (2017) (both publicly available,) as well as HFBaxial Robledo (2010a); Robledo et al. (2018), HFBTri Robledo et al. (2018), and HFB3 Hashimoto (2013). The code HFB3 utilises a basis expansion for two Cartesian directions while employing a Lagrange mesh in coordinate space in the third. Two codes based on relativistic EDFs have recently been used in calculations of self-consistent mean-field configurations as input for modelling fission dynamics: DIRHB Nikšić et al. (2014) and MDC-RMF Lu et al. (2014).

The use of different representations determines the applicability of any given code. For instance, all numerical schemes can efficiently deal with zero-range (Skyrme-like) effective interactions; however, apart from implementations in spherical symmetry Bennaceur (2020), mesh-based discretisation schemes have not been able to employ finite-range effective forces so far. On the other hand, basis-expansion methods have a long history of calculations with (among others) various Gogny interactions Robledo et al. (2018), Coulomb Dobaczewski et al. (1996), and Yukawa Dobaczewski et al. (2009) forces, by expanding the interaction into Gaussian form factors. Extensions to more general finite-range interactions have also been proposed Parrish et al. (2013).

For the time-dependent calculations of fission dynamics, only mesh-based methods have been used so far. They satisfy the demands of extreme deformations encountered at scission, and also allow the description of the nascent fragments beyond scission. While techniques exist to optimise the choice of basis states, mesh-based calculations have the advantage that their numerical precision is essentially independent of the nuclear deformation Ryssens, Heenen and Bender (2015).

A major aspect of the difference between the two numerical schemes is the discretisation of the continuum. The first consequence of this is the treatment of pairing correlations since mesh- and basis-based methods give very different descriptions of positive-energy single-particle states. On the one hand, the coordinate representation takes correctly into account asymptotic properties of single-particle wave functions, but it becomes impractical and time-consuming when dense single-particle states extending to high energies are included. This is required for a full HFB description of the coupling between quasiparticle and quasihole states Dobaczewski et al. (1984). On the other hand, basis-expansion methods can easily manage states of arbitrarily high energy, but fail to reproduce the spatial asymptotics of wave functions, and are thus not appropriate for the description of weakly-bound systems. In view of the importance of pairing for fission applications, see section 3.10, this situation is not satisfactory. One should note, however, that the practical importance of asymptotic properties or high-energy states in the continuum has not yet been evaluated for fission.

Another aspect is the treatment of finite-temperature calculations. When the nucleus is heated the Fermi surface becomes more diffuse and the statistical mixture includes contributions from quasi-bound and unbound single-particle states. While preliminary studies have been reported Bonche et al. (1984); Zhu and Pei (2014); Schuetrumpf et al. (2016), the correct treatment of this degree of freedom is an open problem, see section 3.11. Nevertheless, as in the case of pairing correlations, we can already foresee that the applicability and performance of mesh-based and basis-based approaches will differ significantly.

Table 2 summarises the available codes for modelling deformed nuclei, their collective properties and fission dynamics. A few additional comments are in order: in the code Sky3D all spatial derivatives are evaluated using the finite Fourier transform method; the code of Hashimoto (2013) uses a Lagrangian coordinate-space grid in the direction of the axial-symmetry axis; the code of Jin et al. (2017) uses a complete basis of single-particle states in the solution of the TDHFB; code SkyAx (HFODD) can implement constraints on axial monopole, quadrupole, octupole, and hexadecapole deformations (non-axial deformations up to multipolarity $\lambda =9$) separately.

Even if the codes listed above offer a high level of flexibility and great potentiality, we recommend the development of the following new computing capabilities.

Adaptive meshes

The specific issues of fission dynamics require that nuclear wave functions are computed also in regions where the nuclear density vanishes. For mesh-based implementations, uniformly spaced grids thus include discretised continuum represented by a large number of single-particle states, even at fairly low energies. For this reason, we recommend the development of a new-generation of mesh-based codes that will utilise non-uniform meshes, with lattice points concentrated in space regions where nuclear densities are non-negligible. This method would require self-consistent redefinitions of meshes depending on relative distances, deformations, and relative orientations of nascent fragments.

Adaptive bases

By definition, implementations that utilise two-centre HO bases describe only regions of space where densities are sufficiently different from zero. However, they require adaptive methods to self-consistently define bases corresponding to relative distances, deformations, and relative orientations of nascent fragments, as proposed in Dobaczewski (2019). We recommend the development of the corresponding HO-basis codes. Another direction is to use the multi-resolution techniques with a multi-wavelet basis as in Pei et al. (2014).

General symmetry breaking

The concept of spontaneous symmetry breaking is crucial for a mean-field description of atomic nuclei, but it has not been exploited to its full capabilities yet. Even the most ambitious fission studies to date consider only configurations that still maintain certain self-consistent symmetries. In particular, time-reversal breaking configurations, needed for the description for odd and odd-odd nuclei as well as high-spin physics and multi-quasiparticle configurations, are not included in most available computer codes. To the best of our knowledge, HFODD and Sky3D are the only publicly available codes that include the degrees of freedom necessary for this type of studies.

GCM with time-odd momenta

Along the same line, an extension of current GCM codes to include time-odd collective momenta presents an interesting challenge. In addition to the implementation of time-reversal symmetry breaking, and the development of the actual GCM, this would also require the capability to consistently construct and constrain the relevant conjugate momentum operators.

Overlaps and Kernels

Both GCM and projection methods require multi-reference calculations, that is, determination of overlaps and matrix elements between different paired or unpaired product states. In addition to possible problems related to the density-dependence of the interaction Dobaczewski et al. (2007); Lacroix et al. (2009); Bender et al. (2009); Duguet et al. (2009); Robledo (2010b); Sheikh et al. (2019), such calculations always require a higher degree of symmetry breaking than those for the corresponding single-reference implementations. For example, restoration of the particle symmetry or 3D rotational symmetry implies the time-reversal or simplex symmetry breaking, respectively, even if the single-reference states that are subjected to projection conserve these symmetries. Therefore, it is recommended that new codes are initially developed with a maximum degree of symmetry-breaking capabilities, and then accelerated by implementing conserved symmetries in single-reference calculations. This will ensure that they are automatically portable to multi-reference frameworks.

Matrix Elements

To implement $K$-matrix reaction theory for fission rates, as discussed in section 3.9, the suite of computer programs should be augmented with routines that calculate effective Hamiltonian matrix elements between arbitrary CHF and CHFB configurations. Also, the calculation of decay widths require including momentum operators in the set of constraining fields in the CHF and CHFB codes.

TDDFT, TDGCM, ATDDFT, and QRPA codes

We recommend building new-generation codes for fission dynamics by directly implementing the capabilities of the TDDFT, TDGCM, ATDDFT, and QRPA methods. For example, implementations built on HO bases, proposed in Dobaczewski (2019), can be ported to the time-dependent adaptive bases for TDDFT, or to the iterative solutions of the ATDDFT and QRPA methods.

Mesh-based codes for non-local EDFs

Since it is unlikely that higher accuracy of the calculated nuclear observables can be obtained using local EDFs Kortelainen et al. (2014), current developments are focused on new-generation non-local EDFs. As discussed above, codes based on the HO basis are capable of treating such functionals fairly efficiently. It is, therefore, of paramount importance to develop algorithms for implementing the same capabilities in mesh-based codes. This is certainly a far-reaching goal; presently with no clear ideas on the direction to take. Nevertheless, we recommend that, because of its fundamental importance, a substantial effort should be devoted to attacking this problem.

Algorithmic improvements

All these developments recommended above depend on the efficient and robust generation of self-consistent mean-field solutions with many constraints. Constructing large numbers of these configurations is extremely demanding both in CPU time and in the time required to diagnose convergence issues. Besides the second-order gradient method Robledo and Bertsch (2011), advanced algorithms from the field of non-linear optimisation have been introduced to accelerate the convergence of self-consistent iterations Baran et al. (2008); Ryssens et al. (2019); the potential for further developments could be far greater. Supervised or deep learning techniques could also present significant opportunities to, e.g., optimise basis parameters for better numerical accuracy, perform real-time diagnostic about convergence, or provide good emulators of theoretical models Lasseri et al. (2020).

High-Performance Computing

Already in the late 2000s, nuclear fission was recognised as a scientific grand challenge justifying the development of exascale computer systems Young et al. (2009); Bishop and Messina (2009); Carlson et al. (2016). Many of the various recommendations discussed in this document, from large-scale potential energy surfaces with many degrees of freedom, to the coupling between TDHFB and TDGCM dynamics, will require the power of such facilities. Yet, this will require a serious effort by the community to modify, or re-factor, their codes in order to adapt to choices made at leadership computing facilities. Such choices include hardware architectures (GPU and hybrid chips, memory/core), software libraries (use of abstraction layers, more and more often in C++), or computing policies (limited runtime).

In conclusion, we recommend the development of numerical tools, which (i) target specific requirements of fission dynamics within the single-reference and multi-reference frameworks; (ii) are adapted to modern computing infrastructure; and (iii) build a common code base for fission theory. We advocate to increase the transparency associated with numerical choices by: (i) including a detailed description of numerical procedures in published studies (for benchmarking and an independent reproduction of the results); (ii) making codes publicly available under Open Source license along with their long-term continuous maintenance within, e.g., a git repository; and (iii) writing the codes in a sufficiently modular fashion, so that new advances can be more easily adopted by the community.

As discussed in section 8.1.2, the computational cost to generate and store the mean-field configurations and their related quantities can be extremely high in the context of fission applications. At the same time, many applications only require “integral” quantities related to these configurations. For example, computing spontaneous fission lifetimes in the WKB approximation, or fission fragment distributions within the TDGCM+GOA framework, only requires the HFB PES and the collective inertia tensor. Let us assume for the sake of the argument a 3-dimensional collective space with 200$\times$50$\times$50 = 500 000 points. Lifetimes, barrier penetrabilities, charge and mass distributions of the fission fragments can be computed merely from the knowledge of 1 scalar function of 3 variables (the energy) and a rank-2 tensor function of 3 variables (the inertia tensor). In our example, we would only need a grand total of 3.5$\times {10}^6$ function values, which takes a very small amount of storage. By contrast, storing all the information about the HFB solution across the same PES would take up in excess of 1.6$\times {10}^{13}$ function values (assuming a 40$\times$20$\times$20 box discretisation of 2000 HFB spinors). Having a database of such potential energy surfaces for nuclei, different energy functionals (Skyrme, Gogny, relativistic, different pairing functionals, etc.) and different DFT solver technologies would be very valuable to quantify theoretical uncertainties. Since the cost of generating a PES is high, it would also offer maximum leverage to our small community.