# Building Surrogate Models of Nuclear Density Functional Theory with Gaussian
Processes and Autoencoders

Marc Verriere1, Nicolas Schunck1,∗, Irene Kim2,3, Petar Marević1,4,5, Kevin
Quinlan6, Michelle N. NGo6,7, David Regnier8,9, Raphael David Lasseri4

###### Abstract

From the lightest Hydrogen isotopes up to the recently synthesized Oganesson
($Z=118$), it is estimated that as many as about 3000 atomic nuclei could
exist in nature. Most of these nuclei are too short-lived to be occurring on
Earth, but they play an essential role in astrophysical events such as
supernova explosions or neutron star mergers that are presumed to be at the
origin of most heavy elements in the Universe. Understanding the structure,
reactions, and decays of nuclei across the entire chart of nuclides is an
enormous challenge because of the experimental difficulties in measuring
properties of interest in such fleeting objects and the theoretical and
computational issues of simulating strongly-interacting quantum many-body
systems. Nuclear density functional theory (DFT) is a fully microscopic
theoretical framework which has the potential of providing such a
quantitatively accurate description of nuclear properties for every nucleus in
the chart of nuclides. Thanks to high-performance computing facilities, it has
already been successfully applied to predict nuclear masses, global patterns
of radioactive decay like $\beta$ or $\gamma$ decay, and several aspects of
the nuclear fission process such as, e.g., spontaneous fission half-lives.
Yet, predictive simulations of nuclear spectroscopy – the low-lying excited
states and transitions between them – or of nuclear fission, or the
quantification of theoretical uncertainties and their propagation to basic or
applied nuclear science applications, would require several orders of
magnitude more calculations than currently possible. However, most of this
computational effort would be spent into generating a suitable basis of DFT
wavefunctions. Such a task could potentially be considerably accelerated by
borrowing tools from the field of machine learning and artificial
intelligence. In this paper, we review different approaches to applying
supervised and unsupervised learning techniques to nuclear DFT.

##  1 Keywords:

Nuclear density functional theory, Gaussian process, Deep learning,
Autoencoders, RESNET

##  2 Introduction

Predicting all the properties of every atomic nucleus in the nuclear chart,
from Hydrogen all the way to superheavy elements, remains a formidable
challenge. Density functional theory (DFT) offers a compelling framework to do
so, since the computational cost is, in principle, nearly independent of the
mass of the system Eschrig (1996). Because of our incomplete knowledge of
nuclear forces and of the fact that the nucleus is a self-bound system, the
implementation of DFT in nuclei is slightly different from other systems such
as atoms or molecules and is often referred to as the energy density
functional (EDF) formalism Schunck (2019).

Simple single-reference energy density functional (SR-EDF) calculations of
atomic nuclei can often be done on a laptop. However, large-scale SR-EDF
computations of nuclear properties or higher-fidelity simulations based on the
multi-reference (MR-EDF) framework can quickly become very expensive
computationally. Examples where such computational load is needed range from
microscopic fission theory Schunck and Regnier (2022); Schunck and Robledo
(2016) to parameter calibration and uncertainty propagation Kejzlar et al.
(2020); Schunck et al. (2020a) to calculations at the scale of the entire
chart of nuclides Erler et al. (2012); Ney et al. (2020) relevant, e.g., for
astrophysical simulations Mumpower et al. (2016). Many of these applications
would benefit from a reliable emulator of EDF models.

It may be useful to distinguish two classes of quantities that such emulators
should reproduce. What we may call “integral” quantities are quantum-
mechanical observables such as, e.g., the energy, radius, or spin of the
nucleus, or more complex data such as decay or capture rates. By contrast, we
call “differential” quantities the basic degrees of freedom of the theoretical
model. In this article, we focus on the Hartree-Fock-Bogoliubov (HFB) theory,
which is both the cornerstone of the SR-EDF approach and provides the most
common basis of generator states employed in MR-EDF calculations. In the HFB
theory, all the degrees of freedom are encapsulated into three equivalent
quantities: the quasiparticle spinors, as defined either on some spatial grid
or configuration space; the full non-local density matrix
$\rho(\bm{r}\sigma\tau,\bm{r}^{\prime}\sigma^{\prime}\tau^{\prime})$ and
pairing tensor
$\kappa(\bm{r}\sigma\tau,\bm{r}^{\prime}\sigma^{\prime}\tau^{\prime})$, where
$\bm{r}$ refers to spatial coordinates, $\sigma=\pm 1/2$ to the spin
projection and $\tau=\pm 1/2$ to the isopin projection Perlińska et al.
(2004); the full non-local HFB mean-field and pairing potentials, often
denoted by $h(\bm{r}\sigma\tau,\bm{r}^{\prime}\sigma^{\prime}\tau^{\prime})$
and $\Delta(\bm{r}\sigma\tau,\bm{r}^{\prime}\sigma^{\prime}\tau^{\prime})$.

Obviously, integral quantities have the clearest physical meaning and can be
compared to data immediately. For this reason, they have been the focus of
most of the recent efforts in applying techniques of machine learning and
artificial intelligence (ML/AI) to low-energy nuclear theory, with
applications ranging from mass tables Utama et al. (2016); Utama and
Piekarewicz (2017, 2018); Niu and Liang (2018); Neufcourt et al. (2019);
Lovell et al. (2022); Mumpower et al. (2022), $\beta$-decay rates Niu et al.
(2019), or fission product yields Wang et al. (2019); Lovell et al. (2020).
The main limitation of this approach is that it must be repeated for every
observable of interest. In addition, incorporating correlations between such
observables, for example the fact that $\beta$-decay rates are strongly
dependent on $Q_{\beta}$-values which are themselves related to nuclear
masses, is not easy. This is partly because the behavior of observables such
as the total energy or the total spin is often driven by underlying shell
effects that can lead to very rapid variations, e.g. at a single-particle
crossing. Such effects could be very hard to incorporate accurately in a
statistical model of integral quantities.

This problem can in principle be solved by emulating what we called earlier
differential quantities. For example, single-particle crossings might be
predicted reliably with a good statistical model for the single-particle
spinors themselves. In addition, since differential quantities represent, by
definition, all the degrees of freedom of the SR-EDF theory, any observable of
interest can be computed from them, and the correlations between these
observables would be automatically reproduced. In this sense, an emulator of
differential quantities is truly an emulator for the entire SR-EDF approach.
In the much simpler case of the Bohr collective Hamiltonian, such a strategy
gave promising results Lasseri et al. (2020).

The goal of this paper is precisely to explore the feasibility of training
statistical models to learn the degrees of freedom of the HFB theory. We have
explored two approaches: a simple one based on independent, stationary
Gaussian processes and a more advanced one relying on deep neural networks
with autoencoders and convolutional layers.

In Section 3, we briefly summarize the nuclear EDF formalism with Skyrme
functionals with a focus on the HFB theory preserving axial symmetry. Section
4 presents the results obtained with Gaussian processes. After recalling some
general notions about Gaussian processes, we analyze the results of fitting
HFB potential across a two-dimensional potential energy surface in 240Pu.
Section 5 is devoted to autoencoders. We discuss choices made both for the
network architecture and for the training data set. We quantify the
performance of autoencoders in reproducing canonical wavefunctions across a
potential energy surface in 98Zr and analyze the structure of the latent
space.

##  3 Nuclear Density Functional Theory

In very broad terms, the main assumption of density functional theory (DFT)
for quantum many-body systems is that the energy of the system of interest can
be expressed as a functional of the density of particles Parr and Yang (1989);
Dreizler and Gross (1990); Eschrig (1996). Atomic nuclei are a somewhat
special case of DFT, since the nuclear Hamiltonian is not known exactly and
the nucleus is a self-bound system Engel (2007); Barnea (2007). As a result,
the form of the energy density functional (EDF) is often driven by underlying
models of nuclear forces, and the EDF is expressed as a function of non-local,
symmetry-breaking, intrinsic densities Schunck (2019). In the single-reference
EDF (SR-EDF) approach, the many-body nuclear state is approximated by a simple
product state of independent particles or quasiparticles, possibly with some
constraints reflecting the physics of the problem. We note
$\ket{\Phi({\bm{q}})}$ such as state, with $\bm{q}$ representing a set of
constraints. The multi-reference EDF (MR-EDF) approach builds a better
approximation of the exact many-body state by mixing together SR-EDF states.

###  3.1 Energy Functional

The two most basic densities needed to build accurate nuclear EDFs are the
one-body density matrix $\rho$ and the pairing tensor $\kappa$ (and its
complex conjugate $\kappa^{*}$). The total energy of the nucleus is often
written as

$E[\rho,\kappa,\kappa^{*}]=E_{\rm nuc}[\rho]+E_{\rm Cou}[\rho]+E_{\rm
pair}[\rho,\kappa,\kappa^{*}]\,,$ (1)

where $E_{\rm nuc}[\rho]$ represents the particle-hole, or mean-field,
contribution to the total energy from nuclear forces, $E_{\rm Cou}[\rho]$ the
same contribution from the Coulomb force, and $E_{\rm
pair}[\rho,\kappa,\kappa^{*}]$ the particle-particle contribution to the
energy111The pairing contribution lumps together terms coming from nuclear
forces, Coulomb forces and possibly rearrangement terms.. In this work, we
model the nuclear part of the EDF with a Skyrme-like term

$E_{\rm
nuc}[\rho]=\sum_{t=0,1}\int\operatorname{d}^{3}\\!{\bm{r}}\;\chi_{t}(\bm{r})\,,$
(2)

which includes the kinetic energy term and reads generically

$\chi_{t}(\bm{r})=C_{t}^{\rho\rho}\rho_{t}^{2}+C_{t}^{\rho\tau}\rho_{t}\tau_{t}+C_{t}^{JJ}\mathsf{J}^{2}_{t}+C_{t}^{\rho\Delta\rho}\rho_{t}\Delta\rho_{t}+C_{t}^{\rho\nabla
J}\rho_{t}\bm{\nabla}\cdot\bm{J}_{t}\,.$ (3)

In this expression, the index $t$ refers to the isoscalar ($t=0$) or isovector
$(t=1$) channel and the terms $C_{t}^{uu^{\prime}}$ are the coupling constants
associated with the energy functional. The particle density
$\rho_{t}(\bm{r})$, kinetic energy density $\tau_{t}(\bm{r})$, spin-current
tensor $\mathsf{J}_{t}(\bm{r})$, and vector density $\bm{J}_{t}(\bm{r})$ are
all derived from the full one-body, non-local density
$\rho(\bm{r}\sigma\tau,\bm{r}^{\prime}\sigma^{\prime}\tau^{\prime})$ where
$\bm{r}$ are spatial coordinates, $\sigma$ is the intrinsic spin projection,
$\sigma=\pm 1/2$, and $\tau=\pm 1/2$ is the isospin projection; see Engel et
al. (1975); Dobaczewski and Dudek (1996); Bender et al. (2003); Perlińska et
al. (2004); Lesinski et al. (2007) for their actual definition. Since we do
not consider any proton-neutron mixing, all densities are diagonal in isospin
space. The two remaining terms in (1) are treated in exactly the same way as
in Schunck et al. (2020b). In particular, the pairing energy is derived from a
surface-volume density-dependent pairing force

$V^{(\tau)}(\bm{r},\bm{r}^{\prime})=V_{0}^{(\tau)}\left[1-\frac{1}{2}\frac{\rho(\bm{r})}{\rho_{c}}\right]\delta(\bm{r}-\bm{r}^{\prime})\,,$
(4)

where $\rho_{c}=0.16$ fm-3 is the saturation density of nuclear matter.

###  3.2 Hartree-Fock-Bogoliubov Theory

The actual densities in (3) are obtained by solving the Hartree-Fock-
Bogoliubov (HFB) equation, which derives from applying a variational principle
and imposing that the energy be minimal under variations of the densities
Schunck (2019). The HFB equation is most commonly solved in the form of a non-
linear eigenvalue problem. The eigenfunctions define the quasiparticle (q.p.)
spinors. Without proton-neutron mixing, we can treat neutrons and protons
separately. Therefore, for any one type of particles, the HFB equation giving
the $\mu$th eigenstate reads in coordinate space Dobaczewski et al. (1984)

$\int\operatorname{d}^{3}\\!\bm{r}^{\prime}\sum_{\sigma^{\prime}}\left(\begin{array}[]{cc}h(\bm{r}\sigma,\bm{r}^{\prime}\sigma^{\prime})-\lambda\delta_{\sigma\sigma^{\prime}}&\tilde{h}(\bm{r}\sigma,\bm{r}^{\prime}\sigma^{\prime})\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
\tilde{h}^{*}(\bm{r}\sigma,\bm{r}^{\prime}\sigma^{\prime})&-h(\bm{r}\sigma,\bm{r}^{\prime}\sigma^{\prime})+\lambda\delta_{\sigma\sigma^{\prime}}\end{array}\right)\left(\begin{array}[]{c}U(E_{\mu},\bm{r}^{\prime}\sigma^{\prime})\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
V(E_{\mu},\bm{r}^{\prime}\sigma^{\prime})\end{array}\right)=E_{\mu}\left(\begin{array}[]{c}U(E_{\mu},\bm{r}\sigma)\vskip
6.0pt plus 2.0pt minus 2.0pt\\\ V(E_{\mu},\bm{r}\sigma)\end{array}\right),$
(5)

where $h(\bm{r}\sigma,\bm{r}^{\prime}\sigma^{\prime})$ is the mean field,
$\tilde{h}(\bm{r}\sigma,\bm{r}^{\prime}\sigma^{\prime})$ the pairing field222
Following Dobaczewski et al. (1984, 1996), we employ the ‘russian’ convention
where the pairing field is defined from the pairing density
$\tilde{\rho}(\bm{r}\sigma,\bm{r}^{\prime}\sigma^{\prime})$ rather than the
pairing tensor. The quantity $\tilde{h}$ is related to the more traditional
form of the pairing field $\Delta$ through:
$\tilde{h}(\bm{r}\sigma,\bm{r}^{\prime}\sigma^{\prime})=-2\sigma^{\prime}\Delta(\bm{r}\sigma,\bm{r}^{\prime}-\sigma^{\prime})$.
and $\lambda$ the Fermi energy. Such an eigenvalue problem must be solved for
protons and for neutrons.

For the case of Skyrme energy functionals and zero-range pairing functionals,
both the mean field $h$ and pairing field $\tilde{h}$ become semi-local
functions of $\bm{r}$ (semi-local refers to the fact that these potentials
involve differential operators). We refer to Vautherin and Brink (1972); Engel
et al. (1975) for an outline of the derivations leading to the expressions of
the mean field in the case of Skyrme functionals and to, e.g., Dobaczewski and
Dudek (1997); Bender et al. (2009); Hellemans et al. (2012); Ryssens et al.
(2015a) for the expression of the mean field in terms of coupling constants
rather than the parameters of the Skyrme potential. In the following, we
simply recall the essential formulas needed in the rest of the manuscript.

Expression 5 is written in coordinate space. In configuration space, i.e.,
when the q.p. spinors are expanded on a suitable basis of the single-particle
(s.p.) Hilbert space, the same equation becomes a non-linear eigenvalue
problem that can be written as

$\left(\begin{array}[]{cc}h-\lambda&\tilde{h}\\\
\tilde{h}^{*}&-h^{*}+\lambda\end{array}\right)\left(\begin{array}[]{cc}U&V^{*}\\\
V&U^{*}\end{array}\right)=\left(\begin{array}[]{cc}U&V^{*}\\\
V&U^{*}\end{array}\right)\left(\begin{array}[]{cc}-E&0\\\
0&E\end{array}\right),$ (6)

where $h$, $\tilde{h}$, $U$ and $V$ are now $N_{\rm basis}\times N_{\rm
basis}$ matrices, with $N_{\rm basis}$ the number of basis states. Eigenvalues
are collected in the diagonal $N_{\rm basis}\times N_{\rm basis}$ matrix $E$.
The set of all eigenvectors define the Bogoliubov matrix,

$\mathcal{W}=\left(\begin{array}[]{cc}U&V^{*}\\\ V&U^{*}\end{array}\right),$
(7)

which is a unitary:
$\mathcal{W}\mathcal{W}^{\dagger}=\mathcal{W}^{\dagger}\mathcal{W}=1$. Details
about the HFB theory can be found in the standard references Valatin (1961);
Mang (1975); Blaizot and Ripka (1985); Ring and Schuck (2004).

###  3.3 Mean-field and Pairing Potentials

The mean fields are obtained by functional differentiation of the scalar-
isoscalar energy functional (1) with respect to all relevant isoscalar or
isovector densities, $\rho_{0}$, $\rho_{1}$, $\tau_{0}$, etc. For the case of
a standard Skyrme EDF when time-reversal symmetry is conserved, the
corresponding mean-field potentials in the isoscalar-isovector representation
become semi-local Dobaczewski and Dudek (1995, 1997); Stoitsov et al. (2005);
Hellemans et al. (2012)

$h_{t}(\bm{r})=-\bm{\nabla}M_{t}^{*}(\bm{r})\bm{\nabla}+U_{t}(\bm{r})+\frac{1}{2i}\sum_{\mu\nu}\big{(}\nabla_{\mu}\sigma_{\nu}B_{t,\mu\nu}(\bm{r})+B_{t,\mu\nu}(\bm{r})\nabla_{\mu}\sigma_{\nu}\big{)},$
(8)

