Neural Network Emulation of Spontaneous Fission
Large-scale computations of fission properties are an important ingredient
for nuclear reaction network calculations simulating rapid neutron-capture
process (the r process) nucleosynthesis. Due to the large number of fissioning
nuclei contributing to the r process, a microscopic description of fission
based on nuclear density functional theory (DFT) is computationally
challenging. We explore the use of neural networks (NNs) to construct DFT
emulators capable of predicting potential energy surfaces and collective
inertia tensors across the whole nuclear chart. We use constrained
Hartree-Fock-Boguliubov (HFB) calculations to predict the potential energy and
collective inertia tensor in the axial quadrupole and octupole collective
coordinates, for a set of nuclei in the r-process region. We then employ NNs to
emulate the HFB energy and collective inertia tensor across the considered
region of the nuclear chart. Least-action pathways characterizing spontaneous
fission half-lives and fragment yields are obtained using the nudged elastic
band method. The potential energy predicted by NNs agrees with the DFT value to
within a root-mean-square error of 500 keV, and the collective inertia
components agree to within an order of magnitude. The exit points on the outer
turning line are found to be well emulated. For the spontaneous fission
half-lives the NN emulation provides values that are found to agree with the
DFT predictions within a factor of $10^3$ across more than 70 orders of
magnitude. Neural networks are able to emulate the potential energy and
collective inertia well enough to reasonably predict physical observables.
Future directions of study, such as the inclusion of additional collective
degrees of freedom and active learning, will improve the predictive power of
microscopic theory and further enable large-scale fission studies.
Daniel Lay, Eric Flynn, Samuel A. Giuliani, Witold Nazarewicz, Leó Neufcourt
# Neural Network Emulation of Spontaneous Fission

###### Abstract

**Background:** Large-scale computations of fission properties are an important ingredient for nuclear reaction network calculations simulating rapid neutron-capture process (the \(r\) process) nucleosynthesis. Due to the large number of fissioning nuclei potentially contributing to the \(r\) process, a microscopic description of fission based on nuclear density functional theory (DFT) is computationally challenging.

**Purpose:** We explore the use of neural networks (NNs) to construct DFT emulators capable of predicting potential energy surfaces and collective inertia tensors across the whole nuclear chart, starting from a minimal set of DFT calculations.

**Methods:** We use constrained Hartree-Fock-Boguliubov (HFB) calculations to predict the potential energy and collective inertia tensor in the axial quadrupole and octupole collective coordinates, for a set of nuclei in the \(r\)-process region. We then employ NNs to emulate the HFB energy and collective inertia tensor across the considered region of the nuclear chart. Least-action pathways characterizing spontaneous fission half-lives and fragment yields are then obtained by means of the nudged elastic band method.

**Results:** The potential energy predicted by NNs agrees with the DFT value to within a root-mean-square error of 500 keV, and the collective inertia components agree to within an order of magnitude. These results are largely independent of the NN architecture. The exit points on the outer turning line are found to be well emulated. For the spontaneous fission half-lives the NN emulation provides values that are found to agree with the DFT predictions within a factor of \(10^{3}\) across more than 70 orders of magnitude.

**Conclusions:** Neural networks are able to emulate the potential energy and collective inertia well enough to reasonably predict physical observables. Future directions of study, such as the inclusion of additional collective degrees of freedom and active learning, will improve the predictive power of microscopic theory and further enable large-scale fission studies.

## I Introduction

Large scale calculations of fission properties are an essential ingredient for the modelling of the rapid neutron-capture process (\(r\) process), responsible for the production of roughly half of the nuclei heavier than iron found in nature [1; 2]. Fission determines the range of the heaviest nuclei that can be synthesized during the \(r\)-process, recycles the material during the neutron irradiation phase, and shapes the final abundances [3; 4; 5]. Given the large amount of energy released in this decay, the presence of fissioning nuclei can leave fingerprints in the electromagnetic counterpart produced in neutron star mergers [6; 7]. However, as most of the fissioning nuclei produced during the \(r\) process cannot be measured, theoretical predictions are indispensable to perform accurate nuclear reaction network calculations.