where, as before, $t=0,1$ refers to the isoscalar or isovector channel and the
various contributions are

$\displaystyle M_{t}(\bm{r})$
$\displaystyle=\frac{\hbar^{2}}{2m}+C_{t}^{\rho\tau}\rho_{t},$ (9a)
$\displaystyle U_{t}(\bm{r})$
$\displaystyle=2C_{t}^{\rho\rho}\rho_{t}+C_{t}^{\rho\tau}\tau_{t}+2C_{t}^{\rho\Delta\rho}\Delta\rho_{t}+C_{t}^{\rho\nabla
J}\bm{\nabla}\cdot\bm{J}_{t}+U_{t}^{(\rm rear)},$ (9b) $\displaystyle
B_{t,\mu\nu}(\bm{r})$ $\displaystyle=2C_{t}^{\rho
J}J_{t,\mu\nu}-C_{t}^{\rho\Delta J}\nabla_{\mu}\rho_{t,\nu}.$ (9c)

In these expressions, $\mu,\nu$ label spatial coordinates and $\bm{\sigma}$ is
the vector of Pauli matrices in the chosen coordinate system. For example, in
Cartesian coordinates, $\mu,\nu\equiv x,y,z$ and
$\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$. The term $U_{t}^{(\rm
rear)}$ is the rearrangement potential originating from the density-dependent
part of the energy. The resulting isoscalar and isovector mean-field and
pairing potentials can then recombined to give the neutron and proton
potentials,

$h^{(n)}=h_{0}+h_{1},\qquad h^{(p)}=h_{0}-h_{1}.$ (10)

Note that the full proton potential should also contain the contribution from
the Coulomb potential.

The pairing field is obtained by functional differentiation of the same energy
functional (1), this time with respect to the pairing density. As a result,
one can show that it is simply given by

$\tilde{h}^{(\tau)}(\bm{r})=V_{0}^{(\tau)}\left[1-\frac{1}{2}\frac{\rho_{0}(\bm{r})}{\rho_{c}}\right]\tilde{\rho}^{(\tau)}(\bm{r}).$
(11)

###  3.4 Collective Space

Nuclear fission or nuclear shape coexistence are two prominent examples of
large-amplitude collective motion of nuclei Schunck and Regnier (2022); Heyde
and Wood (2011). Such phenomena can be accurately described within nuclear DFT
by introducing a small-dimensional collective manifold, e.g., associated with
the nuclear shape, where we assume the nuclear dynamics is confined
Nakatsukasa et al. (2016); Schunck (2019). The generator coordinate method
(GCM) and its time-dependent extension (TDGCM) provide quantum-mechanical
equations of motion for such collective dynamics Griffin and Wheeler (1957);
Wa Wong (1975); Reinhard and Goeke (1987); Bender et al. (2003); Verriere and
Regnier (2020). In the GCM, the HFB solutions are generator states, i.e., they
serve as a basis in which the nuclear many-body state is expanded. The choice
of the collective manifold, that is, of the collective variables, depends on
the problem at hand. For shape coexistence or fission, these variables
typically correspond to the expectation value of multipole moment operators on
the HFB state. A pre-calculated set of HFB states with different values for
the collective variables defines a potential energy surface (PES).

In practice, PES are obtained by adding constraints to the solutions of the
HFB equation. This is achieved by introducing a set of constraining operators
$\hat{Q}_{a}$ capturing the physics of the problem at hand. The set of all
such constraints $\bm{q}\equiv(q_{1},\dots,q_{N})$ defines a point in the PES.
In this work, our goal is to design emulators capable of reproducing the HFB
solutions at any given point $\bm{q}$ of a PES. Throughout this article, we
consider exclusively two-dimensional collective spaces made spanned by the
expectation values of the axial quadrupole $\hat{Q}_{20}$ and axial octupole
$\hat{Q}_{30}$ moment operators. In the presence of constraints, the mean-
field potential in the HFB equation is modified as follows

$h(\bm{r}\sigma,\bm{r}^{\prime}\sigma^{\prime})-\lambda\delta_{\sigma\sigma^{\prime}}\quad\rightarrow\quad
h(\bm{r}\sigma,\bm{r}^{\prime}\sigma^{\prime})-\Big{(}\lambda+\sum_{a}\lambda_{a}Q_{a}(\bm{r})\Big{)}\delta_{\sigma\sigma^{\prime}}.\
$ (12)

As well known, the Fermi energies play in fact the role of the Lagrange
parameters $\lambda_{a}$ for the constraints on particle number. When
performing calculations with constraints on the octupole moment, it is also
important to fix the position of the center of mass. This is typically done by
adding a constraint on the dipole moment $\hat{Q}_{10}$. In the following, we
note $q_{\lambda\mu}$ the expectation value of the operator
$\hat{Q}_{\lambda\mu}$ on the quasiparticle vacuum,
$q_{\lambda\mu}=\langle{\Phi({\bm{q}})}|\hat{Q}_{\lambda\mu}|\Phi({\bm{q}})\rangle$.

Potential energy surfaces are a very important ingredient in a very popular
approximation to the GCM called the Gaussian overlap approximation (GOA) Brink
and Weiguny (1968); Onishi and Une (1975); Une et al. (1976). By assuming,
among other things, that the overlap between two HFB states with different
collective variables $\bm{q}$ and $\bm{q}^{\prime}$ is approximately Gaussian,
the GOA allows turning the integro-differential Hill-Wheeler-Griffin equation
of the GCM into a much more tractable Schrödinger-like equation. The time-
dependent version of this equation reads as Verriere and Regnier (2020)

$i\hbar\frac{\partial}{\partial
t}g(\bm{q},t)=\left[-\frac{\hbar^{2}}{2}\sum_{\alpha\beta}\frac{\partial}{\partial
q_{\alpha}}B_{\alpha\beta}(\bm{q})\frac{\partial}{\partial
q_{\beta}}+V(\bm{q})\right]g(\bm{q},t),$ (13)

where $g(\bm{q},t)$ is the probability to be at point $\bm{q}$ of the
collective space at time $t$, $V(\bm{q})$ is the actual PES, typically the HFB
energy as a function of the collective variables $\bm{q}$ (sometimes
supplemented by some zero-point energy correction) and
$B_{\alpha\beta}(\bm{q})$ the collective inertia tensor. In (13), indices
$\alpha$ and $\beta$ run from 1 to the number $N_{\rm col}$ of collective
variables. While the HFB energy often varies smoothly with respect to the
collective variables, the collective inertia tensor can exhibit very rapid
variations near level crossings.

###  3.5 Canonical Basis

The Bloch-Messiah-Zumino theorem states that the Bogoliubov matrix
$\mathcal{W}$ of (7) can be decomposed into a product of three matrices Ring
and Schuck (2004); Bloch and Messiah (1962); Zumino (1962)

$\mathcal{W}=\mathcal{D}\bar{\mathcal{W}}\mathcal{C}=\left(\begin{array}[]{cc}D&0\\\
0&D^{*}\end{array}\right)\left(\begin{array}[]{cc}\bar{U}&\bar{V}\\\
\bar{V}&\bar{U}\end{array}\right)\left(\begin{array}[]{cc}C&0\\\
0&C^{*}\end{array}\right)\,,$ (14)

where $D$ and $C$ are unitary matrices. The matrices $\bar{U}$ and $\bar{V}$
take the very simple canonical form

$\bar{U}=\left(\begin{array}[]{cccccc}0&&&&\\\ &\ddots&&&\\\ &&u_{k}&0&&\\\
&&0&u_{\bar{k}}&&\\\ &&&\ddots&\\\
&&&&0\end{array}\right),\qquad\bar{V}=\left(\begin{array}[]{cccccc}0&&&&\\\
&\ddots&&&\\\ &&0&v_{k}&&\\\ &&v_{\bar{k}}&0&&\\\ &&&\ddots&\\\
&&&&0\end{array}\right)\,.$ (15)

Starting from an arbitrary s.p. basis $({\hat{c}},{\hat{c}}^{\dagger})$ of the
Hilbert space, the transformation characterized by the matrix $\mathcal{D}$
leads to a new basis $({\hat{a}},{\hat{a}}^{\dagger})$ that diagonalizes the
density matrix $\rho$ and puts the pairing tensor $\kappa$ into the canonical
form similar to that of $\bar{V}$. This new basis is called the canonical
basis of the HFB theory. Properties of the canonical basis are discussed in
details in the literature; see, e.g., Ring and Schuck (2004); Schunck (2019).
In the HFB theory, quasiparticles are superpositions of particle operators
${\hat{a}}^{\dagger}$ and hole operators ${\hat{a}}$. Thus, the canonical
basis is transformed according to the matrix $\bar{\mathcal{W}}$ to obtain a
set of quasiparticle operators $({\hat{\alpha}},{\hat{\alpha}}^{\dagger})$.
There is another transformation of these operators associated with the matrix
$\mathcal{C}$. However, the most important property for the purpose of this
paper is that physical observables associated with HFB solutions do not depend
on that last transformation.

In addition to simplifying the calculation of many-body observables, the
canonical basis is also computationally less expensive than the full
Bogoliubov basis333This statement is obviously not true when solving the HFB
equation directly in coordinate space. In the case of the local density
discussed here, the expression
$\rho(\bm{r})=\sum_{\sigma}\sum_{\mu}V_{\mu}(\bm{r},\sigma)V^{*}_{\mu}(\bm{r},\sigma)$
is just as computationally expensive as the canonical basis expression
$\rho(\bm{r})=\sum_{\sigma}\sum_{\mu}v_{\mu}^{2}|\varphi_{\mu}(\bm{r},\sigma)|^{2}$..
As an illustration, let us take the example of the local density
$\rho(\bm{r})$. Assuming the s.p. basis $({\hat{c}},{\hat{c}}^{\dagger})$ is
represented by the basis functions
$\\{\psi_{n}(\bm{r},\sigma)\\}_{n\in\mathbb{N}}$, the local density (for
isospin $\tau$) is obtained from the matrix of the Bogoliubov transformation
by

$\rho(\bm{r})=\sum_{\sigma}\sum_{\mu}\sum_{mn}V_{m\mu}^{*}V_{n\mu}\,\psi_{m}(\bm{r},\sigma)\psi^{*}_{n}(\bm{r},\sigma).$
(16)

Notwithstanding the constraints imposed by the orthonormality of the q.p.
spinors, the number of independent parameters in this expression approximately
scales like $2\times N_{\rm basis}^{2}\times N_{\rm qp}\times N_{\mathbf{r}}$,
where $N_{\rm basis}$ is the size of the s.p. basis, $N_{\rm qp}$ the number
of q.p. states $\mu$ and $N_{\mathbf{r}}$ the total number of points in the
spatial grid $\bm{r}$ (which depends on the symmetries imposed). In the
canonical basis, and assuming that the state
${\hat{a}}_{\mu}^{\dagger}\ket{0}$ is associated with the wavefunction
$\varphi_{\mu}(\bm{r},\sigma)$, the same object is represented by

$\rho(\bm{r})=\sum_{\sigma}\sum_{\mu}v_{\mu}^{2}|\varphi_{\mu}(\bm{r},\sigma)|^{2}.$
(17)

The number of data points now scales like $2\times N_{\rm qp}\times
N_{\mathbf{r}}+N_{\rm qp}$, or about $N_{\rm basis}^{2}$ smaller than before.
For calculations with $N_{\rm basis}\approx 1000$ the compression enabled by
the canonical basis is of the order of $10^{6}$.

###  3.6 Harmonic Oscillator Basis

All calculations in this article were performed with the HFBTHO code Marević
et al. (2022). Recall that HFBTHO works by expanding the solutions on the
axially-deformed harmonic oscillator basis Stoitsov et al. (2005).
Specifically, the HO basis functions are written

$\psi_{\bm{n}}(\bm{r},\sigma)=\psi_{n_{r}}^{\Lambda}(r)\psi_{n_{z}}(z)\frac{e^{i\Lambda\theta}}{\sqrt{2\pi}}\chi_{\Sigma}(\sigma),$
(18)

where $\bm{n}\equiv(n_{r},n_{z},\Lambda,\Omega=\Lambda\pm\Sigma)$ are the
quantum numbers labeling basis states and

$\displaystyle\psi_{n_{r}}^{\Lambda}(r)$
$\displaystyle=N_{n_{r}}\beta_{\perp}\sqrt{2}\eta^{|\Lambda|/2}e^{-\eta/2}L_{n_{r}}^{|\Lambda|}(\eta),$
(19a) $\displaystyle\psi_{n_{z}}(z)$
$\displaystyle=N_{n_{z}}\sqrt{\beta_{z}}e^{-\xi^{2}/2}H_{n_{z}}(\xi),$ (19b)

with $\eta=\beta_{\perp}^{2}r^{2}$ and $\xi=\beta_{z}z$ dimensionless
variables, $L_{n_{r}}^{|\Lambda|}$ the associated Laguerre polynomials of
order $n_{r}$ and $H_{n_{z}}$ the Hermite polynomial of order $n_{z}$. The
oscillator scaling factors $\beta_{\perp}$ and $\beta_{z}$ are the inverse of
the oscillator lengths, i.e., $\beta_{z}=1/b_{z}$.

All integrations are performed by Gauss quadrature, namely Gauss-Hermite for
integrations along the $\xi$-axis of the intrinsic reference frame and Gauss-
Laguerre for integrations along the perpendicular direction characterized by
the variable $\eta$. In the following, we note $N_{z}$ the number of Gauss-
Hermite nodes and $N_{\perp}$ the number of Gauss-Laguerre nodes.

##  4 Supervised Learning with Gaussian Processes

Gaussian processes (GPs) are a simple yet versatile tool for regression that
has found many applications in low-energy nuclear theory over the past few
years, from determining the nuclear equation of state Drischler et al. (2020)
to nuclear cross sections calculations Kravvaris et al. (2020); Acharya and
Bacca (2022) to modeling of neutron stars Pastore et al. (2017). In the
context of nuclear DFT, they were applied to build emulators of $\chi_{2}$
objective functions in the UNEDF project Kortelainen et al. (2010, 2012,
2014); Higdon et al. (2015); McDonnell et al. (2015); Schunck et al. (2020a),
of nuclear mass models Neufcourt et al. (2018, 2019, 2020a, 2020b) or of
potential energy surfaces in actinides Schunck et al. (2020b). In this
section, we test the ability of GPs to learn directly the HFB potentials
across a large, two-dimensional collective space.

###  4.1 Gaussian Processes

Gaussian processes are commonly thought of as the generalization of normally-
distributed random variables (Gaussian distribution) to functions. There
exists a considerable field of applications for GPs and we refer to the
reference textbook by Rasmussen & Williams for a comprehensive review of the
formalism and applications of GPs Rasmussen and Williams (2006). For the
purpose of this work, we are only interested in the ability of GPs to be used
as a regression analysis tool and we very briefly outline below some of the
basic assumptions and formulas.

We assume that we have a dataset of observations
$\\{\mathbf{y}=y_{i}\\}_{i=1,\dots,n}$ and that these data represent $n$
realizations of

$y=f(\mathbf{x})+\epsilon,$ (20)

where $f:\mathbf{x}\mapsto f(\mathbf{x})$ is the unknown function we are
seeking to learn from the data. Saying that a function $f$ is a Gaussian
process means that every finite collection of function values
$\mathbf{f}=(f(\mathbf{x}_{1}),\dots,f(\mathbf{x}_{p}))$ follows a
$p$-dimensional multivariate normal distribution. In other words, we assume
that the unknown function $f$ follows a normal distribution in ‘function
space’. This is denoted by

$f(\mathbf{x})\sim\mathcal{GP}\big{(}m(\mathbf{x}),k(\mathbf{x},\mathbf{x}^{\prime})\big{)},$
(21)

where $m:\mathbf{x}\mapsto m(\mathbf{x})$ is the mean function and
$k:(\mathbf{x},\mathbf{x}^{\prime})\mapsto k(\mathbf{x},\mathbf{x}^{\prime})$
the covariance function, which is nothing but the generalization to functions
of the standard deviation,

$k(\mathbf{x},\mathbf{x}^{\prime})=\mathbb{E}\left[\big{(}f(\mathbf{x})-m(\mathbf{x})\big{)}\big{(}f(\mathbf{x}^{\prime})-m(\mathbf{x}^{\prime})\big{)}\right].$
(22)

Thanks to the properties of Gaussian functions, the mean and covariance
functions have analytical expressions as a function of the test data
$\mathbf{y}$ and covariance $k$; see Eqs.(2.25)-(2.26) in Rasmussen and
Williams (2006).

The covariance function is the central object in GP regression. It is
typically parametrized both with a functional form and with a set of free
parameters called hyperparameters. The hyperparameters are determined from the
observed data by maximizing the likelihood function. In our tests, the
covariance matrix is described by a standard Matérn 5/2 kernel,

$k(\mathbf{x},\mathbf{x}^{\prime})=\left(1+\frac{\sqrt{5}}{\ell}||\mathbf{x}-\mathbf{x}^{\prime}||+\frac{5}{3\ell^{2}}||\mathbf{x}-\mathbf{x}^{\prime}||^{2}\right)\exp\left(-\frac{\sqrt{5}}{\ell}||\mathbf{x}-\mathbf{x}^{\prime}||\right)\,,$
(23)

where $\ell$ is the length-scale that characterizes correlations between
values of the data at different locations. The length-scale is a hyper-
parameter that is optimized in the training phase of the Gaussian process. In
this work, we only considered stationary GPs: the correlation between data
points $\mathbf{x}$ and $\mathbf{x}^{\prime}$ only depends on the distance
$||\mathbf{x}-\mathbf{x}^{\prime}||$ between these points, not on their actual
value. This approximation may be too restrictive.

###  4.2 Study Case

####  4.2.1 HFB Potentials

Section 3.2 showed that the HFB mean-field potential involves several
differential operators. When the HFB matrix is constructed by computing
expectation values of the HFB potential on basis functions, differentiation is
carried over to the basis functions and computed analytically – one of the
many advantages of working with the HO basis. In practice, this means that the
elements of the HFB matrix are computed by multiplying spatial kernels with
different objects representing either the original HO functions or their
derivatives. This means that we cannot consider a single emulator for the
entire HFB potential. Instead, we have to build several different ones for
each of its components: the central potential $U$ (derivative of the EDF with
respect to $\rho$), the $r$\- and $z$-derivatives of the effective mass
$M^{*}$ (derivative with respect to the kinetic density $\tau$), the $r$\- and
$z$-derivatives of the spin-orbit potential $W$, and the pairing field
$\tilde{h}$. There are six such functions for neutrons and another six for
protons. We denote this set of twelve functions as
$\\{f_{i}\\}_{i=1,\dots,12}$.

At any given point $\bm{q}$ of the collective space, these functions are all
local, scalar functions of $\eta$ and $\xi$, $f_{i}(\bm{q})\equiv
f_{i}:(\eta,\xi)\mapsto f_{i}(\eta,\xi;\bm{q})$ where $(\eta,\xi)$ are the
nodes of the Gauss-Laguerre and Gauss-Hermite quadrature grid. We note
generically $f_{ik}(\bm{q})$ the value at point $k$ of the quadrature grid
(linearized) of the sample at point $\bm{q}$ of the function $f_{i}$. When
fitting Gaussian process to reproduce mean-field and pairing potentials, we
consider a quadrature grid of $N_{z}\times N_{\perp}=3200$ points. Our goal is
thus to build 3200 different emulators, one for each point $k$ of that grid,
for each of the 12 local functions characterizing the mean-field and pairing
potentials. This gives a grand total of 38 400 emulators to build. While this
number is large, it is still easily manageable on standard computers. It is
also several orders of magnitude smaller than emulating the full set of
quasiparticle spinors, as we will see in the next section.

In addition, the value of all the Lagrange parameters used to set the
constraints must also be included in the list of data points. In our case, we
have $5$ of them: the two Fermi energies $\lambda_{n}$ and $\lambda_{p}$ and
the three constraints on the value of the dipole, quadrupole and octupole
moments, $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$, respectively.
Finally, we also fit the expectation value of the three constraints on
$\hat{Q}_{10}$, $\hat{Q}_{20}$ and $\hat{Q}_{30}$. We thus have a grand total
of 38 408 functions of $\bm{q}$ to emulate.

Figure 1: Potential energy surface of 240Pu with the SkM* EDF for the grid
$(q_{20},q_{30})\in[0\,{\rm b},300\,{\rm b}]\times[0\,{\rm b}^{3/2},51\,{\rm
b}^{3/2}]$ with steps $\delta q_{20}=6\,{\rm b}$ and $\delta q_{30}=3\,{\rm
b}^{3/2}$. The black crosses are the training points, the white circles the
validation points. Energies indicated by the color bar are in MeV relatively
to -1820 MeV.

####  4.2.2 Training Data and Fitting Procedure

We show in Fig. 1 the potential energy surface that we are trying to
reconstruct. This PES is for the 240Pu nucleus and was generated with the SkM*
parameterization of the Skyrme energy functional Bartel et al. (1982). The
pairing channel is described with the zero-range, density-dependent pairing
force of Eq. (4) that has exactly the same characteristics as in Schunck et
al. (2014).

We imposed constraints on the axial quadrupole and octupole moments such that:
$0\,{\rm b}\leq q_{20}\leq 300\,{\rm b}$ and $0\,{\rm b}^{3/2}\leq q_{30}\leq
51\,{\rm b}^{3/2}$ with steps of $\delta q_{20}=6\,{\rm b}$ and $\delta
q_{30}=3\,{\rm b}^{3/2}$, respectively. The full PES should thus contain 918
collective points. In practice, we obtained $N_{\rm p}=887$ fully converged
solutions. Calculations were performed with the HFBTHO solver by expanding the
solutions on the harmonic oscillator basis with $N_{\rm max}=28$ deformed
shells and a truncation in the number of states of $N_{\rm basis}=1000$. At
each point of the PES, the frequency $\omega_{0}$ and deformation $\beta_{2}$
of the HO basis are set according to the empirical formulas given in Schunck
et al. (2014). Following standard practice, we divided the full $N_{\rm
p}=887$ dataset of points into a training (80% of the points) and validation
(20% of the points) set. The selection was done randomly and resulted in
$N_{\rm train}=709$ training points and $N_{\rm valid.}=178$ validation
points. The training points are marked as small black crosses in Fig. 1 while
the validation points are marked as larger white circles.

Based on the discussion in Section 4.2.1, we fit a Gaussian process to each of
the 38 408 variables needed to characterize completely the HFB matrix. Since
we work in a two-dimensional collective space, we have two features and the
training data is represented by a two-dimensional array $X$ of dimension
$(n_{\rm samples},n_{\rm features})$ with $n_{\rm samples}=N_{\rm p}$ and
$n_{\rm features}=2$. The target values $Y$ (= the value at point $k$ on the
quadrature grid of any of the functions $f_{i}$) are contained in a one-
dimensional array of size $N_{\rm p}$. Prior to the fit, the data is
normalized between 0 and 1. The GP is based on a standard Matérn kernel with
$\nu=2.5$ and length-scale $\ell$. In practice, we use different length-scales
for the $q_{20}$ and $q_{30}$ directions so that $\ell=\bm{\ell}$ is a vector.
We initialized these values at the spacing of the grid $\bm{\ell}=(\delta
q_{20},\delta q_{30})$. We added a small amount of white noise to the Matérn
kernel to account for the global noise level of the data.

####  4.2.3 Performance

Once the GP has been fitted on the training data, we can estimate its
performance on the validation data. For each of the $N_{\rm valid.}=178$
validation points, we used the GP-fitted HFB potentials to perform a single
iteration of the HFB self-consistent loop and extract various observables from
this single iteration. Figure 2 focuses on the total HFB energy and the zero-
point energy correction $\varepsilon_{0}$. Together, these two quantities
define the collective potential energy in the collective Hamiltonian (13) of
the GCM. The left panel of the figure shows the histogram of the error $\Delta
E=E_{\rm HFB}^{\rm(true)}-E_{\rm HFB}^{\rm(GP)}$, where $E_{\rm
HFB}^{\rm(true)}$ is the result from the fully converged HFB solution and
$E_{\rm HFB}^{\rm(GP)}$ is the value predicted by the Gaussian process. The
bin size is 100 keV. Overall, we find that the large majority of the error is
within $\pm 200$ keV. This is a rather good result considering the span of the
PES and the fact that basis truncation errors can easily amount to a few MeV
Schunck (2013a).

Figure 2: Left: Histogram of the error on the GP-predicted total HFB energy
and zero-point energy correction across the validation points. Bin size is 100
keV. Right: Size of the error on the GP-predicted total HFB energy across the
validation set. Gray circles have an error lower than 500 keV and the size of
the markers correspond to energy bins of 100 keV. Black circles have an error
greater than 500 keV and are binned by 400 keV units. Energies indicated by
the color bar are in MeV relatively to -1820 MeV.

To gain additional insight, we draw in the right panel of Fig. 2 each of the
validation points with a marker, the size of which is proportional to the
error of the prediction. To further distinguish between most points and the
few outliers, we show in gray the points for which the absolute value of the
error is less than 500 keV and in black the points for which it is greater
than 500 keV. For the gray points, we use 5 different marker sizes
corresponding to energy bins of 100 keV: the smaller grey symbol corresponds
to an error smaller than 100 keV, the larger one between 400 and 500 keV.
Similarly, the larger black circles have all an error greater than 500 keV and
are ordered by bins of 400 keV (there are only two points for which the error
is larger than 4 MeV). Interestingly, most of the larger errors are
concentrated in the region of small elongation $q_{20}<80\,{\rm b}$ and high
asymmetry $q_{30}>30\,{\rm b}^{3/2}$. This region of the collective space is
very high in energy (more than 100 MeV above the ground state) and plays no
role in the collective dynamics.

Note that the expectation values of the multipole moments themselves are not
reproduced exactly by the GP: strictly speaking, the contour plot in the right
panel of Fig. 2 is drawn based on the requested values of the constraints, not
their actual values as obtained by solving the HFB equation once with the
reconstructed potentials. The histogram in the left panel of Fig. 3 quantifies
this discrepancy. It shows the absolute error $\Delta
q_{\lambda\mu}=q_{\lambda\mu}^{(\rm true)}-q_{\lambda\mu}^{(\rm GP)}$, where
$q_{\lambda\mu}^{(\rm true)}$ is the result from the fully converged HFB
solution and $q_{\lambda\mu}^{(\rm true)}$ is the value predicted by the
Gaussian process. On average, the error remains within $\pm 0.5\,\mathrm{b}$
for $q_{20}$ and $\pm 0.5\,\mathrm{b}^{3/2}$ for $q_{30}$, which is
significantly smaller than the mesh size.

Figure 3: Left: Histogram of the error on the GP-predicted values of the
multipole moments. The bin size is 0.2 $b^{\lambda/2}$ with $\lambda=2$
(quadrupole moment) or $\lambda=3$ (octupole moment). Right: Histogram of the
relative error, in percents, on the GP-predicted values of the components of
the collective inertia tensor. The bin size is 1, corresponding to 1% relative
errors.

The collective potential energy is only one of the two ingredients used to
simulate fission dynamics. As mentioned in Sec. 3.4, see Eq. (13), the
collective inertia tensor is another essential quantity Schunck and Robledo
(2016); Schunck and Regnier (2022). In this work, we computed the collective
inertia at the perturbative cranking approximation Schunck and Robledo (2016).
Since we work in two-dimensional collective spaces, the collective inertia
tensor $\mathsf{B}$ has three independent components, hereafter labeled
$B_{22}$, $B_{33}$ and $B_{32}=B_{23}$. Figure 3 shows the relative error on
these quantities, defined as $\epsilon=(B_{\lambda\lambda^{\prime}}^{(\rm
true)}-B_{\lambda\lambda^{\prime}}^{(\rm
GP)})/B_{\lambda\lambda^{\prime}}^{(\rm true)}$. Overall, the error is more
spread than for the energy but rarely exceeds five percents444Note that
$B_{32}$ vanishes for axially-symmetric shapes. As a result, the relative
error can be artificially large for values of $q_{30}\approx
0\,\mathrm{b}^{3/2}$..

Both the total energy and the collective mass tensor are computed from the HFB
solutions. However, since the GP fit is performed directly on the mean-field
and pairing potentials, one can analyze the error on these quantities
directly. In Fig. 4, we consider two different configurations. The
configuration
$\mathcal{C}_{1}=(q_{20},q_{30})=(198\,\mathrm{b},30\,\mathrm{b}^{3/2})$ is
very well reproduced by the GP with an error in the HFB energy of 4.4 keV and
a relative error on $B_{22}$ of -0.43 % and $B_{22}$ of -0.84 % only. In
contrast, the configuration
$\mathcal{C}_{2}=(q_{20},q_{30})=(138\,\mathrm{b},51\,\mathrm{b}^{3/2})$ is
one of the worst possible cases, with a total error on the HFB energy of 9.0
MeV and relative errors on $B_{22}$ of -71.0 % and $B_{22}$ of -13.7 %. For
each of these two configurations, we look at the central part of the mean-
field potential for protons, the term $U_{p}=U_{0}-U_{1}$ of (9b). The left
side of Fig. 4 shows, respectively, the actual value of $U_{p}(r,z)$ across
the quadrature grid (top panel) and the difference between the true value and
the GP fit (bottom panel) for the configuration $\mathcal{C}_{1}$. The right
side of the figure shows the same quantity for the configuration
$\mathcal{C}_{2}$. In all four plots, the energy scale is in MeV.

Figure 4: Upper left: Central part of the mean-field potential for protons,
$U_{p}(r,z)$ for the configuration
$(q_{20},q_{30})=(198\,\mathrm{b},30\,\mathrm{b}^{3/2})$; bottom left: Error
in the GP fit for that same configuration. Upper right: Central part of the
mean-field potential for protons, $U_{p}(r,z)$ for the configuration
$(q_{20},q_{30})=(138\,\mathrm{b},51\,\mathrm{b}^{3/2})$; bottom right: Error
in the GP fit for that same configuration. For all figures, $i_{\rm GH}$ and
$j_{\rm GL}$ refer to the index $i$ and $j$ on the Gauss-Hermite and Gauss-
Laguerre quadrature grid, and the energy given by the error bar is in MeV.

We see that for the ‘good’ configuration $\mathcal{C}_{1}$, the error is
between -0.6 MeV and 1.0 MeV but is mostly occurring at the surface of the
nucleus and at the edges of the domain. Conversely, the ‘bad’ configuration
$\mathcal{C}_{2}$ actually corresponds to a scissioned configuration: the
mean-field potential (upper right panel) shows two different regions
corresponding to fully separated fragments555This particular scission
configuration corresponds to what is called cluster radioactivity Warda and
Robledo (2011); Warda et al. (2018); Matheson et al. (2019). The heavy
fragment is much larger than the light one. Here, $\langle A_{\rm
H}\rangle=205.6$, $\langle A_{\rm L}\rangle=34.4$. Such a geometric
configuration is very different from the rest of the potential energy surface
shown in Fig. 1, which contains mostly non-scissioned configurations. As a
result, the error in the GP fit is very large in the region between the two
fragments since it predicts this configuration to be non-scissioned. Note that
in HFBTHO, the representation of the potentials on the quadrature points does
not contain the exponential factor
$\exp(-\beta_{z}\xi^{2})\exp(-\beta^{2}_{\perp}\xi^{2})$ which is factored in
in the quadrature weights. Therefore, the large errors at the edges of the
domain, for $i_{\rm GH}\approx 0$, $i_{\rm GH}\approx 80$ or $j_{\rm
GL}\approx 40$ are not significant since they are entirely absorbed by this
exponential factor.

Overall, Gaussian processes seem to provide an efficient way to predict HFB
solutions across potential energy surfaces. Their primary advantage is that
they are very simple to implement, with several popular programming
environments offering ready-to-use, full GP packages, and very fast to train
(a few minutes at most for a few hundreds of samples). As our examples
suggest, GPs are very good at interpolating across a domain where solutions
behave smoothly. In the case of PES, this implies that the training data must
not contain, e.g., scissioned and non-scissioned configurations. More
generally, it should not feature too many discontinuities Dubray and Regnier
(2012). When these conditions are met, GPs can be used to quickly and
precisely densify a PES, e.g., to obtain more precise fission paths in
spontaneous fission half-live calculations Sadhukhan (2020).

However, Gaussian processes are intrinsically limited. In our example, we
treated the value of each potential at each point of the quadrature mesh as an
independent GP. Yet, such data are in reality heavily correlated. To
incorporate such correlations requires generalizing from scalar GPs to vector,
or multi-output GPs Bruinsma et al. (2020). In our example of nuclear
potentials, the output space would be $\mathbb{R}^{D}$ with $D\approx 32008$.
An additional difficulty is related to choosing the kernel that is appropriate
to describe the correlated data and identifying what the prior distribution
should be Álvarez et al. (2012). Yet another deficiency of standard Gaussian
processes, especially in contrast to the deep-learning techniques discussed
below, is that they are not capable to learn a latent representation of the
data. For these reasons, we consider such techniques helpful mostly to densify
existing potential energy surfaces.

##  5 Deep Learning with Autoencoders

Even though self-consistent potential energy surfaces are key ingredients in
the microscopic theory of nuclear fission Bender et al. (2020), we must
overcome two significant obstacles to generate reliable and complete PES.
First, the computational cost of nuclear DFT limits the actual number of
single-particle d.o.f.. When solving the HFB equation with basis-expansion
methods, for example, the basis must be truncated (up to a maximum of about a
few thousand states, typically), making the results strongly basis-dependent
Schunck (2013b); even in mesh-based methods, the size of the box and lattice
spacing also induce truncation effects Ryssens et al. (2015b); Jin et al.
(2017). Most importantly, the number of collective variables that can be
included in the PES is also limited: in spontaneous fission calculations,
which do not require a description of the PES up to scission, up to $N_{\rm
col}=5$ collective variables have been incorporated Sadhukhan (2020); when
simulating the PES up to scission, only 2 collective variables are included
with only rare attempts to go beyondRegnier et al. (2017); Zhao et al. (2021).
As a consequence, the combination of heavily-truncated collective spaces and
the adiabatic hypothesis inherent to such approaches leads to missing regions
in the PES and spurious connections between distinct channels with unknown
effects on physics predictions Dubray and Regnier (2012); Lau et al. (2022).
The field of deep learning may offer an appealing solution to this problem by
allowing the construction of low-dimensional and continuous surrogate
representations of potential energy surfaces. In the following, we test the
ability of autoencoders – a particular class of deep neural networks – to
generate accurate low-dimensional representations of HFB solutions.