During the last decades, several efforts have been devoted to the systematic estimation of fission barriers [8; 9; 10; 11; 12; 13], spontaneous fission half-lives [14; 15; 13], and fragment distributions [16; 17; 18; 19; 20] of \(r\)-process nuclei. However, due to the inherent complexities characterizing the theoretical description of the fission process [21], most of the available calculations resort to phenomenological approaches based on simplified assumptions. This limitation can be overcome by employing the nuclear DFT [22; 23; 24], which is the quantum many-body method based on effective nucleon-nucleon interactions applicable across the whole nuclear landscape. But given its computational costs, using DFT for fission is a daunting task for large-scale studies of \(r\)-process nuclei [25; 21; 26]. As such, the usage of DFT emulators can be an invaluable tool to extend the current reach of microscopic fission calculations.

Machine learning has been used with great success in many areas of nuclear physics (see [27] for a recent review on this topic). In particular, machine learning has been used in many DFT studies to emulate potential energy surfaces (PESs), in both quantum chemistry [28; 29; 30; 31] and in nuclear physics [32; 33]. However, these have generally focused on emulating individual potential energy surfaces, rather than many nuclei across a portion of the nuclear chart (or many related chemical systems in the quantum chemistry case). In an important study, Ref. [34] succeeded in emulating PESs and other quantities using committees of multilayer neural networks.

In this study, we use fully connected, feedforward NNs to emulate the PES and collective inertia tensor, parameterized by the axial quadrupole and octupole moments
\(Q_{20}\) and \(Q_{30}\), between nuclei in the \(r\)-process region of the nuclear chart. The paper is organized as follows: Section II reviews the theoretical approach to spontaneous fission used in this work. Section III describes the characteristics of the employed NNs. Section IV demonstrates the performance of the NNs on the HFB energy and collective inertia tensor, and Sec. V compares the exit points and spontaneous fission half-lives obtained using the DFT inputs and the emulated NN inputs. Finally, conclusions are summarized in Sec. VI.

## II Spontaneous fission within the nuclear density functional theory

Spontaneous fission (SF) is a dynamical process where the nucleus evolves from the ground-state into a split configuration. In the adiabatic approximation, SF is modeled using a finite set of collective variables \(\{q_{i}\}\) usually describing the nuclear shape. The SF half-life can be computed within this approach as \(t_{1/2}=\ln 2/nP_{\text{fs}}\), where \(n\) is the number of assaults on the fission barrier, and \(P_{\text{fs}}\) the fission probability given by the probability of the nucleus to tunnel through the fission barrier, which can be estimated using the semiclassical Wentzel-Kramers-Brillouin (WKB) approach [35]:

\[P_{\text{fs}}=\frac{1}{1+\exp{(2S(L))}}, \tag{1}\]

where \(S(L)\) is the collective action computed along the stationary trajectory \(L[s]\) that minimizes \(S\) in the multi-dimensional space defined by the collective coordinates:

\[S(L[s])=\frac{1}{\hbar}\int_{s_{\text{in}}}^{s_{\text{out}}}\sqrt{2\mathcal{ M}_{\text{eff}}(s)(V(s)-E_{0})}\;ds\,, \tag{2}\]

with \(V\) and \(\mathcal{M}_{\text{eff}}\) being the potential energy and inertia tensor, respectively, computed along the fission path \(L[s]\). The integration limits \(s_{\text{in}}\) and \(s_{\text{out}}\) correspond to the classical inner and outer turning points, respectively, defined by the condition \(V=E_{0}\), where \(E_{0}\) is the collective ground-state zero-point energy stemming from quantum fluctuations in the collective coordinates. While the latter can be estimated from, e.g., the curvature of \(V\) around the ground state (g.s.) configuration, in many SF studies \(E_{0}\) is taken as a fixed positive constant ranging between 0.5 and 2.0 MeV above the ground-state energy [15; 13]. For simplicity, we follow the latter approach and fix \(E_{0}=E_{\text{g.s.}}\). And, throughout this work, we will refer to the collective coordinates at \(s_{\text{out}}\) (in this work, \((Q_{20},Q_{30})\)) as the exit point [36].

From Eq. (2) it can be deduced that the main ingredients required for the estimation of the SF half-lives are the effective potential energy \(V\) and collective inertia \(\mathcal{M}_{\text{eff}}\). In this work, we compute these quantities by employing the self-consistent mean-field method [22; 24] summarized in the following.

Nuclear configurations are obtained by means of the HFB method, where the many-body wave function \(|\Psi\rangle\), described as a generalized quasiparticle product state, is given by the minimization of the mean value of the Routhian:

\[\widehat{\mathcal{H}}^{\prime}=\widehat{\mathcal{H}}_{\text{HFB}}-\sum_{\tau= n,p}\lambda_{\tau}\widehat{N}_{\tau}-\sum_{\mu=1,2,3}\lambda_{\mu}\widehat{Q}_{ \mu 0}\,. \tag{3}\]

In Eq. (3), \(\widehat{\mathcal{H}}_{\text{HFB}}\) is the HFB Hamiltonian, and \(\lambda_{Z}\) and \(\lambda_{N}\) are the Lagrange multipliers fixing the average number of protons and neutrons, respectively. The shape of the nucleus is enforced by constraining the moment operator \(\widehat{Q}_{\mu\nu}\) with multipolarity \(\mu\) and magnetic quantum number \(\nu\). In this work, we explore the evolution of the total energy and collective inertia tensor as a function of the elongation of the nucleus and its mass asymmetry, which are described by the axial quadrupole \(Q_{20}\) and octupole \(Q_{30}\) moment operators, respectively:

\[\widehat{Q}_{20} =\hat{z}^{2}-\frac{1}{2}(\hat{x}^{2}+\hat{y}^{2})\,; \tag{4a}\] \[\widehat{Q}_{30} =\hat{z}^{3}-\frac{3}{2}(\hat{y}^{2}+\hat{x}^{2})\hat{z}\,. \tag{4b}\]

In order to reduce the computational cost, axial symmetry is enforced in all the calculations (\(\langle\widehat{Q}_{\mu\nu}\rangle=0\) for all \(\nu\neq 0\)), and the additional constraint \(\langle\widehat{Q}_{10}\rangle=0\) is imposed to remove the spurious center-of-mass. Finally the nuclear HFB Hamiltonian \(\widehat{\mathcal{H}}_{\text{HFB}}\) is given by the finite-range density-dependent nucleon-nucleon Gogny interaction. We employ the D1S parametrization [37], which has been widely used in nuclear structure studies across the whole nuclear chart including the description of fission properties of heavy and superheavy nuclei [38]. The effective potential is then given by \(V=E-E_{\text{rot}}\), where \(E\) is the energy obtained from the HFB equations for the Routhian (3), and \(E_{\text{rot}}\) is the energy correction related to the restoration of rotational symmetry, computed using the approach of Ref. [39]. Calculations are carried out by employing the HFB solver HFBaxial, which solves the HFB equations by means of a gradient method with an approximate second-order derivative [40]. The quasiparticle wave functions are expanded in an axially-symmetric deformed harmonic oscillator single-particle basis, containing states with \(J_{z}\) quantum number up to 35/2 and up to 26 quanta in the \(z\)-direction. The basis quantum numbers are restricted by the condition \(2n_{\perp}+|m|+n_{z}/q\leq N_{z}^{\text{max}}\), where \(q=1.5\) and \(N_{z}^{\text{max}}=17\). This choice of the basis parameters allows for a proper description of the elongated prolate shapes characteristic of the fission process [41].

The collective inertia tensor \(\mathcal{M}_{\mu\nu}\) is computed within the Adiabatic-Time-Dependent HFB (ATDHFB) approximation using the non-perturbative scheme [42; 43; 44]:

\[\mathcal{M}_{\mu\nu}=\frac{\hbar^{2}}{2\hat{q}_{\mu}\hat{q}_{\nu}}\sum_{\alpha \beta}\frac{F_{\alpha\beta}^{\mu*}F_{\alpha\beta}^{\nu}+F_{\alpha\beta}^{\mu}F _{\alpha\beta}^{\nu*}}{E_{\alpha}+E_{\beta}}\,, \tag{5}\]
where \(q_{i}\) are the collective coordinates and

\[\frac{F^{\mu}}{\dot{q}_{\mu}}=A^{\dagger}\frac{\partial\rho}{\partial q_{\mu}}B^{ *}+A^{\dagger}\frac{\partial\kappa}{\partial q_{\mu}}A^{*}-B^{\dagger}\frac{ \partial\rho^{*}}{\partial q_{\mu}}A^{*}-B^{\dagger}\frac{\partial\kappa^{*}}{ \partial q_{\mu}}B^{*} \tag{6}\]

is given in terms of the matrices of the Bogoliubov transformation \(A\) and \(B\), and the corresponding particle \(\rho\) and pairing \(\kappa\) densities. Then, the effective inertia tensor is given as

\[\mathcal{M}_{\rm eff}=\sum_{\mu\nu}\mathcal{M}_{\mu\nu}\frac{dq_{\mu}}{ds} \frac{dq_{\nu}}{ds}\,. \tag{7}\]