###  5.1 Network Architecture

The term ‘deep learning’ encompasses many different types of mathematical and
computational techniques that are almost always tailored to specific
applications. In this section, we discuss some of the specific features of the
data we seek to encode in a low-dimensional representation, which in turn help
constrain the network architecture. The definition of a proper loss function
adapted to quantum-mechanical datasets is especially important.

####  5.1.1 Canonical States

We aim at building a surrogate model for determining canonical wavefunctions
as a function of a set of continuous constraints. Canonical states are denoted
generically $\varphi^{(\tau)}_{\mu}(\bm{r},\sigma)$ with
$\bm{r}\equiv(r,z,\theta)$ the cylindrical coordinates and $\sigma=\pm 1/2$
the spin. Fully characterizing an HFB state requires the set of canonical
wavefunctions for both neutrons and protons, which are distinguished by their
isospin quantum number $\tau=+1/2$ (neutrons) and $\tau=-1/2$ (protons). As
mentioned in Sec. 3, an HFB solution ${\ket{\Phi({\bm{q}})}}$ is entirely
determined up to a global phase by the set of all canonical states
${\\{\varphi^{(\tau)}_{\mu}(\bm{r},\sigma)\\}_{\mu}}$ and their associated
occupation amplitudes ${\\{v^{(\tau)}_{\mu}\\}_{\mu}}$.

In this work, we restrict ourselves to axially-symmetric configurations. In
that case, the canonical wavefunctions are eigenstates of the projection of
the total angular momentum on the symmetry axis $\hat{j}_{\rm z}$ with
eigenvalue $\Omega$ and acquire the same separable structure (18) as the HO
basis functions,

$\varphi^{(\tau)}_{\mu}(\bm{r},\sigma)=\varphi^{(\tau)}_{\mu}(r,z,\sigma)\frac{e^{i\Lambda\theta}}{\sqrt{2\pi}},$
(24)

where $\varphi^{(\tau)}_{\mu}(r,z,\sigma)$ is the canonical wavefunction at
$\theta=0$. In this initial work, we only consider even-even nuclear systems
and time-reversal symmetric nuclear Hamiltonians. Therefore, Kramer’s
degeneracy ensures that paired particles in the canonical basis are time-
reversal partners of each other:
$\varphi^{(\tau)}_{\bar{\mu}}(\bm{r},\sigma)=2\sigma\varphi^{(\tau)*}_{\mu}(\bm{r},-\sigma)$.
This guarantees that the canonical wavefunction at $\theta=0$ can be chosen
purely real. Incidentally, it also means that we only need to describe one
wavefunction per pair of particles. Using these properties, we can completely
describe a canonical wavefunction in our model by only predicting a single
pair of real-valued functions (one for each spin projection $\sigma$).

As shown by (8), (11) and (9a)-(9c), all mean-field and pairing potentials are
functions of the Skyrme densities. The kinetic energy density $\tau(r,z)$,
spin-current tensor $\mathsf{J}(r,z)$, and vector density $\bm{J}(r,z)$
involve derivatives of the quasiparticle spinors or, in the canonical basis,
of the canonical wavefunctions on the quadrature grid Stoitsov et al. (2005).
We compute these derivatives by first extracting the coefficients
$\alpha^{(\tau)}_{\bm{n}\mu}$ of the expansion of the canonical wavefunctions
$\varphi^{(\tau)}_{\mu}(\bm{r},\sigma)$ in the HO basis

$\varphi^{(\tau)}_{\mu}(\bm{r},\sigma)=\sum_{\bm{n}}\alpha^{(\tau)}_{\bm{n}\mu}\psi_{\bm{n}}(\bm{r},\sigma)\Rightarrow\alpha^{(\tau)}_{\bm{n}\mu}=\int\operatorname{d}^{3}\\!\bm{r}\,\psi^{*}_{\bm{n}}(\bm{r},\sigma)\varphi^{(\tau)}_{\mu}(\bm{r},\sigma),$
(25)

using Gauss-Laguerre and Gauss-Hermite quadrature. Since all the derivatives
of the HO functions can be computed analytically, the expansion (25) makes it
very easy to compute partial derivatives with respect to $r$ or $z$, for
example,

$\frac{\partial\varphi^{(\tau)}_{\mu}}{\partial
z}(\bm{r},\sigma)=\sum_{\bm{n}}\alpha^{(\tau)}_{\bm{n}\mu}\frac{\partial\psi_{\bm{n}}}{\partial
z}(\bm{r},\sigma).$ (26)

####  5.1.2 Structure of the Predicted Quantity

In the ideal case, the canonical wavefunctions evolve smoothly with the
collective variables. The resulting continuity of the many-body state with
respect to collective variables is a prerequisite for a rigorous description
of the time evolution of fissioning systems, yet it is rarely satisfied in
practical calculations. We discuss below the three possible sources of
discontinuity of the canonical wavefunctions in potential energy surfaces.

First, the canonical wavefunctions are invariant through a global phase. Since
the quantity we want to predict is real, the orbitals can be independently
multiplied by an arbitrary sign. Even though this type of discontinuity does
not impact the evolution of global observables as a function of deformation,
it affects the learning of the model: since we want to obtain continuous
functions, a flipping of the sign would be seen by the neural network as a
discontinuity in the input data. We address this point through the choice of
the loss, as discussed in Sec. 5.1.3, and through the determination of the
training set, as detailed in Sec. 5.2.

Second, we work within the adiabatic approximation, which consists in building
PES by selecting the q.p. vacuum that minimizes the energy at each point. When
the number $N_{\rm col}$ of collective variables of the PES is small, this
approximation may lead to discontinuities Dubray and Regnier (2012). These
discontinuities correspond to missing regions of the collective space and are
related to the inadequate choice of collective variables. Since we want to
obtain a continuous description of the fission path, we must give our neural
network the ability to choose the relevant degrees of freedom. This could be
achieved with autoencoders. Autoencoders are a type of neural networks
analogous to the zip/unzip programs. They are widely used and greatly
successful for representation learning – the field of Machine Learning that
attempts to find a more meaningful representation of complex data Baldi
(2012); Burda et al. (2015); Chen et al. (2012); Gong et al. (2019); Bengio et
al. (2013); Zhang et al. (2014); Yu et al. (2017) and can be viewed as a non-
linear generalization of principal component analysis (PCA). As illustrated in
Fig. 5, an autoencoder $\Xi$ typically consists of two components. The encoder
$E(T^{({\varphi})})$ encodes complex and/or high-dimensional data
$T^{({\varphi})}$ to a typically lower-dimensional representation
$\bm{v}^{({\varphi})}$. The latent space is the set of all possible such
representations. The decoder $D(\bm{v}^{({\varphi})})$ takes the low-
dimensional representation of the encoder and uncompresses it into a tensor
$T^{({\phi})}$ as close as possible to $T^{({\varphi})}$. Such architectures
are trained using a loss function that quantifies the discrepancy between the
initial input and the reconstructed output,

$\mathcal{L}_{\rm rec.}(T^{({\phi})})=d(T^{({\varphi})},T^{({\phi})}),$ (27)

where $d(.,.)$ defines the metric in the space of input data. We discuss the
choice of a proper loss in more details in Sec. 5.1.3.

Figure 5: An autoencoder is the association of two blocks. The first one, on
the left, compresses the input data into a lower-dimensional representation,
or code, in the latent space. The second one, on the right, decompresses the
code back into the original input.

Third, the evolution of the q.p. wavefunctions as a function of the collective
variables $\bm{q}$ may lead to specific values $\bm{q}_{i}$ where the q.p.
solutions are degenerate. These degeneracies form a sub-manifold of dimension
at most $D-2$, where $D$ is the number of collective d.o.f.s. As a
consequence, they cannot appear in one-dimensional PES: q.p. solutions with
the same symmetry “cannot cross” (the famous no-crossing rule von Neuman and
Wigner (1929)). In multi-dimensional spaces, this rule does not hold anymore:
when following a closed-loop trajectory around such a degeneracy, the sign of
the q.p. wavefunctions is flipped, in a similar manner that we flip side when
winding around a Moebius strip Teller (1937); Longuet-Higgins et al. (1958);
Longuet-Higgins (1975). In the field of quantum chemistry, such degeneracies
are referred to as diabolical points or conical intersections Domcke et al.
(2011); Larson et al. (2020). The practical consequence of conical
intersections for deep learning is that the manifold of all the q.p.
wavefunctions cannot be embedded in a $D$-dimensional latent space. Such
singularities can be treated in two ways: (i) by using a latent space of
higher dimension than needed or (ii) by implementing specific neural network
layers capable of handling such cases. For now, we do not consider these
situations.

####  5.1.3 Loss Functions and Metrics

As already discussed in Sec. 5.1.2, autoencoders are trained through the
minimization of a loss function that contains a reconstruction term of the
form (27). As suggested by its name, this term ensures that the autoencoder
can correctly reconstruct the input tensor $T^{({\varphi})}$ from its
compressed representation. It depends on a definition for the metric $d(.,.)$
used to compare the different elements of the input space. Since our canonical
wavefunctions $\varphi_{\mu}$ are expanded on the axial harmonic oscillator
basis of Sec. 3.6, they are discretized on the Gauss quadrature mesh without
any loss of information. Therefore, both the input and output tensors of our
surrogate model are a rank-3 tensor $T^{({\phi})}\equiv T^{({\varphi})}\equiv
T_{ijk}$ of dimensions $N_{\perp}\times N_{z}\times 2$, where $i$ is the index
of the Gauss-Laguerre node along the $r$-axis, $j$ the index of the Gauss-
Hermite node along the $z$-axis, and $k$ the index of the spin component.

A standard loss used with autoencoders is the mean-square-error (MSE). Because
of the structure of our input data, see Sec. 5.1.1, the MSE loss reads in our
case

$d_{\rm MSE}(T^{({\varphi})},T^{({\phi})})=\frac{1}{N_{\perp}\times
N_{z}\times
2}\sum_{i=0}^{N_{\perp}-1}\sum_{j=0}^{N_{z}-1}\sum_{i=0}^{1}\big{(}T^{({\phi})}_{ijk}-T^{({\varphi})}_{ijk}\big{)}^{2}.$
(28)

The MSE is very general and can be thought of, quite simply, as the mean
squared “distance” between the initial and reconstructed data. However, this
generality implies that it does not contain any information about the
properties of the data one tries to reconstruct.

Indeed, we can define a metric that is better suited to the physics we aim to
describe. Let us recall that our goal is to compute potential energy surfaces
that can be used, e.g., for fission simulations. These PES are nothing but
generator states for the (TD)GCM mentioned in Sec. 3.4. The GCM relies on the
norm kernel $\mathcal{N}(\bm{q},\bm{q}^{\prime})$ and the Hamiltonian kernel
$\mathcal{H}(\bm{q},\bm{q}^{\prime})$, which are defined as

$\displaystyle\mathcal{N}(\bm{q},\bm{q}^{\prime})$
$\displaystyle=\braket{\Phi({\bm{q}})}{\Phi({\bm{q}^{\prime}})},$ (29)
$\displaystyle\mathcal{H}(\bm{q},\bm{q}^{\prime})$
$\displaystyle=\braket{\Phi({\bm{q}})}{\hat{H}}{\Phi({\bm{q}^{\prime}})}.$
(30)

Since the norm kernel involves the standard inner product in the many-body
space, it represents the topology of that space. Therefore, it should be
advantageous to use for the loss a metric induced by the same inner product
that defines the norm kernel.

In our case, we want to build an AE where the encoder
$\bm{v}^{({\varphi})}=E(T^{({\varphi})})$ compresses the single-particle,
canonical orbitals $\\{\varphi_{\mu}\\}_{\mu}$ associated with $\ket{\Phi}$
into a low-dimensional vector $\bm{v}^{({\varphi})}$ and where the decoder
$T^{({\phi})}=D(\bm{v}^{({\varphi})})$ is used to compute the set of
reconstructed canonical orbitals $\\{\phi_{\mu}\\}_{\mu}$. Most importantly,
this reconstruction should be such that the reconstructed many-body state
$\ket{\Psi}$ is as close as possible to the original state $\ket{\Phi}$. In
other words, we need to use a loss that depends on the norm overlap (between
many-body states) but since we work with single-particle wavefunctions, we
must have a way to relate the norm overlap to these s.p. wavefunctions. This
can be achieved with Equations (5.4) and (5.6) of Haider and Gogny (1992),
which relate the inner product $\braket{\Phi}{\Psi}$ in the many-body space
with the inner product (overlap) $\braket{\varphi_{\mu}}{\phi_{\nu}}$ between
the related canonical orbitals $\varphi_{\mu}$ and $\phi_{\nu}$,

$\braket{\varphi_{\mu}}{\phi_{\nu}}\equiv\tau^{(\varphi\phi)}_{\mu\nu}=\\{{\hat{a}}^{(\varphi){\dagger}}_{\mu},{\hat{a}}^{(\phi)}_{\nu}\\}=\sum_{\sigma}\int\operatorname{d}^{3}\\!{\bm{r}}\,\varphi_{\mu}^{*}(\bm{r},\sigma)\phi_{\nu}(\bm{r},\sigma)$
(31)

and with the occupation amplitudes. However, it assumes that the canonical
wavefunctions of each many-body state are orthogonal. This property is not
guaranteed for our reconstructed canonical wavefunctions. In fact, because of
this lack of orthogonality, the reconstructed wavefunctions cannot be
interpreted as representing the canonical basis of the Bloch-Messiah-Zumino
decomposition of the quasiparticle vacuum and the Haider & Gogny formula
cannot be applied ‘as is’. However, we show in Appendix A that it is possible
to find a set of transformations of the reconstructed wavefunctions that
allows us to define such as genuine canonical basis.

We want the loss function to depend only on the error associated with the
reconstructed orbital $\phi_{\mu}$. Therefore, we should in principle consider
the many-body state $\ket{\tilde{\Phi}_{\mu}}$ where only the orbital
$\varphi_{\mu}$ is substituted by its reconstruction $\phi_{\mu}$. We can then
compute the inner product between $\ket{\Phi}$ and $\ket{\tilde{\Phi}_{\mu}}$
using Appendix A and deduce any induced metric $f$

$d_{\rm
exact}^{f}(T^{(\varphi)},T^{(\phi)})=f\left(\frac{\braket{\Phi}{\tilde{\Phi}_{\mu}}}{\sqrt{\braket{\tilde{\Phi}_{\mu}}{\tilde{\Phi}_{\mu}}}}\right).$
(32)

However, computing this metric is too computationally involved to be carried
out explicitly for each training data at each epoch. Instead, we keep this
metric for comparing a posteriori the performance of our model.

Instead of explicitly determining $d_{\rm
exact}^{f}(T^{(\varphi)},T^{(\phi)})$, we focus on reproducing canonical
orbitals using the metrics of the one-body Hilbert space. In practice we
considered the distance noted $d_{\circ}^{(0)}$ that is induced by the inner
product between normalized functions in the one-body Hilbert space, that is,

$d_{\circ}^{(0)}(\varphi,\phi)\equiv\left(\frac{\bra{\varphi}}{\sqrt{\braket{\varphi}{\varphi}}}-\frac{\bra{\phi}}{\sqrt{\braket{\phi}{\phi}}}\right)\left(\frac{\ket{\varphi}}{\sqrt{\braket{\varphi}{\varphi}}}-\frac{\ket{\phi}}{\sqrt{\braket{\phi}{\phi}}}\right)\,,$
(33)

which is nothing but

$d_{\circ}^{(0)}(\varphi,\phi)=\sum_{\sigma}\int\operatorname{d}^{3}\\!{\bm{r}}\big{|}\varphi(\bm{r},\sigma)-\phi(\bm{r},\sigma)\big{|}^{2}\,,$
(34)

where the $\varphi(\bm{r},\sigma)$ and $\phi(\bm{r},\sigma)$ have been
normalized. Since all wavefunctions are discretized on the Gauss quadrature
mesh, this distance reads

$d_{\circ}^{(0)}(\varphi,\phi)=\sum_{n_{\perp}n_{z}n_{\sigma}}W_{n_{\perp}n_{z}}\big{|}T^{(\varphi)}_{n_{\perp}n_{z}n_{\sigma}}-T^{(\phi)}_{n_{\perp}n_{z}n_{\sigma}}\big{|}^{2},$
(35)

where the weights $W$ are given by

$W_{n_{\perp}n_{z}}=\frac{w^{\rm GL}_{n_{\perp}}}{2b_{\perp}^{2}}\times
2\pi\times\frac{w^{\rm GH}_{n_{z}}}{b_{z}}.$ (36)

These weights, which depend on the indices $n_{\perp}$ and $n_{z}$ in the
summation, are the only difference between the squared distance loss (35) and
the MSE loss (28). Although the distance (35) is norm-invariant, it still
depends on the global phase of each orbital. We have explored other possible
options for the loss based on norm- and phase-invariant distances; see
Appendix B for a list. However, we found in our tests that the distance
$d_{\circ}^{(0)}$ systematically outperformed the other ones and, for this
reason, only show results obtained with this one.