It is important to remark that the \(\mathcal{M}_{\mu\nu}\) components can suffer from rapid oscillations in the presence of single-particle level crossings near the Fermi surface. Such abrupt changes of occupied single-particle configurations produce variations in the derivatives of the densities in Eq. (6), resulting in pronounced peaks of \(\mathcal{M}_{\rm eff}\) along the fission path [43; 44; 45].

The least action paths are computed using the nudged elastic band method (NEB) [36]. Due to the large number of paths that must be explored, NEB parameters cannot be tuned by hand. Instead, multiple NEB runs are started, with initial paths ending at various points along the outer turning line. The NEB algorithm depends on two parameters, \(k\) and \(\kappa\), which adjust spring and harmonic restoring forces, respectively. Not varying \(k\) and \(\kappa\) will, on occasion, miss some LAPs, akin to skipping over a narrow minimum in an optimization routine. Different runs are started for \(k\) and \(\kappa\) in the range \(0.05-10\), for each initial path. These runs converge to a number of different stationary paths. Typically, there is some component of the path that travels along the outer turning line. To select the final tunneling path, the paths are interpolated using 500 points, and truncated when near the outer turning line and within an energy tolerance of 0.5 MeV. The unique paths are chosen based on the clustering of the exit point using the mean shift algorithm as implemented in scikit-learn [46], and the path corresponding to a given exit point with the least action is chosen as the LAP.

## III Neural networks

In this work, we use feedforward NNs as our emulators. We train separate NNs on the potential energy \(V\) and the components of \(\mathcal{M}\). Each NN takes as input \((A,Z,Q_{20},Q_{30})\), specifying the nucleus and deformation, then outputs the value (either \(V\) or one of \(\mathcal{M}_{\mu\nu}\)) at that point. As discussed in Sec. IV.1, to further improve NN performance, we rescale the NN inputs to lie between zero and one.

We train NNs with a number of hidden layers varying between 2 and 7, with 200 hidden nodes in the first layer, and a decreasing number of nodes in each subsequent layer. We use the RELU activation function, and train to minimize the root-mean-square error in the desired quantity. For each variant on the NN depth, we train multiple NNs, forming a committee of NNs. We then combine the predictions from each NN in the committee in a weighted average, to further reduce the error on the prediction.

To train the NNs, we have computed PESs and the collective inertia for 194 nuclei, each on a regular grid of 4 b for \(0\leq Q_{20}\leq 248\) b, and 6 b\({}^{3/2}\) for \(0\leq Q_{30}\leq 60\) b\({}^{3/2}\). These nuclei are then labeled as either training, combining, or validation, with the latter two sampled randomly across the chart. These different datasets are indicated in Fig. 1. For each nucleus, the entire grid is used in the training/combining/validation. The nuclei in the training set are used to train individual NNs, the nuclei in the combining dataset are used to combine predictions from the committee members in a weighted average, and the nuclei in the validation set are used for validation of the NN predictions. The weights for each committee member are chosen to minimize the root-mean-square error on the nuclei in the combining dataset. As can be seen, most of the nuclei (about 70%) are used for training, with the remaining 30% split equally between the combining and validation datasets. In general, the NN performance is not sensitive to the distribution of training data, provided the NN does not attempt to extrapolate across the nuclear chart. No detailed optimization of the choice of training nuclei was carried out.

As mentioned in Sec. II, it is known that the collective inertia tensor can develop discontinuities and rapid variations due to level crossings. This makes emulation of the tensor challenging since the tensor components can span many orders of magnitude as a function of deformation. If the NN is trained on the inertia tensor components by themselves, the network predictions are poor. However, while these problems are features of the approximations used to calculate the inertia, the NN can still learn certain features of the inertia tensor by carrying out the eigenvalue decomposition of the inertia tensor,

\[\mathcal{M}=U\Sigma U^{T} \tag{8}\]