####  5.1.4 Physics-Informed Autoencoder

From a mathematical point of view, deep neural networks can be thought of as a
series of compositions of functions. Each composition operation defines a new
layer in the network. Networks are most often built with alternating linear
and nonlinear layers. The linear part is a simple matrix multiplication.
Typical examples of nonlinear layers include sigmoid, tanh, Rectified Linear
Unit (ReLU) functions. In addition to these linear and nonlinear layers, there
could be miscellaneous manipulations of the model for more specific purposes,
such as adding batch normalization layers Ioffe and Szegedy (2015), applying
dropout Srivastava et al. (2014) to some linear layers, or skip connection He
et al. (2016) between layers.

Figure 6: Schematic example of a convolutional layer. For any $2\times 2$
chunk $C$ of the input image on the left, this convolutional layer performs
the point-wise multiplication of $C$ with the filter $F$ followed by the
addition of all elements. This compresses the initial chunk of the image into
a single integer.

Our data is a smooth function defined over a $N_{\perp}\times N_{z}=60\times
40$ grid and is analogous to a small picture. For this reason, we chose a 2D
convolutional network architecture. Convolutional layers are popular for image
analysis, because they incorporate the two-dimensional pixel arrangement in
the construction of the weights of the network. These two-dimensional weights,
or filters, capture local shapes and can model the dependent structure in
nearby pixels of image data. Given a 2D $m\times m$ input array, a 2D filter
$F$ is a $n\times n$ matrix, usually with $n\ll m$. If we note $I^{n}$ the
space of $n\times n$ integer-valued matrices, then the convolutional layer
$\mathcal{C}$ is an operation of the
$\mathcal{C}:(I^{n},I^{n})\rightarrow\mathbb{N}$ that is applied to all pairs
$(F,C)$ where $C$ is any $n\times n$ chunk of the input image; see Fig. 6 for
an example. This way, the resulting output summarizes the strength and
location of that particular filter shape within the image. As the model gets
trained, the filter parameters are fitted to a shape that is learned to be
important in the training data. Convolutional neural network are very
effective for image analysis and are currently widely used Krizhevsky et al.
(2012); Zeiler and Fergus (2014); Sermanet et al. (2013); Szegedy et al.
(2017)

In this work, we used the Resnet 18 model as our encoder and constructed the
decoder from a transposed convolution architecture of the Resnet 18. The
Resnet 18 model was first introduced by He et al. as a convolutional neural
network for image analysis He et al. (2016). It was proposed as a solution to
the degradation of performance as the network depth increases. Resnet branches
an identity-function addition layer to sub-blocks (some sequential layers of
composition) of a given network. While a typical neural network sub-block
input and output could be represented by $x$ and $f(x)$, respectively, a
Resnet sub-block would output $f(x)+x$ for the same input $x$, as in Figure 2
of He et al. (2016). This architecture is called ‘skip connection’ and was
shown to be helpful for tackling multiple challenges in training deep neural
network such as vanishing gradient problem and complex loss function Li et al.
(2018); He et al. (2016). Since then, the Resnet architecture has been widely
successful, often being used as a baseline for exploring new architectures
Zhang et al. (2022); Radosavovic et al. (2020) or as the central model for
many analyses Cubuk et al. (2020); Yun et al. (2019); Zhang et al. (2017). In
a few cases, it was also combined with autoencoders for feature learning from
high-dimensional data Wickramasinghe et al. (2021).

Figure 7: Schematic illustration of the 2D-transposed convolution. Each input
value, e.g., 55, 57, etc., is multiplied by the entire kernel resulting in a
$3\times 3$ matrix. These matrices are then added to one another in a sliding
and overlapping way. Figure 8: Representation of our modified Resnet 18
architecture for the encoder (left) and the decoder (right). Large numbers on
the left of each side label the different layers. Numbers such as 64, 128,
etc. refer to the size of the filer; see text for a discussion of some of the
main layers.

For the decoder part, we designed a near-mirror image of the encoder using
transposed convolution. Transposed convolution is essentially the opposite
operation to convolution in terms of input and output dimensions. Here the
meaning of transpose refers to the form of the filter matrix when the
convolution layer is represented by a 1D vector input obtained from
linearizing the 2D input. Note that the mirror-located filters in the decoder
are independent parameters and not the actual transposed filter matrix of the
encoder. Such a construction ensures symmetrical encoder and decoder models,
making the decoder model close to the inverse shape of the encoder model.
Figure 7 illustrates the operation: one input value is multiplied by the
entire kernel (filter) and is added to the output matrix at its corresponding
location. The corresponding output location for each colored input number are
color-coded and show how the addition is done.

The first and the last layer of the Resnet architecture are mostly for
resizing and were minimally modified from the original Resnet 18 model since
the size of our input data is significantly smaller than typical image sizes
used for Resnet image classification analyses. We also modified the number of
input channels of the first layer of the encoder to be 2 (for each of the spin
components of the nuclear wavefunction) instead of the usual number 3 (for the
RGB colors of colored images) or 1 (for black-and-white images). The spin
components are closely related to each other with covariance structure,
similar to how colors interact within an image. Therefore, we treat a pair of
spin components as a single sample and treat each component as an input
channel. The same applies to the output channel of the decoder.

The full network is represented schematically in Fig. 8. Parametric Rectified
Linear Unit, or PRELU, layers were added to impose nonlinearity in the model
He et al. (2015). PRELU layers are controlled by a single hyperparameter that
is trained with the data. Batch normalization is a standardizing layer that is
applied to each batch by computing its mean and standard deviation. It is
known to accelerate training by helping with optimization steps Ioffe and
Szegedy (2015). The average pooling layer (bottom left) averages each local
batch of the input and produces a downsized output. The upsampling layer (top
right) upsamples the input using a bilinear interpolation.

###  5.2 Training

As mentioned in Section 5.1.4, the loss is the discrepancy between the input
of the encoder and the output of the decoder. The minimization of the loss
with respect to all the model parameters $w$, such as the filter parameters,
is the training process. We used the standard back-propagation algorithm to
efficiently compute the gradient of the loss function with respect to the
model parameters. The gradient computation is done with the chain rule,
iterating from the last layer in the backward direction. We combined this with
the mini-batch gradient descent algorithm: ideally, one would need the entire
dataset to estimate the gradient at the current model parameter value.
However, with large datasets, this becomes computationally inefficient.
Instead we use a random subset of the entire data, called mini-batch, to
approximate the gradient, and expedite the convergence of the optimization.
For each mini-batch, we update each parameter $w$ by taking small steps of
gradient descent,
$w_{t+1}=w_{t}-\alpha\frac{\displaystyle\partial\mathcal{L}}{\displaystyle\partial
w_{t}}$. At step $t$, or at $t^{\mathrm{th}}$ mini-batch, the average loss
$\mathcal{L}$ and the gradient with respect to current model parameter $w_{t}$
are computed. Then $\alpha$-sized gradient descent step is taken to update the
model parameters. Instead of using the current gradient for the update, one
can use a weighted average of past gradients. We employed the well-known Adam
algorithm, which uses the exponential moving average of current and past
gradients Kingma and Ba (2014).

Iterating over the entire dataset once, using multiple mini-batches, is called
an epoch. Typically a deep neural network needs hundreds to thousands of
epochs for the algorithm to converge. Parameters such as the batch size or
learning rate, the parameters of the optimizer itself (Adam’s or other), and
the number of epochs are hyper-parameters that must be tuned for model
fitting. For our training, we used the default initialization method in
PyTorch for the model parameters. The linear layers were initialized with a
random uniform distribution over $[-1/k,1/k]$, where $k$ is the size of the
weight. For example, if there are 2 input channels and 3$\times$3 convolution
filters are used, $k=2\times 3\times 3$. PRELU layers were initialized with
their default PyTorch value of $0.25$. We proceeded with mini-batches of size
32 with the default $\beta_{1}=0.9,\beta_{2}=0.999$ and $\epsilon=10^{-8}$:
all these numbers refer to the PyTorch implementation of the Adam’s optimizer.
For $\alpha$, we used $0.001$ as starting value and used a learning rate
scheduler, which reduces the $\alpha$ value by a factor of 0.5 when there is
no improvement in the loss for 15 epochs. After careful observation of the
loss curves, we have estimated that at least 1000 epochs are needed to achieve
convergence.

To mitigate the problem of the global phase invariance of the canonical
wavefunctions discussed in Sec. 5.1.2, we doubled the size of the dataset: at
each point $\bm{q}$ of the collective space (=the sample), we added to each
canonical wavefunction $\varphi_{\mu}(\bm{r},\sigma)$ the same function with
the opposite sign $-\varphi_{\mu}(\bm{r},\sigma)$. The resulting dataset was
then first split into three components, training, validation and test
datasets, which represent 70%, 15%, and 15% of the entire data respectively.
Training data is used for minimizing the loss with respect to the model
parameters as explained above. Then we choose the model at the epoch that
performs the best with the validation dataset as our final model. Finally, the
model performance is evaluated using the test data.

###  5.3 Results

In this section, we summarize some of the preliminary results we have obtained
after training several variants of the AE. In Sec. 5.3.1, we give some details
about the training data and the quality of the reconstructed wavefunctions. We
discuss some possible tools to analyze the structure of the latent space in
Sec. 5.3.2. In these two sections, we only present results obtained for latent
spaces of dimension $D=20$. In Sec. 5.3.3, we use the reconstructed
wavefunctions to recalculate HFB observables with the code HFBTHO. We show the
results of this physics validation for both $D=20$ and $D=10$.

####  5.3.1 Performance of the Network

Figure 9 shows the initial potential energy surface in 98Zr used in this work.
Using the HFBTHO solver, we performed a total of 548 HFB calculations with
constraints on the axial quadrupole, $q_{20}=\langle\hat{Q}_{20}\rangle$, and
axial octupole moment, $q_{30}=\langle\hat{Q}_{30}\rangle$. The mesh was:
$-12.5\,\mathrm{b}\leq q_{20}\leq 25.0\,\mathrm{b}$ with steps $\delta
q_{20}=1\,\mathrm{b}$ and $0.0\,\mathrm{b}{3/2}\leq q_{30}\leq
3.0\,\mathrm{b}^{3/2}$ with step $\delta q_{30}=0.125\,\mathrm{b}^{3/2}$. The
black dots in Fig. 9 indicate the location of the converged solutions. For
each solution, the $n_{p}=60$ highest-occupation proton and $n_{n}=87$
highest-occupation neutron canonical wavefunctions were used as training data
for the network666Since time-reversal symmetry is conserved, the Fermi energy
is located around states with indices $\mu_{p}\approx 20$ and $\mu_{n}\approx
29$. Therefore, our choice implies that in our energy window, about 1/3 of all
states are below the Fermi level and about 2/3 of them are above it..

Figure 9: Potential energy surface of 98Zr in the $(q_{20},q_{30})$ plane.
Converged HFBTHO solutions are represented by black dots. Energies given by
the color bar are in MeV relatively to the ground state.

For each of the losses discussed in Appendix B, we trained the AE with the
slightly modified Resnet 18 architecture described in Sec. 5.1.4. It is
important to keep in mind that the value of these losses should not be
compared with one another. The only rigorous method to compare the performance
of both networks would be to compute the many-body norm overlap across all the
points in each case – or to perform a posteriori physics validation with the
reconstructed data, as will be shown in Sec. 5.3.3.

To give an idea of the quality of the AE, we show in Fig. 10 one example of
the original and reconstructed canonical wavefunctions. Specifically, we
consider the configuration
$(q_{20},q_{30})=(-7.0\,\mathrm{b},-0.25\,\mathrm{b}^{3/2})$ in the collective
space and look at the neutron wavefunction with occupation number
$v_{\mu}^{2}=0.945255$, which is located near the Fermi surface. This example
was obtained for an AE trained with the $d_{o}^{(0)}$ loss and compressed to
$D=20$. The figure shows, in the left panel, the logarithm of the squared norm
of the original wavefunction across the quadrature mesh,
$\ln|\varphi_{\mu}|^{2}\equiv\ln|T^{(\varphi_{\mu})}_{n_{\perp}n_{z}n_{\sigma}}|^{2}$,
in the middle panel, the same quantity for the reconstructed wavefunction, and
in the right panel the logarithm of the difference between the two. On this
example, the AE can reconstruct the wavefunction with about 3% error.

Figure 10: Left: Contour plot of the logarithm of the squared norm of the
neutron canonical wavefunction with occupation number $v_{\mu}^{2}=0.945255$
(without the exponential factor). Middle: Same for the reconstructed
wavefunction. Right: Logarithm of the difference between the squared norm of
the original and reconstructed wavefunctions.

####  5.3.2 Structure of the Latent Space

One of the advantages of AEs is the existence of a low-dimensional
representation of the data. In principle, any visible structure in this latent
space would be the signal that the network has properly learned, or encoded,
some dominant features of the dataset. Here, our latent space has dimension
$D=20$. This means that every canonical wavefunction, which is originally a
matrix of size $n=N_{\perp}\times N_{z}$, is encoded into a single vector of
size $D$. From a mathematical point of view, the encoder is thus a function

$\begin{array}[]{rl}\hat{E}:\mathbb{R}^{n}&\longrightarrow\mathbb{R}^{D}\\\\[-2.0pt]
\varphi&\longmapsto\bm{v}=\hat{E}(\varphi)\end{array}$ (37)

Let us consider some (scalar) quantity $P$ associated with the many-body state
$\ket{\Phi(\bm{q})}$ at point $\bm{q}$. Such a quantity could be an actual
observable such as the total energy but it could also be an auxiliary object
such as the expectation value of the multipole moment operators. In fact, $P$
could also be a quantity associated with the individual degrees of freedom at
point $\bm{q}$, for example the q.p. energies. In general terms, we can think
of $P$ as the output value of the function

$\begin{array}[]{rl}\hat{P}:\mathbb{R}^{n}&\longrightarrow\mathbb{R}\\\\[-2.0pt]
\varphi&\longmapsto P=\hat{P}(\varphi)\end{array}$ (38)

For example, if $P$ represents the s.p. canonical energies, then the function
$\hat{P}$ is the one that associates with each canonical wavefunction its s.p.
energy. Therefore, for every canonical wavefunction, there is a different
value of $P$. Conversely, if $P=\braket{\hat{Q}_{20}}$, there is a single
value for all the canonical wavefunctions at point $\bm{q}$. Since there is a
vector in the latent space for each canonical wavefunction, and there is also
a value for the quantity $P$ for each such function, we can then define the
new function $\hat{\mathcal{P}}$ acting on vectors of the latent space and
defined as

$\begin{array}[]{rl}\hat{\mathcal{P}}:\mathbb{R}^{D}&\longrightarrow\mathbb{R}\\\\[-2.0pt]
\bm{v}&\longmapsto P=\hat{\mathcal{P}}(\bm{v})\end{array}$ (39)

and it is straightforward to see that:
$\hat{P}=\hat{\mathcal{P}}\circ\hat{E}$. Our goal is now to try to analyze
where various quantities $P$ are located in the latent space and whether one
can identify some specific features of these locations.

Since we have a total of $n_{t}=147$ wavefunctions for each of the $N_{p}=552$
points in the collective space, the encoder yields a set of $n_{t}\times
N_{p}$ vectors of dimension $D$. This means that, in the latent space, every
quantity $P$ above is also represented by a cloud of $n_{t}\times N_{p}$ such
vectors. This is obviously impossible to visualize. For this reason, we
introduce the following analysis. First, we perform a linear regression in the
$D$-dimensional latent space of a few select quantities of interest $P$, that
is, we write

$P=\bm{\alpha}\cdot\bm{v}+b\,,$ (40)

where $\bm{\alpha}$ is a $D$-dimensional vector, $\bm{v}$ is the vector
associated with the quantity $P$ in the latent space and $b\in\mathbb{R}$. The
unit vector $\bm{u}=\bm{\alpha}/||\bm{\alpha}||$ can be interpreted as
representing the leading direction in the latent space. The quantity
$u=\bm{u}\cdot\bm{v}$ is a scalar which we obtain easily from the result of
the linear regression. We can thus plot the function $P:u\mapsto P(u)$.
Examples of such functions are shown in Fig. 11. Each point in the figures
represents the value $P=\hat{\mathcal{P}}(\bm{v})$ of some characteristic
quantity at point $u=\bm{u}\cdot\bm{v}$.

Figure 11: One-dimensional projections of the $D$-dimensional linear fit for
the total energy $E_{\rm HFB}$ (left panel), the projection $\Omega$ of the
canonical state (middle panel) and the neutron Fermi energy $\lambda_{n}$
(right panel). Each point represents one of these quantities for a canonical
wavefunction $\mu$ and a point $\bm{q}$ in the collective space.

The three cases shown in Fig. 11 illustrate that the network has not always
identified relevant features. The case of $\Omega$, middle panel, is the
cleanest: there is a clear slope as a function of $u$: if one sets $u=1$, for
example, then only values of $7/2\leq\Omega\leq 15/2$ are possible.
Conversely, the AE has not really discovered any feature in the neutron Fermi
energy (right panel): for any given value of $u$, there is a large range of
possible values of Fermi energies. In the case of the total energy (left
panel), the situation is somewhat intermediate: there is a faint slope
suggesting a linear dependency of the energy as a function of $u$.