where \(U\) is the \(2\times 2\) matrix of eigenvectors and \(\Sigma\) is the diagonal matrix of eigenvalues. Since \(U\) is an orthogonal matrix, we can represent \(U\) as an element of the set SO(2) parameterized by Euler angle \(\theta\). In this representation, \(\mathcal{M}\) is completely parameterized by its eigenvalues and the Euler angle \(\theta\). So, the NN is trained on \(\theta\) and the log of the eigenvalues at each point (\(Q_{20}\), \(Q_{30}\)). Training on this representation of the tensor is similar to normalizing the network inputs, as both put NN inputs/outputs on a similar scale. Additionally, this forces the tensor predictions to be positive semi-definite. We also transform \(\theta\) to the range \((-\pi/2,\pi/2)\), so that the angles are mostly clustered near zero (on the interval \((0,\pi)\), there are two clusters: one at 0 and one at \(\pi\), which the NN has difficulties learning).

Once the NNs are trained, PESs and inertias are com
puted for the same grid of deformations as the original DFT calculations. While the NNs can be evaluated at arbitrary \((Q_{20},Q_{30})\), it is less computationally expensive to use a standard cubic spline interpolator on the grid predicted by the NN. Moreover, the LAPs computed using the NN evaluations and the spline interpolator agree well with each other. Due to the relatively large number of LAP calculations, we report the LAPs computed using the spline-interpolated NN predictions, rather than using the NN predictions themselves.

## IV Neural Network Quality

Here, we examine the quality of the NNs, on both the PES and the collective inertia. In general, we observe that the NN is able to reproduce both the PES and the collective inertia for most of the nuclei under consideration. Moreover, the quality of the NN is relatively stable across the different architectures considered.

Throughout this section, we will refer to the PES and collective inertia computed using DFT as the reference data, and the PES and inertia computed using the NN as the NN reconstruction.

### Potential energy surfaces

For a single nucleus, we define the root-mean-square error (RMSE) \(\Delta V(A,Z)\) in energy over the collective domain considered as

\[\Delta V(A,Z)^{2}=\frac{1}{n}\sum_{Q_{20},Q_{30}}[V^{\rm DFT}(Q_{ 20},Q_{30},A,Z)\\ -V^{\rm NN}(Q_{20},Q_{30},A,Z)]^{2}, \tag{9}\]

where \(n=693\) is the number of gridpoints evaluated in the PES. A similar quantity can be defined for the components of \(\mathcal{M}_{\rm eff}\), although there \(n\) varies slightly from nucleus to nucleus.

Figure 1 shows \(\Delta V(A,Z)\) across the region of the nuclear chart considered, for the deepest NN (7 hidden layers, with 200-175-150-125-100-75-50 hidden units), with rescaled inputs. As can be seen, for most nuclei, \(\Delta V(A,Z)\lesssim 0.5\) MeV. Exceptions occur, with most remaining below 1.5 MeV. For some nuclei, such as \({}^{308}\)Cf, \({}^{314}\)Fm, and \({}^{318}\)No, relatively poor performance may be expected: these nuclei are on the outer edge of the region of the nuclear chart considered, and hence the NN is extrapolating from the training region to reach them. For other nuclei, such as \({}^{232}\)Th and \({}^{280}\)Cm, poor performance is unexpected: these nuclei are surrounded by training nuclei, and so should be emulated fairly well. As such, it seems unlikely that poor performance is due solely to the location of the nucleus on the nuclear chart relative to the training data.

To understand the reduced performance, we examine the nuclei in question. Figure 2 shows both the reference PES and its NN reconstruction for \({}^{280}\)Cm. This nucleus is chosen because it has \(\Delta V=2.15\) MeV, which is the largest of all nuclei in the validation set. While the PESs are not identical, features such as the ground state and outer turning line locations, as well as the general shape of the PES, agree quite well. Large discrepancies tend to be limited to the high-energy region of \(Q_{20}\approx 50\) b, \(Q_{30}\gtrsim 30\) b\({}^{3/2}\), which is hardly relevant to fission. We conclude therefore that even for nuclei with larger RMSE, NNs could provide a very reasonable description of the fission path. This aspect will be examined further in Sec. V.

To assess the sensitivity of our results with respect to the NN architecture, we repeated our calculations employing different NN sizes and rescaling the inputs. Figure 3 shows the RMSE (now averaged across all nuclei in a given dataset) across the different datasets for a variety of NN depths. As is generally expected, the training dataset has a monotonically decreasing RMSE as the NN depth increases; this is simply due to the increasing number of tunable parameters in the NN. On the other hand, the RMSE for the combining and validation sets is fairly stable with respect to the number of hidden layers of the NN.

A general improvement is observed when normalizing the inputs \((A,Z,Q_{20},Q_{30})\) to be between 0 and 1. This is due to two factors. First,