####  5.3.3 Physics Validation

The results presented in the Section 5.3.1 suggest the AE has the ability to
reproduce the canonical wavefunctions with good precision. To test this
hypothesis, we recalculated the HFB solution at all the training, validation
and testing points by substituting in the HFBTHO binary files the original
canonical wavefunctions by the ones reconstructed by the AE. Recall that only
the lowest $n_{t}$ wavefunctions with the largest occupation were encoded in
the AE ($n_{n}=87$ for neutrons and $n_{p}=60$ for protons); the remaining
ones were unchanged. In practice, their occupation is so small that their
contribution to nuclear observables is very small ($<$ 10 keV for the total
energy, for example).

Figure 12: Left: Potential energy surface in the $(q_{20},q_{30})$ plane for
98Zr obtained after replacing the first $n_{n}=87$ and $n_{p}=60$ highest-
occupation canonical wavefunctions by their values reconstructed by the AE for
a latent space of dimension $D=20$. The black dots show the location of the
training points only, the white circles the location of the validation points.
Right: same figure for a latent space of dimension $D=10$. For both figures,
energies are given in MeV.

Figure 12 shows the error on the potential energy across the $(q_{20},q_{30})$
collective space obtained with the reconstructed canonical wavefunctions for
latent spaces of dimension $D=20$ (left) and $D=10$ (right). In each case, we
only show results obtained when using the $d_{o}^{(0)}$ loss, which gives the
best results. The black crosses denote the location of all the original
points; the white circles show the location of the validation points. Overall,
the results are very encouraging. In both cases, most of the error is
concentrated near regions of the PES where there are discontinuities (hence,
the lack of converged solutions). Everywhere else, the error is small and
mostly randomly distributed across the PES, that is, it is not systematically
larger at the validation points. As expected, the quality of the
reconstruction is a little worse when $D=10$: one can notice about a dozen of
points for which the error is significantly larger, in absolute value.
Examples include $(q_{20},q_{30})=(-5.0\,\rm{b},1.125\,\rm{b}^{3/2})$ or
$(q_{20},q_{30})=(+8.0\,\rm{b},2.5\,\rm{b}^{3/2})$ in the validation set, and
$(q_{20},q_{30})=(0.0\,\rm{b},1.75\,\rm{b}^{3/2})$ or the region around
$1\,\rm{b}\leq q_{20}\leq 4\,\rm{b}$ and $1.5\,\rm{b}^{3/2}\leq q_{30}\leq
2.25\,\rm{b}^{3/2}$ in the training set. These may suggest that for $D=10$,
the loss may not have fully converged yet. Because of the existence of
discontinuities near these points, this could also be the manifestation that
our continuous AE cannot build a continuous representation of the data
everywhere. However, the fact that an increase of the compression by a factor
2, from $D=20$ to $D=10$, does not substantially degrade the performance of
the AE is very promising.

Figure 13: Left: Histogram of the difference in total HFB energy between the
original HFBTHO calculation and the result obtained by computing the energy in
the canonical basis with the reconstructed wavefunctions (see text for
details). Calculations were performed both for a $D=20$ and $D=10$ latent
space. Right: Similar histogram for the expectation value of the axial
quadrupole moment.

The two histograms in Fig. 13 give another measure of the quality of the AE.
The histogram in the left shows the distribution of the error on the HFB
energy for two sizes of the latent space, $D=20$ and $D=10$. In both cases,
most of the error is less than $\pm 100$ keV, which is rather good. As
mentioned before, the points with the higher error $\epsilon_{E}<-150$ keV or
$\epsilon_{E}>250$ keV do not correspond only to testing points. The histogram
on the right shows the distribution of the error for $q_{20}$ and $q_{30}$ (in
units of b and $b^{3/2}$, respectively).

##  6 Conclusion

Extreme-scale calculations based on nuclear density functional methods
relevant for, e.g., nuclear astrophysics simulations or uncertainty
quantification remain computationally expensive and fraught with formal and
practical issues associated with self-consistency or reduced collective
spaces. In this article, we have analyzed two different techniques to build
fast, efficient and accurate surrogate models, or emulators, or DFT objects.

We first showed that Gaussian processes could reproduce reasonably well the
values of the mean-field and pairing-field potentials of the HFB theory across
a large two-dimensional potential energy surface. The absolute error on the
total energy was within $\pm 100$ keV and the relative errors on the
collective inertia tensor smaller than 5%. However, GPs require the training
data to be “smoothly-varying”, i.e., they should not include phenomena such as
nuclear scission or, more generally, discontinuities in the PES. It is well
known that GPs are not reliable for extrapolation: such a technique can thus
be very practical to densify (=interpolate) an existing potential energy
surface but must not be applied outside its training range.

Although standard versions of GPs are fast and simple to use, incorporating
more of the many of the existing correlations between the values of the HFB
potentials may considerably increase the complexity of the emulator. In such a
case, it is more natural to use directly deep-learning techniques. In this
work, we reported the first application of autoencoders to emulate the
canonical wavefunctions of the HFB theory. Autoencoders are a form of deep
neural network that compresses the input data, here the canonical
wavefunctions, into a small-dimensional space called the latent
representation. The encoder is trained simultaneously with a decoder by
enforcing that the training data is left invariant after compression followed
by decompression. In practice, the measure of such “invariance” is set by what
is called the loss of the network. We discussed possible forms of the loss
that are best adapted to learning quantum-mechanical wavefunctions of many-
body systems such as nuclei. We showed that such an AE could successfully
reduce the data into a space of dimension $D=10$ while keeping the total error
on the energy lower than $\Delta E=150$ keV (on average). The analysis of the
latent space revealed well-identified structures in a few cases, which
suggests the network can learn some of the physics underlying the data. This
exploratory study suggests that AE could serve as reliable canonical
wavefunctions generators. The next step will involve learning a full sequence
of such wavefunctions, i.e., an ordered list, in order to emulate the full HFB
many-body state.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

M. Verriere led the study of autoencoders for canonical wavefunctions: design
of the overall architecture of the autoencoders, development of the software
stack and analysis of the results. N. Schunck led the study of Gaussian
process for mean-field potentials – training and validation runs, code
development and HFBTHO calculations – and supervised the whole project. I. Kim
implemented, tested and trained different architectures of autoencoders and
helped with the analysis of the results. P. Marević developed an HFBTHO module
to use canonical wavefunctions in HFB calculations. K. Quinlan provided
technical expertise about Gaussian processes and supervised M. N’Go, who
implemented and fitted Gaussian processes to mean-field potentials. D. Regnier
and R.D. Lasseri provided technical consulting on deep neural network and the
architecture of autoencoders. All authors contributed to the writing of the
manuscript, and read and approved the submitted version.

## Funding

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.
Computing support came from the Lawrence Livermore National Laboratory (LLNL)
Institutional Computing Grand Challenge program.

##  Appendix A Overlap Between Non-Orthogonal Canonical Bases

In Haider and Gogny (1992) the norm overlap between two different HFB vacuua
is expressed as a function of the single-particle overlap between the two
respective sets of canonical wavefunctions and occupation numbers. This
formula implicitly assumes that the canonical wavefunctions are orthornormal.
When using canonical wavefunctions reconstructed by the AE, this property may
not hold anymore and the Haider & Gogny formula should not be used ‘as is’
when evaluating the norm overlap. In this appendix, we show how to generalize
it.

We recall that for a system with conserved time-reversal symmetry, the
quasiparticle vacuum can be written Ring and Schuck (2004); Schunck (2019)

$\ket{\Phi}=\prod_{\mu_{>}}\beta^{\dagger}_{\mu}\beta^{\dagger}_{\bar{\mu}}\ket{0},$
(41)

where the quasiparticle operators $(\hat{\beta},\hat{\beta}^{\dagger})$ are
obtained from a single-particle basis $({\hat{c}},{\hat{c}}^{\dagger})$ by the
Bogoliubov transformation $W$ of (14) in Sec. 3.5. The Bloch-Messiah-Zumino
decomposition of the matrix $W$ is a consequence of the fact that the
quasiparticle operators should obey the same anticommutation relations as the
particle operators. In the resulting canonical basis, the quasiparticle vacuum
can be written in the BCS form,

$\ket{\Phi}=\prod_{\mu_{>}}\left(u_{\mu}+v_{\mu}{\hat{a}}_{\mu}^{\dagger}{\hat{a}}_{\bar{\mu}}^{\dagger}\right)\ket{0},$
(42)

where $\ket{0}$ is the particle vacuum, $\bar{\mu}$ refers to the time-
reversed partner of the state with index $\mu$, and the summation runs only
over “positive” indices $\mu$. Implicit in this expression is the fact that
the canonical wavefunctions associated with the operators
$({\hat{a}},{\hat{a}}^{\dagger})$ are orthonormal or, equivalently, that the
operators $({\hat{a}},{\hat{a}}^{\dagger})$ anticommute. If these conditions
are not verified, the form (42) is not valid and the formulas for the norm
overlap given in Haider and Gogny (1992) cannot apply. In our case, the fit of
the AE gives a set of reconstructed canonical wavefunctions which we noted
$\\{\phi_{\mu}(\bm{r},\sigma)\\}_{\mu}$ and are associated with a set of
single-particle creation and annihilation operators
$({\hat{f}},{\hat{f}}^{\dagger})$. Although we call these objects canonical
orbitals, this is somewhat a misnomer since the wavefunctions are not
necessarily orthonormal. As a consequence, one cannot define a BCS state (42)
with these operators. Our goal is to find a transformation of these single-
particle operators that allows us to define a BCS state.

Following the notations of Haider and Gogny (1992), we thus define the single-
particle overlaps

$\tilde{\tau}{}_{\mu\nu}^{(ff)}=\\{{\hat{f}}_{\mu}^{\dagger},{\hat{f}}_{\nu}\\}=\sum_{\sigma}\int\operatorname{d}^{3}\\!\bm{r}\,\phi_{\mu}^{*}(\bm{r},\sigma)\phi_{\nu}(\bm{r},\sigma).$
(43)

The set of all such overlaps define the overlap matrix
$\tilde{\tau}{}^{(ff)}$. This overlap matrix is block diagonal as it satisfies
the relations

$\tilde{\tau}{}^{(ff)}_{\mu\nu}={\tilde{\tau}{}^{(ff)}_{\bar{\mu}\bar{\nu}}}^{*},\qquad\tilde{\tau}{}^{(ff)}_{\mu\bar{\nu}}=\tilde{\tau}{}^{(ff)}_{\bar{\mu}\nu}=0.$
(44)

From the s.p. operators $({\hat{f}},{\hat{f}}^{\dagger})$, we can introduce a
new set of q.p. operators $({\hat{\chi}},{\hat{\chi}}^{\dagger})$ through

$\displaystyle{\hat{\chi}}_{\mu}$
$\displaystyle=\tilde{u}_{\mu}{\hat{f}}_{\mu}-\tilde{v}_{\mu}{\hat{f}}_{\bar{\mu}}^{\dagger},$
(45a) $\displaystyle{\hat{\chi}}_{\bar{\mu}}$
$\displaystyle=\tilde{u}_{\mu}{\hat{f}}_{\bar{\mu}}+\tilde{v}_{\mu}{\hat{f}}_{\mu}^{\dagger}.$
(45b)

It is easy to see that these q.p. operators do not obey the Fermion
anticommutation relation. In other words, the q.p. spinors associated with
these operators are not orthogonal. We thus introduce the overlap matrix
$\tilde{\tau}{}^{(\chi\chi)}_{\mu\nu}$ between any two such spinors
$\mu,\nu>0$. Owing to (45a)-(45b) and (43) it is straightforward to show that
it is given by

$\tilde{\tau}{}^{(\chi\chi)}_{\mu\nu}=\tilde{u}_{\mu}\tilde{u}_{\nu}\tilde{\tau}{}^{(ff)}_{\mu\nu}+\tilde{v}_{\mu}\tilde{v}_{\nu}{\tilde{\tau}{}^{(ff)}_{\mu\nu}}^{*},$
(46)

and verify the same properties (44) as the single-particle overlap. We
symmetrically orthogonalize the q.p. basis by eigendecomposing
$\tilde{\tau}{}^{(\chi\chi)}_{\mu\nu}$ 777In fact, we can limit ourself to
compute the Cholesky decomposition of $\tilde{\tau}{}^{(\chi\chi)}_{\mu\nu}$,
but we use the eigenvalues to check the rank and invert the matrix.

$\tilde{\tau}{}^{(\chi\chi)}_{\mu\nu}=Q\Sigma^{2}Q^{\dagger},$ (47)

where $Q$ is unitary and
$\Sigma=\operatorname{diag}(\sigma_{0},\sigma_{1},\dots)$ with
$\sigma_{\mu}>0$. We then construct a new orthogonal q.p. basis

${\hat{\gamma}}_{\mu}^{\dagger}=\sum_{k}{\hat{\chi}}_{k}^{\dagger}\big{(}Q\Sigma^{-1}\big{)}_{k\mu}$
(48)

that satisfies the fermion commutation relations,
$\\{{\hat{\gamma}}_{\mu}^{\dagger},{\hat{\gamma}}_{\nu}\\}=\delta_{\mu\nu}$.
We can associate with these new q.p. operators
$({\hat{\gamma}}^{\dagger},{\hat{\gamma}})$ a quasiparticle vacuum of the type
(41). We now need to find the Bogoliubov transformation $\tilde{W}$ (and its
Bloch-Messiah decomposition) that relate the
$({\hat{\gamma}}^{\dagger},{\hat{\gamma}})$ to a properly orthonormal s.p.
basis. To this end, we first diagonalize the single-particle overlap matrix

$\tilde{\tau}{}^{(ff)}_{\mu\nu}=R\tilde{\Sigma}^{2}R^{\dagger},$ (49)

which defines a new set of particle operators
$({\hat{b}},{\hat{b}}^{\dagger})$ through the relations

$\displaystyle{\hat{b}}_{i}^{\dagger}$
$\displaystyle=\sum_{k}{\hat{f}}_{k}^{\dagger}\big{(}R\tilde{\Sigma}^{-1}\big{)}_{ki}$
(50a) $\displaystyle{\hat{b}}_{\bar{i}}$
$\displaystyle=\sum_{k}{\hat{f}}_{\bar{k}}\big{(}R\tilde{\Sigma}^{-1}\big{)}_{ki},$
(50b)

By construction these new particle operators also satisfy the Fermion anti-
commutation relations,
$\\{{\hat{b}}_{i}^{\dagger},{\hat{b}}_{j}\\}=\delta_{ij}$. By inverting
relations (50a)-(50b), using the expression (45a)-(45b) relating the
$({\hat{\chi}}^{\dagger},{\hat{\chi}})$ to the
$({\hat{f}}^{\dagger},{\hat{f}})$ and using (48), these new particle operators
can be related to the q.p. operators
$({\hat{\gamma}}^{\dagger},{\hat{\gamma}})$ through

$\displaystyle{\hat{\gamma}}_{\mu}^{\dagger}$
$\displaystyle=\sum_{l}{\hat{b}}_{l}^{\dagger}\left[\tilde{\Sigma}R^{\dagger}\tilde{u}Q\Sigma^{-1}\right]_{l\mu}-{\hat{b}}_{\bar{l}}\left[\tilde{\Sigma}R^{\dagger}\tilde{v}Q\Sigma^{-1}\right]_{l\mu}$
(51a) $\displaystyle{\hat{\gamma}}_{\bar{\mu}}^{\dagger}$
$\displaystyle=\sum_{l}{\hat{b}}_{\bar{l}}^{\dagger}\left[\tilde{\Sigma}R^{\top}\tilde{u}Q^{*}\Sigma^{-1}\right]_{l\mu}+{\hat{b}}_{l}\left[\tilde{\Sigma}R^{\top}\tilde{v}Q^{*}\Sigma^{-1}\right]_{l\mu}.$
(51b)

These two equations are the main result of this appendix. They show that we
can extract from the non-orthogonal reconstructed, canonical wavefunctions a
set of quasiparticle operators that obey the Fermion anticommutation relation,
define a quasiparticle vacuum and are related to an orthonormal basis of the
single-particle Hilbert space through the following Bogoliubov transformation

$\tilde{W}=\left(\begin{array}[]{cc}\tilde{\Sigma}R^{\top}\tilde{u}{}Q^{*}\Sigma^{-1}&\tilde{\Sigma}R^{\dagger}\tilde{v}{}Q\Sigma^{-1}\\\
\tilde{\Sigma}R^{\top}\tilde{v}{}Q^{*}\Sigma^{-1}&\tilde{\Sigma}R^{\dagger}\tilde{u}{}Q\Sigma^{-1}\end{array}\right)=\left(\begin{array}[]{cc}\tilde{U}&\tilde{V}^{*}\\\
\tilde{V}&\tilde{U}^{*}\end{array}\right)\,.$ (52)

This matrix only depends on the initial canonical occupations $\tilde{u}{}$
and $\tilde{v}{}$, as well as on the eigenvalues and eigenvectors of both the
s.p. overlap matrix (43) and the q.p. overlap matrix (46).

From the new transformation $\tilde{W}$, we can define the one-body density
matrix

$\rho=\tilde{V}^{*}\tilde{V}^{\top}=\tilde{\Sigma}R^{\dagger}\tilde{v}{\tilde{\tau}{}^{(\chi\chi)}}^{-1}\tilde{v}R\tilde{\Sigma}.$
(53)

and put it into canonical form by diagonalizing it

$\rho=\tilde{D}v^{2}\tilde{D}^{\dagger}.$ (54)

The transformation $\tilde{D}$ defines the new canonical basis. By
construction, these new canonical wavefunctions are expressed in the
$({\hat{b}},{\hat{b}}^{\dagger})$ basis, which is itself related to the
original, non-orthogonal basis of the reconstructed wavefunctions
$({\hat{f}},{\hat{f}}^{\dagger})$ through (50a)-(50b). One can easily show
that we have

$\tilde{D}_{n}(\bm{r},\sigma)=\sum_{k}\phi_{k}(\bm{r},\sigma)\big{(}R\tilde{\Sigma}^{-1}\tilde{D}\big{)}_{kn}\,.$
(55)

At this point, we have obtained a set of genuine canonical wavefunctions
$\tilde{D}_{n}(\bm{r},\sigma)$ that are orthonormal and are associated with
the new occupations $v_{n}$ defined by (54). The relation between these
canonical wavefunctions and the wavefunctions reconstructed by the AE is given
by (55). Thanks to this expression, we can now apply the Haider & Gogny
formulas for the norm overlap between two many-body states $\ket{\Phi}$ and
$\ket{\Psi}$. We find

$\braket{\Phi}{\Psi}=\det\big{(}{\Sigma^{(\Phi)}}^{-1}\big{)}\det\big{(}{\Sigma^{(\Psi)}}^{-1}\big{)}\det\big{(}\tau^{\Phi\Psi}\big{)}\det\big{(}\tilde{Z}^{(\Phi\Psi)}\big{)},$
(56)

where

$\displaystyle\tau_{mn}^{(\Phi\Psi)}$
$\displaystyle=\sum_{\sigma}\int\operatorname{d}^{3}\\!{\bm{r}}\,\tilde{D}_{m}^{(\Phi)*}(\bm{r},\sigma)\tilde{D}_{n}^{(\Psi)}(\bm{r},\sigma),$
(57a) $\displaystyle Z^{(\Phi\Psi)}$
$\displaystyle=u^{(\Phi)}\big{(}\tau^{(\Phi\Psi)\dagger}\big{)}^{-1}u^{(\Psi)}+v^{(\Phi)}\tau^{(\Phi\Psi)}v^{(\Psi)},$
(57b) $\displaystyle u^{(\Phi/\Psi)}$
$\displaystyle=\sqrt{1-v^{(\Phi/\Psi)2}}.$ (57c)

##  Appendix B Metric Induced by an Inner Product

We present in this section the different notions of distance associated with
an inner product. We note $\braket{\bm{a}}{\bm{b}}$ the inner product between
two vectors $\bm{a}$ and $\bm{b}$. There are many examples of inner product in
nuclear physics, such as the overlap $\braket{\Phi}{\Psi}$ between two many-
body states $\ket{\Phi}$ and $\ket{\Psi}$ or the overlap between single-
particle orbitals $\varphi(\bm{r},\sigma)$ and $\phi(\bm{r},\sigma)$ defined
as

$\braket{\varphi}{\phi}=\sum_{\sigma}\int\operatorname{d}^{3}\\!{\bm{r}}\,\varphi^{*}(\bm{r},\sigma)\phi(\bm{r},\sigma).$
(58)

Let us first recall some standard mathematics notations. The norm induced by
the inner product $\braket{\bm{a}}{\bm{b}}$ is defined in the usual way as

$\|\bm{a}\|=\sqrt{\braket{\bm{a}}{\bm{a}}}.$ (59)

We can then introduce the distance induced by the inner product as

$d_{\rm I}(\bm{a},\bm{b})=\|\bm{a}-\bm{b}\|.$ (60)

Note that this distance depends on the possible phase and norm of $\bm{a}$ and
$\bm{b}$. However, quantum-mechanical observables do not depend on either of
them. The norm-independent distance thus reads

$d_{\circ}(\bm{a},\bm{b})=\left\|\frac{\bm{a}}{\|\bm{a}\|}-\frac{\bm{b}}{\|\bm{b}\|}\right\|,$
(61)

which can be rewritten as a function of the inner product between the two
normed vectors

$d_{\circ}(\bm{a},\bm{b})=\sqrt{2}\sqrt{1-\Re\left[\Braket{\frac{\bm{a}}{\lVert\bm{a}\rVert}}{\frac{\bm{b}}{\lVert\bm{b}\rVert}}\right]}=\sqrt{2}\sqrt{1-\left|\Braket{\frac{\bm{a}}{\lVert\bm{a}\rVert}}{\frac{\bm{b}}{\lVert\bm{b}\rVert}}\right|\cos\Theta},$
(62)

where $\Theta=\arg{\braket{\bm{a}}{\bm{b}}}$. Another choice for a norm-
invariant distance is the Great-Circle distance, also known as orthodromic or
spherical distance, that is defined as Deza and Deza (2009)

$d_{\perp}(\bm{a},\bm{b})=\arccos\left(\Re\left[\Braket{\frac{\bm{a}}{\lVert\bm{a}\rVert}}{\frac{\bm{b}}{\lVert\bm{b}\rVert}}\right]\right)=\arccos\left(\left|\Braket{\frac{\bm{a}}{\lVert\bm{a}\rVert}}{\frac{\bm{b}}{\lVert\bm{b}\rVert}}\right|\cos\Theta\right)$
(63)

The orthodromic distance is defined on the manifold of unit vectors. In the
case of real vector spaces, it can be interpreted as the angle between
$\bm{a}$ and $\bm{b}$. Equations (62) and (63) clearly show that both
distances $d_{\circ}$ and $d_{\perp}$ still depend on the phase $\Theta$
between $\bm{a}$ and $\bm{b}$. To remove this dependency, we minimize each
distance $d_{\circ}$ and $d_{\perp}$ over $\Theta$. This gives the following
two norm- and phase-independent distances 888Note that they are distances over
the projective space $P(\mathcal{H})$, not over the 1-body Hilbert space
$\mathcal{H}$,

$\displaystyle D(\bm{a},\bm{b})$
$\displaystyle=\sqrt{2}\sqrt{1-\left|\Braket{\frac{\bm{a}}{\lVert\bm{a}\rVert}}{\frac{\bm{b}}{\lVert\bm{b}\rVert}}\right|},$
(64a)

and

$\displaystyle D_{\perp}(\bm{a},\bm{b})$
$\displaystyle=\arccos\left|\Braket{\frac{\bm{a}}{\lVert\bm{a}\rVert}}{\frac{\bm{b}}{\lVert\bm{b}\rVert}}\right|.$
(65a)

The distance $D_{\perp}$ is an intrinsic metric and is named the Fubini–Study
metric. It is a generalization of the Bloch sphere. Table 1 presents all the
distances on the 1-body Hilbert space between orbitals that we have considered
in this work.

|  | Invariance |
---|---|---|---
Notation | Space | Norm | Phase | Definition
$d_{\rm I}$ | $\mathcal{H}$ | No | No | (60)
$d_{\circ}$ | Unit vectors of $\mathcal{H}$ | Yes | No | (61)
$d_{\perp}$ | Unit vectors of $\mathcal{H}$ | Yes | No | (63)
$D$ | Riemann sphere of $\mathcal{H}$ | Yes | Yes | (64)
$D_{\perp}$ | Riemann sphere of $\mathcal{H}$ | Yes | Yes | (65)
Table 1: Different metrics can be defined on the set of orbitals.

All these distances are defined on the one-body Hilbert space of s.p.
wavefunctions. As a result, they do not depend on the occupation probability
of canonical orbitals, in contrast to the many-body state which takes the BCS
form. As already mentioned in the main text, determining such a dependency
exactly from (32) is not trivial and computationally demanding. Instead, we
can adopt the approximation that the dependency should be proportional to some
power $p$ of the occupation number $v_{\mu}^{2}$ associated with the current
orbital,

$d^{(p)}(\varphi,\phi)=\big{(}v_{\mu}^{2}\big{)}^{p}\times d(\varphi,\phi).$
(66)