Figure 1: \(\Delta V(A,Z)\) (in MeV) for the deepest NN. The different shapes indicate which dataset each nucleus belongs to.

Figure 2: The reference PES (left) and the NN reconstruction (right) for \({}^{280}\)Cm, in MeV. The ground state is marked with a \(\times\) symbol.
invariant. An input much larger than 1 is equivalent to a large initial weight, with a normalized input. Because the final NN weights are expected to be small (hence initializations following e.g. the Xavier initialization [47], as in this work), the initial weights are far from the final values, and convergence slows. Second, the optimization method itself is not scale-invariant: non-normalized inputs correspond to an ill-conditioned Hessian matrix, in which case gradient descent (and related methods) converge slowly [48; 49].

We conclude that the NN performance in predicting the PES is relatively stable with respect to the NN architecture; Sec. V will demonstrate that performance on this level is adequate for predicting SF observables.

### Collective Inertia

Since the components of \(\mathcal{M}\) vary across multiple orders of magnitude and the network is trained on the log of the eigenvalue decomposition, a loss function such as the root-mean-squared error is not an adequate measure of the performance of the NN. Instead, Fig. 4 shows the reference inertia components, plotted against the NN reconstructions, for all nuclei considered. The NN used is the 7-hidden-layer NN with rescaled inputs, with the number of hidden units as described in Sec. IV.1. The diagonal components \(\mathcal{M}_{22}\) and \(\mathcal{M}_{33}\) are predicted fairly well, as the distributions align roughly along the diagonal. It is worth noting that the distributions are slightly misaligned, in all data sets considered, indicating that the NN tends to underpredict relatively large values, and overpredict relatively small values. This, in turn, shows that the NN is slightly biased towards the mean value of the inertia.

However, the off-diagonal component \(|\mathcal{M}_{23}|\) is not aligned along the diagonal, except for large values. This is because this component actually varies across almost 10 orders of magnitude (compared to the 4 orders of magnitude for \(\mathcal{M}_{22}\) and \(\mathcal{M}_{33}\)), and so the NN is biased towards predicting the larger values more accurately, resulting in a general overprediction of \(\mathcal{M}_{23}\). In terms of the angle \(\theta\) that is actually determined by the NN, it is difficult to predict both small and large angles, and because \(\theta\) is allowed to be negative, a logarithm transform is not possible. Nevertheless, one obtains a reasonable-looking distribution above \(|\mathcal{M}_{23}|\gtrsim 10^{-4}\,\text{MeV}^{-1}\,\text{b}^{-5/2}\), indicating that some learning has indeed taken place. And, the poorly-learned values below \(10^{-4}\,\text{MeV}^{-1}\,\text{b}^{-5/2}\) are truncated at values \(10^{-6}-10^{-2}\,\text{MeV}^{-1}\,\text{b}^{-5/2}\).

When changing the depth of the NN, performance is similar. For shallow networks, predictions on the training dataset show a larger bias: the distribution of points on the inertia plot is less aligned with the diagonal for the \(\mathcal{M}_{22}\) and \(\mathcal{M}_{33}\) components. In other words, the larger reference values are underestimated, and the smaller reference values are overestimated. The validation dataset is aligned similar to the deepest network, shown in Fig. 4. As the depth of the network is increased, the training data points are aligned closer with the diagonal. This is indicative of the NN tending to overfit on the training data as the number of variational parameters increases. The distribution of \(\mathcal{M}_{23}\) values remains approximately

Figure 3: The RMSE for a variety of different NN sizes. The dashed line shows the same depth NN, but with input variables normalized to the range \([0,1]\).
the same when increasing NN depth, with a slight improvement on the truncated \(\mathcal{M}_{23}\) values. In general, the NN performance on the validation dataset is mostly stable when varying the NN depth. The overarching question is whether this performance is sufficient for predicting observable quantities of interest. As with the PES, this question can be directly answered by looking at NN predictions of physical observables.

## V Impact on Observable Quantities

While encouraging, the results discussed in Sec. IV do not give a perfectly clear estimation of the performance of the NNs. For instance, the NN reconstruction of the PES for \({}^{280}\)Cm may be adequate for reproducing fission observables - especially SF fragment yields and half-lives, despite the poor RMSE, since the largest deviations occur at deformations that will not be explored by LAPs. Similarly, the NN commonly fails to reproduce the off-diagonal component of the collective inertia, \(\mathcal{M}_{23}\), but primarily for small values of \(\mathcal{M}_{23}\).

Here, we examine the performance of the NN on the lifetime-weighted exit point, as a proxy for the fragment yield [20, 50], and the half-life of the nucleus. For both quantities, we compare three sets of data: the quantity computed using (i) the reconstructed PES and the identity inertia; (ii) using the reference PES and the reconstructed inertia; and (iii) the reconstructed PES and inertia. In this way, we can isolate the impact of the PES and inertia emulations separately, and combine them to assess the overall error of the emulator. In this section, we use the 7-hidden-layer NN with rescaled inputs, with a number of hidden units described in Sec. IV.1. Based on the relative insensitivity to the depth of the NN shown in Sec. IV, the overall performance is expected to be similar for different NN depths.

Similar to Sec. IV, exit points and SF half-lives computed using only DFT inputs will be referred to as reference quantities; those with any NN input will be referred to as reconstructed quantities.

### Exit Points

As demonstrated in [20, 50], the location of the exit points is sufficient for roughly estimating the fission fragment yields. For this reason, we can consider exit points as reasonable proxies for the fragment yields. When multiple fission channels exist, the combined fragment yields are attained by adding the yields of each channel, weighted by the probability of populating a particular channel. Thus, agreement of the lifetime-weighted exit point indicates strong agreement in the fission fragment yields (and, by necessity, indicates that the dominant fission mode is also in agreement between the reference data and the NN reconstruction).

Figure 5 shows the difference in the octupole moment of the lifetime-weighted exit point, for configuration (iii) mentioned above. The octupole moment is chosen because it is critical for explaining multimodality in SF. The \(Q_{30}\) error is similar for the other configurations, and the quadrupole moment is typically within \(\pm 1\,\mathrm{b}\) for all configurations.

The agreement is good between the reference exit point and the NN reconstruction: at \(\pm 1\)\(\mathrm{b}^{3/2}\), we expect the fragment yields to agree well (within the hybrid method of Refs. [20, 50]). This agreement is mainly due to the accurate PES reconstruction, as previous studies have shown that the exit point location is fairly robust with respect to variations in the collective inertia [50, 51, 52, 45]. This agreement holds even for nuclei whose PES reconstruction has a large error, such as \({}^{280}\)Cm, indicating that the qualitative features shown in Fig. 2 are reconstructed well enough to describe multimodality in SF.

Notice, however, that the exit point locations are not reproduced perfectly for some nuclei, especially in the thorium (\(Z=90\)) chain, where the difference can be as much as \(5\,\mathrm{b}^{3/2}\). This is not due to the PES reconstruction: Figure 1 shows that the thorium isotopes have RMSE \(\Delta V(Z=90)\lesssim 100\) keV, and the exit point reconstruction when considering configuration (i) is within \(1\,\mathrm{b}^{3/2}\) of the reference value. Additionally, a side-by-side comparison of the collective inertia components does not show a systematic deviation between the reference inertia and the NN reconstruction.

Nevertheless, the error is due to the inaccurate collective inertia reconstruction. However, it is not a systematic error. Rather, random error is present for every deformation considered, and it is the accumulation of this random error that causes the discrepancy. While the location of any individual exit point is not sensitive to the collective inertia, the probability of tunneling to a particular point depends on the probability given in Eq. (1). Because the probability is exponentially dependent on the action (and therefore exponentially dependent on the collective inertia reconstruction), comparatively small errors can add up and actually switch the

Figure 5: The \(Q_{30}\) component of the reconstructed lifetime-weighted exit point, minus the reference \(Q_{30}\) component, in \(\mathrm{b}^{3/2}\). These results were computed using configuration (iii), i.e. the NN was used to reconstruct both the PES and the collective inertia.
dominant exit point, from asymmetric to symmetric and vice versa. This is especially important for nuclei with a wide fission barrier, as the cumulative error along the path is large.

In general, we observe that both the PES and the collective inertia are emulated well enough to predict exit points that agree with the reference data. And, for most nuclei, the dominant mode is also in agreement. Together, this means that the SF fragment yields are in agreement between the reference data and the NN reconstruction for most nuclei under consideration.

### Spontaneous fission half-lives

In this section we examine the performance of the NN when predicting the SF half-life, \(t_{1/2}^{\text{sf}}\). For the sake of simplicity, we do not include triaxiality and pairing correlations as collective degrees of freedom, despite their large impact on the predicted \(t_{1/2}^{\text{sf}}\)[53; 54; 55; 56; 51].