## References

  * Eschrig (1996) Eschrig R.  Fundamentals of Density Functional Theory (Leipzig: Teubner) (1996).
  * Schunck (2019) Schunck N.  Energy Density Functional Methods for Atomic Nuclei..  IOP Expanding Physics (Bristol, UK: IOP Publishing) (2019).
  * Schunck and Regnier (2022) Schunck N, Regnier D.  Theory of nuclear fission.  Prog. Part. Nucl. Phys. 125 (2022) 103963.  10.1016/j.ppnp.2022.103963.
  * Schunck and Robledo (2016) Schunck N, Robledo LM.  Microscopic theory of nuclear fission: A review.  Rep. Prog. Phys. 79 (2016) 116301.  10.1088/0034-4885/79/11/116301.
  * Kejzlar et al. (2020) Kejzlar V, Neufcourt L, Nazarewicz W, Reinhard PG.  Statistical aspects of nuclear mass models.  J. Phys. G: Nucl. Part. Phys. 47 (2020) 094001.  10.1088/1361-6471/ab907c.
  * Schunck et al. (2020a) Schunck N, O’Neal J, Grosskopf M, Lawrence E, Wild SM.  Calibration of energy density functionals with deformed nuclei.  J. Phys. G: Nucl. Part. Phys. 47 (2020a) 074001\.  10.1088/1361-6471/ab8745.
  * Erler et al. (2012) Erler J, Birge N, Kortelainen M, Nazarewicz W, Olsen E, Perhac AM, et al.  The limits of the nuclear landscape.  Nature 486 (2012) 509.  10.1038/nature11188.
  * Ney et al. (2020) Ney EM, Engel J, Li T, Schunck N.  Global description of $\beta^{-}$ decay with the axially deformed Skyrme finite-amplitude method: Extension to odd-mass and odd-odd nuclei.  Phys. Rev. C 102 (2020) 034326.  10.1103/PhysRevC.102.034326.
  * Mumpower et al. (2016) Mumpower MR, Surman R, McLaughlin GC, Aprahamian A.  The impact of individual nuclear properties on r-process nucleosynthesis.  Prog. Part. Nucl. Phys. 86 (2016) 86.  10.1016/j.ppnp.2015.09.001.
  * Perlińska et al. (2004) Perlińska E, Rohoziński SG, Dobaczewski J, Nazarewicz W.  Local density approximation for proton-neutron pairing correlations: Formalism.  Phys. Rev. C 69 (2004) 014316.  10.1103/PhysRevC.69.014316.
  * Utama et al. (2016) Utama R, Piekarewicz J, Prosper HB.  Nuclear mass predictions for the crustal composition of neutron stars: A Bayesian neural network approach.  Phys. Rev. C 93 (2016) 014311.  10.1103/PhysRevC.93.014311.
  * Utama and Piekarewicz (2017) Utama R, Piekarewicz J.  Refining mass formulas for astrophysical applications: A Bayesian neural network approach.  Phys. Rev. C 96 (2017) 044308.  10.1103/PhysRevC.96.044308.
  * Utama and Piekarewicz (2018) Utama R, Piekarewicz J.  Validating neural-network refinements of nuclear mass models.  Phys. Rev. C 97 (2018) 014306.  10.1103/PhysRevC.97.014306.
  * Niu and Liang (2018) Niu ZM, Liang HZ.  Nuclear mass predictions based on Bayesian neural network approach with pairing and shell effects.  Phys. Lett. B 778 (2018) 48.  10.1016/j.physletb.2018.01.002.
  * Neufcourt et al. (2019) Neufcourt L, Cao Y, Nazarewicz W, Olsen E, Viens F.  Neutron Drip Line in the Ca Region from Bayesian Model Averaging.  Phys. Rev. Lett. 122 (2019) 062502.  10.1103/PhysRevLett.122.062502.
  * Lovell et al. (2022) Lovell AE, Mohan AT, Sprouse TM, Mumpower MR.  Nuclear masses learned from a probabilistic neural network.  Phys. Rev. C 106 (2022) 014305.  10.1103/PhysRevC.106.014305.
  * Mumpower et al. (2022) Mumpower MR, Sprouse TM, Lovell AE, Mohan AT.  Physically interpretable machine learning for nuclear masses.  Phys. Rev. C 106 (2022) L021301.  10.1103/PhysRevC.106.L021301.
  * Niu et al. (2019) Niu ZM, Liang HZ, Sun BH, Long WH, Niu YF.  Predictions of nuclear $\beta$-decay half-lives with machine learning and their impact on r-process nucleosynthesis.  Phys. Rev. C 99 (2019) 064307.  10.1103/PhysRevC.99.064307.
  * Wang et al. (2019) Wang ZA, Pei J, Liu Y, Qiang Y.  Bayesian Evaluation of Incomplete Fission Yields.  Phys. Rev. Lett. 123 (2019) 122501.  10.1103/PhysRevLett.123.122501.
  * Lovell et al. (2020) Lovell AE, Mohan AT, Talou P.  Quantifying uncertainties on fission fragment mass yields with mixture density networks.  J. Phys. G: Nucl. Part. Phys. 47 (2020) 114001.  10.1088/1361-6471/ab9f58.
  * Lasseri et al. (2020) Lasseri RD, Regnier D, Ebran JP, Penon A.  Taming Nuclear Complexity with a Committee of Multilayer Neural Networks.  Phys. Rev. Lett. 124 (2020) 162502.  10.1103/PhysRevLett.124.162502.
  * Parr and Yang (1989) Parr R, Yang W.  Density Functional Theory of Atoms and Molecules.  International Series of Monographs on Chemistry (New York: Oxford University Press) (1989).
  * Dreizler and Gross (1990) Dreizler R, Gross E.  Density Functional Theory: An Approach to the Quantum Many-Body Problem (Springer-Verlag) (1990).  10.1007/978-3-642-86105-5.
  * Engel (2007) Engel J.  Intrinsic-density functionals.  Phys. Rev. C 75 (2007) 014306.  10.1103/PhysRevC.75.014306.
  * Barnea (2007) Barnea N.  Density functional theory for self-bound systems.  Phys. Rev. C 76 (2007) 067302.  10.1103/PhysRevC.76.067302.
  * Engel et al. (1975) Engel YM, Brink DM, Goeke K, Krieger SJ, Vautherin D.  Time-dependent Hartree-Fock theory with Skyrme’s interaction.  Nucl. Phys. A 249 (1975) 215.  10.1016/0375-9474(75)90184-0.
  * Dobaczewski and Dudek (1996) Dobaczewski J, Dudek J.  Time-Odd Components in the Rotating Mean Field and Identical Bands.  Acta Phys. Pol. B 27 (1996) 45.
  * Bender et al. (2003) Bender M, Heenen PH, Reinhard PG.  Self-consistent mean-field models for nuclear structure.  Rev. Mod. Phys. 75 (2003) 121.  10.1103/RevModPhys.75.121.
  * Lesinski et al. (2007) Lesinski T, Bender M, Bennaceur K, Duguet T, Meyer J.  Tensor part of the Skyrme energy density functional: Spherical nuclei.  Phys. Rev. C 76 (2007) 014312.  10.1103/PhysRevC.76.014312.
  * Schunck et al. (2020b) Schunck N, Quinlan KR, Bernstein J.  A Bayesian analysis of nuclear deformation properties with Skyrme energy functionals.  J. Phys. G: Nucl. Part. Phys. 47 (2020b) 104002\.  10.1088/1361-6471/aba4fa.
  * Dobaczewski et al. (1984) Dobaczewski J, Flocard H, Treiner J.  Hartree-Fock-Bogolyubov description of nuclei near the neutron-drip line.  Nucl. Phys. A 422 (1984) 103.  10.1016/0375-9474(84)90433-0.
  * Dobaczewski et al. (1996) Dobaczewski J, Nazarewicz W, Werner TR, Berger JF, Chinn CR, Dechargé J.  Mean-field description of ground-state properties of drip-line nuclei: Pairing and continuum effects.  Phys. Rev. C 53 (1996) 2809.  10.1103/PhysRevC.53.2809.
  * Vautherin and Brink (1972) Vautherin D, Brink DM.  Hartree-Fock Calculations with Skyrme’s Interaction. I. Spherical Nuclei.  Phys. Rev. C 5 (1972) 626.  10.1103/PhysRevC.5.626.
  * Dobaczewski and Dudek (1997) Dobaczewski J, Dudek J.  Solution of the Skyrme-Hartree-Fock equations in the Cartesian deformed harmonic oscillator basis I. The method.  Comput. Phys. Commun. 102 (1997) 166.  10.1016/S0010-4655(97)00004-0.
  * Bender et al. (2009) Bender M, Bennaceur K, Duguet T, Heenen PH, Lesinski T, Meyer J.  Tensor part of the Skyrme energy density functional. II. Deformation properties of magic and semi-magic nuclei.  Phys. Rev. C 80 (2009) 064302.  10.1103/PhysRevC.80.064302.
  * Hellemans et al. (2012) Hellemans V, Heenen PH, Bender M.  Tensor part of the Skyrme energy density functional. III. Time-odd terms at high spin.  Phys. Rev. C 85 (2012) 014326.  10.1103/PhysRevC.85.014326.
  * Ryssens et al. (2015a) Ryssens W, Hellemans V, Bender M, Heenen PH.  Solution of the Skyrme-HF+BCS equation on a 3D mesh, II: A new version of the Ev8 code.  Comput. Phys. Commun. 187 (2015a) 175.  10.1016/j.cpc.2014.10.001.
  * Valatin (1961) Valatin JG.  Generalized Hartree-Fock Method.  Phys. Rev. 122 (1961) 1012.  10.1103/PhysRev.122.1012.
  * Mang (1975) Mang HJ.  The self-consistent single-particle model in nuclear physics.  Phys. Rep. 18 (1975) 325.  10.1016/0370-1573(75)90012-5.
  * Blaizot and Ripka (1985) Blaizot JP, Ripka G.  Quantum Theory of Finite Systems (Cambridge: The MIT Press) (1985).
  * Ring and Schuck (2004) Ring P, Schuck P.  The Nuclear Many-Body Problem.  Texts and Monographs in Physics (Springer) (2004).
  * Dobaczewski and Dudek (1995) Dobaczewski J, Dudek J.  Time-odd components in the mean field of rotating superdeformed nuclei.  Phys. Rev. C 52 (1995) 1827.  10.1103/PhysRevC.52.1827.
  * Stoitsov et al. (2005) Stoitsov MV, Dobaczewski J, Nazarewicz W, Ring P.  Axially deformed solution of the Skyrme-Hartree-Fock-Bogolyubov equations using the transformed harmonic oscillator basis. The program HFBTHO (v1.66p).  Comput. Phys. Commun. 167 (2005) 43.  10.1016/j.cpc.2005.01.001.
  * Heyde and Wood (2011) Heyde K, Wood JL.  Shape coexistence in atomic nuclei.  Rev. Mod. Phys. 83 (2011) 1467.  10.1103/RevModPhys.83.1467.
  * Nakatsukasa et al. (2016) Nakatsukasa T, Matsuyanagi K, Matsuo M, Yabana K.  Time-dependent density-functional description of nuclear dynamics.  Rev. Mod. Phys. 88 (2016) 045004.  10.1103/RevModPhys.88.045004.
  * Griffin and Wheeler (1957) Griffin JJ, Wheeler JA.  Collective Motions in Nuclei by the Method of Generator Coordinates.  Phys. Rev. 108 (1957) 311.  10.1103/PhysRev.108.311.
  * Wa Wong (1975) Wa Wong C.  Generator-coordinate methods in nuclear physics.  Phys. Rep. 15 (1975) 283.  10.1016/0370-1573(75)90036-8.
  * Reinhard and Goeke (1987) Reinhard PG, Goeke K.  The generator coordinate method and quantised collective motion in nuclear systems.  Rep. Prog. Phys. 50 (1987) 1.  10.1088/0034-4885/50/1/001.
  * Verriere and Regnier (2020) Verriere M, Regnier D.  The Time-Dependent Generator Coordinate Method in Nuclear Physics.  Front. Phys. 8 (2020) 1.  10.3389/fphy.2020.00233.
  * Brink and Weiguny (1968) Brink DM, Weiguny A.  The generator coordinate theory of collective motion.  Nucl. Phys. A 120 (1968) 59.  10.1016/0375-9474(68)90059-6.
  * Onishi and Une (1975) Onishi N, Une T.  Local Gaussian Approximation in the Generator Coordinate Method.  Prog. Theor. Phys. 53 (1975) 504.  10.1143/PTP.53.504.
  * Une et al. (1976) Une T, Ikeda A, Onishi N.  Collective Hamiltonian in the Generator Coordinate Method with Local Gaussian Approximation.  Prog. Theor. Phys. 55 (1976) 498.  10.1143/PTP.55.498.
  * Bloch and Messiah (1962) Bloch C, Messiah A.  The canonical form of an antisymmetric tensor and its application to the theory of superconductivity.  Nucl. Phys. 39 (1962) 95.  10.1016/0029-5582(62)90377-2.
  * Zumino (1962) Zumino B.  Normal Forms of Complex Matrices.  J. Math. Phys. 3 (1962) 1055–1057.  10.1063/1.1724294.
  * Marević et al. (2022) Marević P, Schunck N, Ney EM, Navarro Pérez R, Verriere M, O’Neal J.  Axially-deformed solution of the Skyrme-Hartree-Fock-Bogoliubov equations using the transformed harmonic oscillator basis (iv) HFBTHO (v4.0): A new version of the program.  Comput. Phys. Commun. 276 (2022) 108367.  10.1016/j.cpc.2022.108367.
  * Drischler et al. (2020) Drischler C, Furnstahl RJ, Melendez JA, Phillips DR.  How Well Do We Know the Neutron-Matter Equation of State at the Densities Inside Neutron Stars? A Bayesian Approach with Correlated Uncertainties.  Phys. Rev. Lett. 125 (2020) 202702.  10.1103/PhysRevLett.125.202702.
  * Kravvaris et al. (2020) Kravvaris K, Quinlan KR, Quaglioni S, Wendt KA, Navrátil P.  Quantifying uncertainties in neutron-$\alpha$ scattering with chiral nucleon-nucleon and three-nucleon forces.  Phys. Rev. C 102 (2020) 024616.  10.1103/PhysRevC.102.024616.
  * Acharya and Bacca (2022) Acharya B, Bacca S.  Gaussian process error modeling for chiral effective-field-theory calculations of $np\leftrightarrow d\gamma$ at low energies.  Physics Letters B 827 (2022) 137011.  10.1016/j.physletb.2022.137011.
  * Pastore et al. (2017) Pastore A, Shelley M, Baroni S, Diget CA.  A new statistical method for the structure of the inner crust of neutron stars.  J. Phys. G: Nucl. Part. Phys. 44 (2017) 094003.  10.1088/1361-6471/aa8207.
  * Kortelainen et al. (2010) Kortelainen M, Lesinski T, Moré J, Nazarewicz W, Sarich J, Schunck N, et al.  Nuclear energy density optimization.  Phys. Rev. C 82 (2010) 024313.  10.1103/PhysRevC.82.024313.
  * Kortelainen et al. (2012) Kortelainen M, McDonnell J, Nazarewicz W, Reinhard PG, Sarich J, Schunck N, et al.  Nuclear energy density optimization: Large deformations.  Phys. Rev. C 85 (2012) 024304.  10.1103/PhysRevC.85.024304.
  * Kortelainen et al. (2014) Kortelainen M, McDonnell J, Nazarewicz W, Olsen E, Reinhard PG, Sarich J, et al.  Nuclear energy density optimization: Shell structure.  Phys. Rev. C 89 (2014) 054314.  10.1103/PhysRevC.89.054314.
  * Higdon et al. (2015) Higdon D, McDonnell JD, Schunck N, Sarich J, Wild SM.  A Bayesian approach for parameter estimation and prediction using a computationally intensive model.  J. Phys. G: Nucl. Part. Phys. 42 (2015) 034009.  10.1088/0954-3899/42/3/034009.
  * McDonnell et al. (2015) McDonnell JD, Schunck N, Higdon D, Sarich J, Wild SM, Nazarewicz W.  Uncertainty Quantification for Nuclear Density Functional Theory and Information Content of New Measurements.  Phys. Rev. Lett. 114 (2015) 122501.  10.1103/PhysRevLett.114.122501.
  * Neufcourt et al. (2018) Neufcourt L, Cao Y, Nazarewicz W, Viens F.  Bayesian approach to model-based extrapolation of nuclear observables.  Phys. Rev. C 98 (2018) 034318.  10.1103/PhysRevC.98.034318.
  * Neufcourt et al. (2020a) Neufcourt L, Cao Y, Giuliani S, Nazarewicz W, Olsen E, Tarasov OB.  Beyond the proton drip line: Bayesian analysis of proton-emitting nuclei.  Phys. Rev. C 101 (2020a) 014319.  10.1103/PhysRevC.101.014319.
  * Neufcourt et al. (2020b) Neufcourt L, Cao Y, Giuliani SA, Nazarewicz W, Olsen E, Tarasov OB.  Quantified limits of the nuclear landscape.  Phys. Rev. C 101 (2020b) 044307.  10.1103/PhysRevC.101.044307.
  * Rasmussen and Williams (2006) Rasmussen CE, Williams CKI.  Gaussian Processes for Machine Learning.  Adaptive Computation and Machine Learning (Cambridge, Mass: MIT Press) (2006).
  * Bartel et al. (1982) Bartel J, Quentin P, Brack M, Guet C, Håkansson HB.  Towards a better parametrisation of Skyrme-like effective forces: A critical study of the SkM force.  Nucl. Phys. A 386 (1982) 79.  10.1016/0375-9474(82)90403-1.
  * Schunck et al. (2014) Schunck N, Duke D, Carr H, Knoll A.  Description of induced nuclear fission with Skyrme energy functionals: Static potential energy surfaces and fission fragment properties.  Phys. Rev. C 90 (2014) 054305.  10.1103/PhysRevC.90.054305.
  * Schunck (2013a) Schunck N.  Density Functional Theory Approach to Nuclear Fission.  Acta Phys. Pol. B 44 (2013a) 263.  10.5506/APhysPolB.44.263.
  * Warda and Robledo (2011) Warda M, Robledo LM.  Microscopic description of cluster radioactivity in actinide nuclei.  Phys. Rev. C 84 (2011) 044608.  10.1103/PhysRevC.84.044608.
  * Warda et al. (2018) Warda M, Zdeb A, Robledo LM.  Cluster radioactivity in superheavy nuclei.  Phys. Rev. C 98 (2018) 041602(R).  10.1103/PhysRevC.98.041602.
  * Matheson et al. (2019) Matheson Z, Giuliani SA, Nazarewicz W, Sadhukhan J, Schunck N.  Cluster radioactivity of ${}_{118}^{294}\mathrm{Og}_{176}$.  Phys. Rev. C 99 (2019) 041304.  10.1103/PhysRevC.99.041304.
  * Dubray and Regnier (2012) Dubray N, Regnier D.  Numerical search of discontinuities in self-consistent potential energy surfaces.  Comput. Phys. Commun. 183 (2012) 2035.  10.1016/j.cpc.2012.05.001.
  * Sadhukhan (2020) Sadhukhan J.  Microscopic Theory for Spontaneous Fission.  Front. Phys. 8 (2020) 567171.  10.3389/fphy.2020.567171.
  * Bruinsma et al. (2020) Bruinsma W, Perim E, Tebbutt W, Hosking S, Solin A, Turner R.  Scalable Exact Inference in Multi-Output Gaussian Processes.  Proceedings of the 37th International Conference on Machine Learning (PMLR) (2020), Proceedings of Machine Learning Research, vol. 119, 1190.
  * Álvarez et al. (2012) Álvarez MA, Rosasco L, Lawrence ND.  Kernels for Vector-Valued Functions: A Review.  MAL 4 (2012) 195.  10.1561/2200000036.
  * Bender et al. (2020) Bender M, Bernard R, Bertsch G, Chiba S, Dobaczewski J, Dubray N, et al.  Future of nuclear fission theory.  J. Phys. G: Nucl. Part. Phys. 47 (2020) 113002.  10.1088/1361-6471/abab4f.
  * Schunck (2013b) Schunck N.  Microscopic description of induced fission.  J. Phys.: Conf. Ser. 436 (2013b) 012058.  10.1088/1742-6596/436/1/012058.
  * Ryssens et al. (2015b) Ryssens W, Heenen PH, Bender M.  Numerical accuracy of mean-field calculations in coordinate space.  Phys. Rev. C 92 (2015b) 064318.  10.1103/PhysRevC.92.064318.
  * Jin et al. (2017) Jin S, Bulgac A, Roche K, Wlazłowski G.  Coordinate-space solver for superfluid many-fermion systems with the shifted conjugate-orthogonal conjugate-gradient method.  Phys. Rev. C 95 (2017) 044302.  10.1103/PhysRevC.95.044302.
  * Regnier et al. (2017) Regnier D, Dubray N, Schunck N, Verrière M.  Microscopic description of fission dynamics: Toward a 3D computation of the time dependent GCM equation.  EPJ Web Conf. 146 (2017) 04043.  10.1051/epjconf/201714604043.
  * Zhao et al. (2021) Zhao J, Nikšić T, Vretenar D.  Microscopic self-consistent description of induced fission: Dynamical pairing degree of freedom.  Phys. Rev. C 104 (2021) 044612.  10.1103/PhysRevC.104.044612.
  * Lau et al. (2022) Lau NWT, Bernard RN, Simenel C.  Smoothing of one- and two-dimensional discontinuities in potential energy surfaces.  Phys. Rev. C 105 (2022) 034617.  10.1103/PhysRevC.105.034617.
  * Baldi (2012) Baldi P.  Autoencoders, unsupervised learning, and deep architectures.  Proceedings of ICML workshop on unsupervised and transfer learning (JMLR Workshop and Conference Proceedings) (2012), 37–49.
  * Burda et al. (2015) Burda Y, Grosse R, Salakhutdinov R.  Importance weighted autoencoders.  arXiv preprint arXiv:1509.00519 (2015).
  * Chen et al. (2012) Chen M, Xu Z, Weinberger K, Sha F.  Marginalized denoising autoencoders for domain adaptation.  arXiv preprint arXiv:1206.4683 (2012).
  * Gong et al. (2019) Gong D, Liu L, Le V, Saha B, Mansour MR, Venkatesh S, et al.  Memorizing normality to detect anomaly: Memory-augmented deep autoencoder for unsupervised anomaly detection.  Proceedings of the IEEE/CVF International Conference on Computer Vision (2019), 1705–1714.
  * Bengio et al. (2013) Bengio Y, Courville A, Vincent P.  Representation learning: A review and new perspectives.  IEEE transactions on pattern analysis and machine intelligence 35 (2013) 1798–1828.
  * Zhang et al. (2014) Zhang F, Du B, Zhang L.  Saliency-guided unsupervised feature learning for scene classification.  IEEE Transactions on Geoscience and Remote Sensing 53 (2014) 2175–2184.
  * Yu et al. (2017) Yu J, Hong C, Rui Y, Tao D.  Multitask autoencoder model for recovering human poses.  IEEE Transactions on Industrial Electronics 65 (2017) 5060–5068.
  * von Neuman and Wigner (1929) von Neuman J, Wigner E.  Uber merkwürdige diskrete Eigenwerte. Uber das Verhalten von Eigenwerten bei adiabatischen Prozessen.  Physikalische Zeitschrift 30 (1929) 467–470.
  * Teller (1937) Teller E.  The Crossing of Potential Surfaces.  J. Phys. Chem. 41 (1937) 109.  10.1021/j150379a010.
  * Longuet-Higgins et al. (1958) Longuet-Higgins HC, Öpik U, Pryce MHL, Sack RA.  Studies of the Jahn-Teller Effect. II. The Dynamical Problem.  Proc. R. Soc. Lond. A 244 (1958) 1.  10.1098/rspa.1958.0022.
  * Longuet-Higgins (1975) Longuet-Higgins H.  The intersection of potential energy surfaces in polyatomic molecules.  Proc. R. Soc. Lond. A 344 (1975) 147.  10.1098/rspa.1975.0095.
  * Domcke et al. (2011) Domcke W, Yarkony D, Köppel H, editors.  Conical Intersections: Theory, Computation and Experiment.  Advanced Series in Physical Chemistry: Volume 17 (World Scientific) (2011).  10.1142/7803.
  * Larson et al. (2020) Larson J, Sjöqvist E, Öhberg P.  Intersections in Physics. An Introduction to Synthetic Gauge Theories.  Lecture Notes in Physics (Springer) (2020).  10.1007/978-3-030-34882-3.
  * Haider and Gogny (1992) Haider Q, Gogny D.  Microscopic approach to the generator coordinate method with pairing correlations and density-dependent forces.  J. Phys. G: Nucl. Part. Phys. 18 (1992) 993.  10.1088/0954-3899/18/6/003.
  * Ioffe and Szegedy (2015) Ioffe S, Szegedy C.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  International conference on machine learning (PMLR) (2015), 448–456.
  * Srivastava et al. (2014) Srivastava N, Hinton G, Krizhevsky A, Sutskever I, Salakhutdinov R.  Dropout: a simple way to prevent neural networks from overfitting.  The journal of machine learning research 15 (2014) 1929–1958.
  * He et al. (2016) He K, Zhang X, Ren S, Sun J.  Deep residual learning for image recognition.  Proceedings of the IEEE conference on computer vision and pattern recognition (2016), 770–778.
  * Krizhevsky et al. (2012) Krizhevsky A, Sutskever I, Hinton GE.  Imagenet classification with deep convolutional neural networks.  Advances in neural information processing systems 25 (2012).
  * Zeiler and Fergus (2014) Zeiler MD, Fergus R.  Visualizing and understanding convolutional networks.  European conference on computer vision (Springer) (2014), 818–833.
  * Sermanet et al. (2013) Sermanet P, Eigen D, Zhang X, Mathieu M, Fergus R, LeCun Y.  Overfeat: Integrated recognition, localization and detection using convolutional networks.  arXiv preprint arXiv:1312.6229 (2013).