Figure 6 shows \(t_{1/2}^{\text{sf}}\) computed using the reference data vs. \(t_{1/2}^{\text{sf}}\) computed using the NN reconstruction, for configurations (i) and (iii) mentioned above (results for configuration (ii) are similar to those of (iii)). As can be seen, the \(t_{1/2}^{\text{sf}}\) predictions agree well, typically within 3 orders of magnitude across the approximately 80 orders of magnitude under consideration.

Figure 6(a) demonstrates that the PES reconstruction is sufficient to predict \(t_{1/2}^{\text{sf}}\) values that agree well with the reference values. As with the SF fragment yields, this is true even for nuclei with a large \(\Delta V\), e.g. \({}^{280}\)Cm, once again demonstrating that the PES emulation quality is indeed sufficient to reproduce SF observables.

Figure 6(b) includes the collective inertia emulation. As can be seen, the reproduced \(t_{1/2}^{\text{sf}}\) values agree less well with the reference values, although the disagreement is still within 3 orders of magnitude for most nuclei. This is not unexpected: Sec. V.1 shows that the collective inertia emulation, while sufficient for most nuclei, is not accurate enough for all nuclei.

Similar to Sec. V.1, the reason for the disagreement in \(t_{1/2}^{\text{sf}}\) is the accumulation of random errors when the fission pathway goes across the fission barrier. Now, rather than changing the dominant fission mode, \(t_{1/2}^{\text{sf}}\) is simply changed from the reference value in a more-or-less random manner. The effect is most prominent for long-lived nuclei, where errors in the collective inertia add up to a fairly large value as the pathway traverses a wider fission barrier.

While it may be desirable in principle to improve the emulation, the nuclei whose \(t_{1/2}^{\text{sf}}\) values are reproduced with a large error are those predicted to be stable to SF, within the \((Q_{20},Q_{30})\) collective space. As such, errors in the SF observables have little effect on results that are further dependent on \(t_{1/2}^{\text{sf}}\), such as \(r\)-process network calculations.

The inset panels in Fig. 6 magnify the range \(10^{-5}-10^{10}\) s, to highlight the relevant \(r\)-process range. As can be seen, almost all nuclei within this range are reproduced nicely within three orders of magnitude. Therefore, we conclude that NNs are able to reproduce both the PES and the collective inertia well enough that \(t_{1/2}^{\text{sf}}\) is reproduced within 3 orders of magnitude for nuclei for which SF is relevant in the \(r\)-process region.

## VI Conclusions

In this work, we have shown that fully connected feed-forward NNs are able to emulate both the potential energy and the collective inertia across a region of the nuclear chart, in the collective space consisting of the axial quadrupole and octupole moments. In general, the emulation error on the potential energy is about 500 keV, and the largest discrepancies are found in high-energy regions far from the fission path. The inertia tensor is reproduced within roughly an order of magnitude. We find that the NN performance is stable with respect to changes in the architecture, while the rescaling of input variables produces a general improvement overall. Most of the exit points predicted by the NN agree with the DFT predictions within a \((\Delta Q_{20},\Delta Q_{30})=(2\,\text{b},1\,\text{b}^{3/2})\) range. The SF half-lives are usually reproduced within a factor \(10^{3}\) over a span of more than 70 orders of magnitude. We find that the largest source of discrepancies is the emulation of the collective inertia tensor, due to the rapid changes of the inertia tensor in regions where single-particle level crossings are present. For some very long-lived nuclei, the associated error accumulates along the wider fission barrier. Conversely, in nuclei where fission can be a major decay mode, the emulations are in

Figure 6: The half-life predicted using the DFT reference data, \(t_{1/2}^{\text{sf}\text{-DFT}}\), plotted against the half-life computed using the NN reconstruction, \(t_{1/2}^{\text{sf}\text{-NN}}\). Panel (a) shows configuration (i), in which only the PES is emulated; panel (b) shows configuration (iii), in which both the PES and the collective inertia are emulated. The black line marks the diagonal: \(t_{1/2}^{\text{sf}\text{-DFT}}=t_{1/2}^{\text{sf}\text{-NN}}\). Gray bars are drawn at \(t_{1/2}^{\text{sf}\text{-DFT}}\times 10^{\pm 3}\), i.e. 3 orders of magnitude above and below the diagonal. Insets show the range \(10^{-5}-10^{10}\,\text{s}\), to highlight the relevant \(r\)-process range.