# Neural Network Emulation of Spontaneous Fission

Daniel Lay  Department of Physics and Astronomy and FRIB Laboratory, Michigan
State University, East Lansing, Michigan 48824, USA  Eric Flynn  Department of
Physics and Astronomy and FRIB Laboratory, Michigan State University, East
Lansing, Michigan 48824, USA  Samuel A. Giuliani  Departamento de Física
Teórica and CIAFF, Universidad Autónoma de Madrid, Madrid 28049, Spain
Department of Physics, Faculty of Engineering and Physical Sciences,
University of Surrey, Guildford, Surrey GU2 7XH, United Kingdom.  Witold
Nazarewicz  Department of Physics and Astronomy and FRIB Laboratory, Michigan
State University, East Lansing, Michigan 48824, USA  Leó Neufcourt  FRIB
Laboratory, Michigan State University, East Lansing, Michigan 48824, USA

###### Abstract

Background

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

Purpose

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

Methods

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

Results

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

Conclusions

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

##  I Introduction

Large scale calculations of fission properties are an essential ingredient for
the modelling of the rapid neutron-capture process ($r$ process), responsible
for the production of roughly half of the nuclei heavier than iron found in
nature Burbidge  _et al._ (1957); Cameron (1957). Fission determines the range
of the heaviest nuclei that can be synthesized during the $r$-process,
recycles the material during the neutron irradiation phase, and shapes the
final abundances Horowitz  _et al._ (2019); Kajino  _et al._ (2019); Cowan
_et al._ (2021). Given the large amount of energy released in this decay, the
presence of fissioning nuclei can leave fingerprints in the electromagnetic
counterpart produced in neutron star mergers Zhu  _et al._ (2018); Wu  _et
al._ (2019). However, as most of the fissioning nuclei produced during the $r$
process cannot be measured, theoretical predictions are indispensable to
perform accurate nuclear reaction network calculations.

During the last decades, several efforts have been devoted to the systematic
estimation of fission barriers Aboussir  _et al._ (1992); Myers and Swiatecki
(1999); Goriely  _et al._ (2009); Möller  _et al._ (2015); Agbemava  _et al._
(2017); Giuliani  _et al._ (2018), spontaneous fission half-lives Staszczak
_et al._ (2013); Baran  _et al._ (2015); Giuliani  _et al._ (2018), and
fragment distributions Kelic  _et al._ (2009); Schmidt  _et al._ (2016);
Mumpower  _et al._ (2020); Lemaître  _et al._ (2021); Sadhukhan  _et al._
(2022) of $r$-process nuclei. However, due to the inherent complexities
characterizing the theoretical description of the fission process Bender  _et
al._ (2020), most of the available calculations resort to phenomenological
approaches based on simplified assumptions. This limitation can be overcome by
employing the nuclear DFT Bender  _et al._ (2003); Duguet (2014); Schunck
(2019), which is the quantum many-body method based on effective nucleon-
nucleon interactions applicable across the whole nuclear landscape. But given
its computational costs, using DFT for fission is a daunting task for large-
scale studies of $r$-process nuclei Schunck and Robledo (2016); Bender  _et
al._ (2020); Schunck and Regnier (2022). As such, the usage of DFT emulators
can be an invaluable tool to extend the current reach of microscopic fission
calculations.

Machine learning has been used with great success in many areas of nuclear
physics (see Boehnlein  _et al._ (2022) for a recent review on this topic). In
particular, machine learning has been used in many DFT studies to emulate
potential energy surfaces (PESs), in both quantum chemistry Manzhos and
Carrington (2021); Nagai  _et al._ (2020); Jiang  _et al._ (2020); Dral  _et
al._ (2020) and in nuclear physics Akkoyun  _et al._ (2013); Verriere  _et
al._ (2022). However, these have generally focused on emulating individual
potential energy surfaces, rather than many nuclei across a portion of the
nuclear chart (or many related chemical systems in the quantum chemistry
case). In an important study, Ref. Lasseri  _et al._ (2020) succeeded in
emulating PESs and other quantities using committees of multilayer neural
networks.

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

##  II Spontaneous fission within the nuclear density functional theory

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

$P_{\textrm{fis}}=\frac{1}{1+\exp\left(2S(L)\right)},$ (1)

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

$S(L[s])=\frac{1}{\hbar}\int_{s_{\textrm{in}}}^{s_{\textrm{out}}}\sqrt{2\mathcal{M_{\textrm{eff}}}(s)(V(s)-E_{0})}\,\,ds\,,$
(2)

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

From Eq. (2) it can be deduced that the main ingredients required for the
estimation of the SF half-lives are the effective potential energy $V$ and
collective inertia $\mathcal{M_{\textrm{eff}}}$. In this work, we compute
these quantities by employing the self-consistent mean-field method Bender
_et al._ (2003); Schunck (2019) summarized in the following.

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

$\widehat{\mathcal{H}}^{\prime}=\widehat{\mathcal{H}}_{\textup{HFB}}-\sum_{\tau=n,p}\lambda_{\tau}\widehat{N}_{\tau}-\sum_{\mu=1,2,3}\lambda_{\mu}\widehat{Q}_{\mu
0}\,.$ (3)

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

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

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

The collective inertia tensor $\mathcal{M_{\mu\nu}}$ is computed within the
Adiabatic-Time-Dependent HFB (ATDHFB) approximation using the non-perturbative
scheme Yuldashbaeva  _et al._ (1999); Baran  _et al._ (2011); Giuliani and
Robledo (2018):

$\mathcal{M_{\mu\nu}}=\frac{\hbar^{2}}{2\dot{q}_{\mu}\dot{q}_{\nu}}\sum_{\alpha\beta}\frac{F^{\mu*}_{\alpha\beta}F^{\nu}_{\alpha\beta}+F^{\mu}_{\alpha\beta}F^{\nu*}_{\alpha\beta}}{E_{\alpha}+E_{\beta}}\,,$
(5)

where $q_{i}$ are the collective coordinates and

$\frac{F^{\mu}}{\dot{q}_{\mu}}=A^{\dagger}\frac{\partial\rho}{\partial
q_{\mu}}B^{*}+A^{\dagger}\frac{\partial\kappa}{\partial
q_{\mu}}A^{*}-B^{\dagger}\frac{\partial\rho^{*}}{\partial
q_{\mu}}A^{*}-B^{\dagger}\frac{\partial\kappa^{*}}{\partial q_{\mu}}B^{*}$ (6)

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

$\mathcal{M_{\textrm{eff}}}=\sum_{\mu\nu}\mathcal{M_{\mu\nu}}\frac{dq_{\mu}}{ds}\frac{dq_{\nu}}{ds}\,.$
(7)

It is important to remark that the $\mathcal{M_{\mu\nu}}$ components can
suffer from rapid oscillations in the presence of single-particle level
crossings near the Fermi surface. Such abrupt changes of occupied single-
particle configurations produce variations in the derivatives of the densities
in Eq. (6), resulting in pronounced peaks of $\mathcal{M_{\textrm{eff}}}$
along the fission path Baran  _et al._ (2011); Sadhukhan  _et al._ (2013);
Giuliani and Robledo (2018).

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

##  III Neural networks

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

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

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

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

$\displaystyle\mathcal{M}=U\Sigma U^{T}$ (8)

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

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

##  IV Neural Network Quality

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

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

###  IV.1 Potential energy surfaces

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

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

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

Figure 1 shows $\Delta V(A,Z)$ across the region of the nuclear chart
considered, for the deepest NN (7 hidden layers, with
200-175-150-125-100-75-50 hidden units), with rescaled inputs. As can be seen,
for most nuclei, $\Delta V(A,Z)\lesssim 0.5$ MeV. Exceptions occur, with most
remaining below 1.5 MeV.

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

For some nuclei, such as 308Cf, 314Fm, and 318No, relatively poor performance
may be expected: these nuclei are on the outer edge of the region of the
nuclear chart considered, and hence the NN is extrapolating from the training
region to reach them. For other nuclei, such as 232Th and 280Cm, poor
performance is unexpected: these nuclei are surrounded by training nuclei, and
so should be emulated fairly well. As such, it seems unlikely that poor
performance is due solely to the location of the nucleus on the nuclear chart
relative to the training data.

Figure 2: The reference PES (left) and the NN reconstruction (right) for
280Cm, in MeV. The ground state is marked with a
$\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\bm{\times}$
symbol.

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

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

A general improvement is observed when normalizing the inputs
$(A,Z,Q_{20},Q_{30})$ to be between 0 and 1. This is due to two factors.
First, the weights are not scale-invariant. An input much larger than 1 is
equivalent to a large initial weight, with a normalized input. Because the
final NN weights are expected to be small (hence initializations following
e.g. the Xavier initialization Glorot and Bengio (2010), as in this work), the
initial weights are far from the final values, and convergence slows. Second,
the optimization method itself is not scale-invariant: non-normalized inputs
correspond to an ill-conditioned Hessian matrix, in which case gradient
descent (and related methods) converge slowly Curry (1944); Rigler  _et al._
(1991).

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

Figure 3: The RMSE for a variety of different NN sizes. The dashed line shows
the same depth NN, but with input variables normalized to the range $[0,1]$.

###  IV.2 Collective Inertia

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

Figure 4: The reference components of $\mathcal{M}$, plotted against the NN
reconstructions, for all nuclei considered. The black line is the diagonal,
$\mathcal{M}_{\mu\nu}^{\textrm{DFT}}=\mathcal{M}_{\mu\nu}^{\textrm{NN}}$. The
blue squares/green circles/orange triangles correspond to the
training/combining/validation datasets. Note that, for use on a log scale, we
plot the absolute value of $\mathcal{M}_{23}$ (the other components are
nonnegative). $\mathcal{M}_{22}$ is in MeV-1b-2, $\mathcal{M}_{23}$ is in
MeV-1b-5/2, and $\mathcal{M}_{33}$ is in MeV-1b-3.

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

When changing the depth of the NN, performance is similar. For shallow
networks, predictions on the training dataset show a larger bias: the
distribution of points on the inertia plot is less aligned with the diagonal
for the $\mathcal{M}_{22}$ and $\mathcal{M}_{33}$ components. In other words,
the larger reference values are underestimated, and the smaller reference
values are overestimated. The validation dataset is aligned similar to the
deepest network, shown in Fig. 4. As the depth of the network is increased,
the training data points are aligned closer with the diagonal. This is
indicative of the NN tending to overfit on the training data as the number of
variational parameters increases. The distribution of $\mathcal{M}_{23}$
values remains approximately the same when increasing NN depth, with a slight
improvement on the truncated $\mathcal{M}_{23}$ values. In general, the NN
performance on the validation dataset is mostly stable when varying the NN
depth. The overarching question is whether this performance is sufficient for
predicting observable quantities of interest. As with the PES, this question
can be directly answered by looking at NN predictions of physical observables.

##  V Impact on Observable Quantities

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

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

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

###  V.1 Exit Points

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

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

Figure 5: The $Q_{30}$ component of the reconstructed lifetime-weighted exit
point, minus the reference $Q_{30}$ component, in b3/2. These results were
computed using configuration (iii), i.e. the NN was used to reconstruct both
the PES and the collective inertia.

The agreement is good between the reference exit point and the NN
reconstruction: at $\pm 1$ b3/2, we expect the fragment yields to agree well
(within the hybrid method of Refs. Sadhukhan  _et al._ (2020, 2022)). This
agreement is mainly due to the accurate PES reconstruction, as previous
studies have shown that the exit point location is fairly robust with respect
to variations in the collective inertia Sadhukhan  _et al._ (2013, 2014);
Matheson  _et al._ (2019); Sadhukhan  _et al._ (2020). This agreement holds
even for nuclei whose PES reconstruction has a large error, such as 280Cm,
indicating that the qualitative features shown in Fig. 2 are reconstructed
well enough to describe multimodality in SF.

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

Nevertheless, the error is due to the inaccurate collective inertia
reconstruction. However, it is not a systematic error. Rather, random error is
present for every deformation considered, and it is the accumulation of this
random error that causes the discrepancy. While the location of any individual
exit point is not sensitive to the collective inertia, the probability of
tunneling to a particular point depends on the probability given in Eq. (1).
Because the probability is exponentially dependent on the action (and
therefore exponentially dependent on the collective inertia reconstruction),
comparatively small errors can add up and actually switch the dominant exit
point, from asymmetric to symmetric and vice versa. This is especially
important for nuclei with a wide fission barrier, as the cumulative error
along the path is large.

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

###  V.2 Spontaneous fission half-lives

In this section we examine the performance of the NN when predicting the SF
half-life, $t_{1/2}^{\textrm{sf}}$. For the sake of simplicity, we do not
include triaxiality and pairing correlations as collective degrees of freedom,
despite their large impact on the predicted $t_{1/2}^{\textrm{sf}}$ Giuliani
and Robledo (2013); Rodríguez-Guzmán and Robledo (2014); Sadhukhan  _et al._
(2014); Giuliani  _et al._ (2014); Zhao  _et al._ (2016).

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

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

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

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

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

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

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

##  VI Conclusions

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

In future studies we plan to include more collective coordinates. Dynamic
pairing fluctuations have been shown to be important for accurately predicting
SF half-lives Giuliani and Robledo (2013); Rodríguez-Guzmán and Robledo
(2014); Sadhukhan  _et al._ (2014); Giuliani  _et al._ (2014); Zhao  _et al._
(2016), and other multipole moments have a significant impact on both the
collective motion from the ground state to the fission isomer Sheikh  _et al._
(2009); Stazczak  _et al._ (2011) and the fragment yields Warda  _et al._
(2002); Younes and Gogny (2009); Zdeb  _et al._ (2021). The NNs described in
this study are expected to have similar performance when including more
collective coordinates. However, these NNs require a large amount of training
data, and generating such data may be prohibitively expensive for many
collective coordinates. One attractive option, as demonstrated in Ref. Lasseri
_et al._ (2020), is to select nuclei for use in training in an iterative
manner based on the performance of the NN, so as to reduce the required number
of DFT calculations.

Another option is to consider different emulators that may require a smaller
amount of training data. Reduced order modeling techniques have recently been
used in nuclear physics to study both non-relativistic and relativistic DFT
Bonilla  _et al._ (2022); Giuliani  _et al._ (2023), especially for use in
uncertainty quantification. We intend to explore the usefulness of these
techniques in emulating DFT results.

##  VII Acknowledgements

This work was supported by the U.S. Department of Energy under Award Numbers
DOE-DE-NA0004074 (NNSA, the Stewardship Science Academic Alliances program),
DE-SC0013365 (Office of Science), and DE-SC0023175 (Office of Science, NUCLEI
SciDAC-5 collaboration); and by the Spanish Agencia Estatal de Investigación
(AEI) of the Ministry of Science and Innovation (MCIN) under grant agreements
No. PID2021-127890NB-I00 and No. RYC2021-031880-I funded by
MCIN/AEI/10.13039/501100011033 and the European Union-“NextGenerationEU”.

## References

  * Burbidge  _et al._ (1957) E. M. Burbidge, G. R. Burbidge, W. A. Fowler, and F. Hoyle, “Synthesis of the Elements in Stars,” Rev. Mod. Phys. 29, 547–650 (1957).
  * Cameron (1957) A. G. W. Cameron, “Nuclear Reactions in Stars and Nucleogenesis,” Publ. Astron. Soc. Pacific 69, 201 (1957).
  * Horowitz  _et al._ (2019) C. J. Horowitz, A. Arcones, B. Côté, I. Dillmann, W. Nazarewicz, I. U. Roederer, H. Schatz, A. Aprahamian, D. Atanasov, A. Bauswein, T. C. Beers, J. Bliss, M. Brodeur, J. A. Clark, A. Frebel, F. Foucart, C. J. Hansen, O. Just, A. Kankainen, G. C. McLaughlin, J. M. Kelly, S. N. Liddick, D. M. Lee, J. Lippuner, D. Martin, J. D. J. Mendoza-Temis, B. D. Metzger, M. R. Mumpower, G. Perdikakis, J. Pereira, B. W. O’Shea, R. Reifarth, A. M. Rogers, D. M. Siegel, A. Spyrou, R. Surman, X. Tang, T. Uesaka, and M. Wang, “r-process nucleosynthesis: connecting rare-isotope beam facilities with the cosmos,” J. Phys. G Nucl. Part. Phys. 46, 083001 (2019), 1805.04637 .
  * Kajino  _et al._ (2019) T. Kajino, W. Aoki, A. Balantekin, R. Diehl, M. Famiano, and G. Mathews, “Current status of r-process nucleosynthesis,” Prog. Part. Nucl. Phys. 107, 109–166 (2019).
  * Cowan  _et al._ (2021) J. J. Cowan, C. Sneden, J. E. Lawler, A. Aprahamian, M. Wiescher, K. Langanke, G. Martínez-Pinedo, and F.-K. Thielemann, “Origin of the heaviest elements: The rapid neutron-capture process,” Rev. Mod. Phys. 93, 015002 (2021), 1901.01410 .
  * Zhu  _et al._ (2018) Y.-L. Zhu, R. T. Wollaeger, N. Vassh, R. Surman, T. M. Sprouse, M. R. Mumpower, P. Möller, G. C. McLaughlin, O. Korobkin, T. Kawano, P. J. Jaffke, E. M. Holmbeck, C. L. Fryer, W. P. Even, A. J. Couture, and J. Barnes, “Californium-254 and Kilonova Light Curves,” Astrophys. J. 863, L23 (2018), 1806.09724 .
  * Wu  _et al._ (2019) M.-R. Wu, J. Barnes, G. Martínez-Pinedo, and B. D. Metzger, “Fingerprints of Heavy-Element Nucleosynthesis in the Late-Time Lightcurves of Kilonovae,” Phys. Rev. Lett. 122, 062701 (2019).
  * Aboussir  _et al._ (1992) Y. Aboussir, J. Pearson, A. Dutta, and F. Tondeur, “Thomas-Fermi approach to nuclear-mass formula,” Nucl. Phys. A 549, 155–179 (1992).
  * Myers and Swiatecki (1999) W. D. Myers and W. J. Swiatecki, “Thomas-Fermi fission barriers,” Phys. Rev. C 60, 014606 (1999).
  * Goriely  _et al._ (2009) S. Goriely, S. Hilaire, A. J. Koning, M. Sin, and R. Capote, “Towards a prediction of fission cross sections on the basis of microscopic nuclear inputs,” Phys. Rev. C 79, 024612 (2009).
  * Möller  _et al._ (2015) P. Möller, A. J. Sierk, T. Ichikawa, A. Iwamoto, and M. R. Mumpower, “Fission barriers at the end of the chart of the nuclides,” Phys. Rev. C 91, 024310 (2015).
  * Agbemava  _et al._ (2017) S. Agbemava, A. V. Afanasjev, D. Ray, and P. Ring, “Assessing theoretical uncertainties in fission barriers of superheavy nuclei,” Phys. Rev. C 95, 054324 (2017).
  * Giuliani  _et al._ (2018) S. A. Giuliani, G. Martínez-Pinedo, and L. M. Robledo, “Fission properties of superheavy nuclei for r-process calculations,” Phys. Rev. C 97, 034323 (2018), 1704.00554 .
  * Staszczak  _et al._ (2013) A. Staszczak, A. Baran, and W. Nazarewicz, “Spontaneous fission modes and lifetimes of superheavy elements in the nuclear density functional theory,” Phys. Rev. C 87, 024320 (2013).
  * Baran  _et al._ (2015) A. Baran, M. Kowal, P.-G. Reinhard, L. M. Robledo, A. Staszczak, and M. Warda, “Fission barriers and probabilities of spontaneous fission for elements with Z$\geq$100,” Nucl. Phys. A 944, 442–470 (2015).
  * Kelic  _et al._ (2009) A. Kelic, M. V. Ricciardi, and K.-H. Schmidt, “ABLA07 - towards a complete description of the decay channels of a nuclear system from spontaneous fission to multifragmentation,” arXiv:0906.4193 (2009).
  * Schmidt  _et al._ (2016) K.-H. Schmidt, B. Jurado, C. Amouroux, and C. Schmitt, “General description of fission observables: GEF model code,” Nucl. Data Sheets 131, 107–221 (2016).
  * Mumpower  _et al._ (2020) M. R. Mumpower, P. Jaffke, M. Verriere, and J. Randrup, “Primary fission fragment mass yields across the chart of nuclides,” Phys. Rev. C 101, 054607 (2020).
  * Lemaître  _et al._ (2021) J.-F. Lemaître, S. Goriely, A. Bauswein, and H.-T. Janka, “Fission fragment distributions and their impact on the $r$-process nucleosynthesis in neutron star mergers,” Phys. Rev. C 103, 025806 (2021).
  * Sadhukhan  _et al._ (2022) J. Sadhukhan, S. A. Giuliani, and W. Nazarewicz, “Theoretical description of fission yields: Toward a fast and efficient global model,” Phys. Rev. C 105, 014619 (2022).
  * Bender  _et al._ (2020) M. Bender, R. Bernard, G. Bertsch, S. Chiba, J. Dobaczewski, N. Dubray, S. A. Giuliani, K. Hagino, D. Lacroix, Z. Li, P. Magierski, J. Maruhn, W. Nazarewicz, J. Pei, S. Péru, N. Pillet, J. Randrup, D. Regnier, P.-G. Reinhard, L. M. Robledo, W. Ryssens, J. Sadhukhan, G. Scamps, N. Schunck, C. Simenel, J. Skalski, I. Stetcu, P. Stevenson, S. Umar, M. Verriere, D. Vretenar, M. Warda, and S. Åberg, “Future of nuclear fission theory,” J. Phys. G 47, 113002 (2020).
  * Bender  _et al._ (2003) M. Bender, P.-H. Heenen, and P.-G. Reinhard, “Self-consistent mean-field models for nuclear structure,” Rev. Mod. Phys. 75, 121–180 (2003).
  * Duguet (2014) T. Duguet, “The Nuclear Energy Density Functional Formalism,” in  _Euroschool Exot. Beams, Vol. IV_, Vol. 879, edited by C. Scheidenberger and M. Pfützner (Springer Berlin Heidelberg, 2014) Chap. 7, pp. 293–350.
  * Schunck (2019) N. Schunck,  _Energy Density Functional Methods for Atomic Nuclei_ (IOP Publishing, 2019).
  * Schunck and Robledo (2016) N. Schunck and L. M. Robledo, “Microscopic theory of nuclear fission: a review,” Rep. Prog. Phys. 79, 116301 (2016).
  * Schunck and Regnier (2022) N. Schunck and D. Regnier, “Theory of nuclear fission,” Prog. Part. Nucl. Phys. 125, 103963 (2022).
  * Boehnlein  _et al._ (2022) A. Boehnlein, M. Diefenthaler, N. Sato, M. Schram, V. Ziegler, C. Fanelli, M. Hjorth-Jensen, T. Horn, M. P. Kuchera, D. Lee, W. Nazarewicz, P. Ostroumov, K. Orginos, A. Poon, X.-N. Wang, A. Scheinker, M. S. Smith, and L.-G. Pang, “Colloquium: Machine learning in nuclear physics,” Rev. Mod. Phys. 94, 031003 (2022).
  * Manzhos and Carrington (2021) S. Manzhos and T. J. Carrington, “Neural network potential energy surfaces for small molecules and reactions,” Chem. Rev. 121, 10187–10217 (2021).
  * Nagai  _et al._ (2020) R. Nagai, R. Akashi, and O. Sugino, “Completing density functional theory by machine learning hidden messages from molecules,” npj Comput. Mater. 6 (2020), 10.1038/s41524-020-0310-0.
  * Jiang  _et al._ (2020) B. Jiang, J. Li, and H. Guo, “High-fidelity potential energy surfaces for gas-phase and gas–surface scattering processes from machine learning,” J. Phys. Chem. Lett. 11, 5120–5131 (2020).
  * Dral  _et al._ (2020) P. O. Dral, A. Owens, A. Dral, and G. Csányi, “Hierarchical machine learning of potential energy surfaces,” J. Chem. Phys. 152, 204110 (2020).
  * Akkoyun  _et al._ (2013) S. Akkoyun, T. Bayram, S. O. Kara, and N. Yildiz, “Consistent empirical physical formulas for potential energy curves of 38-66Ti isotopes by using neural networks,” Phys. Part. Nucl. 10, 528–534 (2013).
  * Verriere  _et al._ (2022) M. Verriere, N. Schunck, I. Kim, P. Marević, K. Quinlan, M. N. Ngo, D. Regnier, and R. D. Lasseri, “Building surrogate models of nuclear density functional theory with gaussian processes and autoencoders,” Front. in Phys. 10 (2022), 10.3389/fphy.2022.1028370.
  * Lasseri  _et al._ (2020) R.-D. Lasseri, D. Regnier, J.-P. Ebran, and A. Penon, “Taming nuclear complexity with a committee of multilayer neural networks,” Phys. Rev. Lett. 124, 162502 (2020).
  * Brack  _et al._ (1972) M. Brack, J. Damgaard, A. S. Jensen, H. C. Pauli, V. M. Strutinsky, and C. Y. Wong, “Funny Hills: The Shell-Correction Approach to Nuclear Shell Effects and Its Applications to the Fission Process,” Rev. Mod. Phys. 44, 320–405 (1972).
  * Flynn  _et al._ (2022) E. Flynn, D. Lay, S. Agbemava, P. Giuliani, K. Godbey, W. Nazarewicz, and J. Sadhukhan, “Nudged elastic band approach to nuclear fission pathways,” Phys. Rev. C 105, 054302 (2022).
  * Berger  _et al._ (1984) J. Berger, M. Girod, and D. Gogny, “Microscopic analysis of collective dynamics in low energy fission,” Nucl. Phys. A 428, 23–36 (1984).
  * Robledo  _et al._ (2019) L. M. Robledo, T. R. Rodríguez, and R. R. Rodríguez-Guzmán, “Mean field and beyond description of nuclear structure with the Gogny force: a review,” J. Phys. G: Nucl. Part. Phys. 46, 013001 (2019), 1807.02518 .
  * Egido and Robledo (2004) J. Egido and L. Robledo, “10 angular momentum projection and quadrupole correlations effects in atomic nuclei,” in  _Extended Density Functionals in Nuclear Structure Physics_, edited by G. A. Lalazissis, P. Ring, and D. Vretenar (Springer Berlin Heidelberg, Berlin, Heidelberg, 2004) pp. 269–302.
  * Robledo and Bertsch (2011) L. M. Robledo and G. F. Bertsch, “Application of the gradient method to Hartree-Fock-Bogoliubov theory,” Phys. Rev. C 84, 014312 (2011), 1104.5453 .
  * Warda  _et al._ (2002) M. Warda, J. L. Egido, L. M. Robledo, and K. Pomorski, “Self-consistent calculations of fission barriers in the Fm region,” Phys. Rev. C 66, 014310 (2002).
  * Yuldashbaeva  _et al._ (1999) E. Yuldashbaeva, J. Libert, P. Quentin, and M. Girod, “Mass parameters for large amplitude collective motion: A perturbative microscopic approach,” Phys. Lett. B 461, 1–8 (1999).
  * Baran  _et al._ (2011) A. Baran, J. A. Sheikh, J. Dobaczewski, W. Nazarewicz, and A. Staszczak, “Quadrupole collective inertia in nuclear fission: Cranking approximation,” Phys. Rev. C 84, 054321 (2011).
  * Giuliani and Robledo (2018) S. A. Giuliani and L. M. Robledo, “Non-perturbative collective inertias for fission: A comparative study,” Phys. Lett. B 787, 134–140 (2018).
  * Sadhukhan  _et al._ (2013) J. Sadhukhan, K. Mazurek, A. Baran, J. Dobaczewski, W. Nazarewicz, and J. A. Sheikh, “Spontaneous fission lifetimes from the minimization of self-consistent collective action,” Phys. Rev. C 88, 064314 (2013).
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine learning in Python,” J. Mach. Learn. Res. 12, 2825–2830 (2011).
  * Glorot and Bengio (2010) X. Glorot and Y. Bengio, “Understanding the difficulty of training deep feedforward neural networks,” in  _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_, Proceedings of Machine Learning Research, Vol. 9, edited by Y. W. Teh and M. Titterington (PMLR, Chia Laguna Resort, Sardinia, Italy, 2010) pp. 249–256.
  * Curry (1944) H. B. Curry, “The method of steepest descent for non-linear minimization problems,” Q. of Appl. Math. 2, 258–261 (1944).
  * Rigler  _et al._ (1991) A. Rigler, J. Irvine, and T. Vogl, “Rescaling of variables in back propagation learning,” Neural Netw. 4, 225–229 (1991).
  * Sadhukhan  _et al._ (2020) J. Sadhukhan, S. A. Giuliani, Z. Matheson, and W. Nazarewicz, “Efficient method for estimation of fission fragment yields of $r$-process nuclei,” Phys. Rev. C 101, 065803 (2020).
  * Sadhukhan  _et al._ (2014) J. Sadhukhan, J. Dobaczewski, W. Nazarewicz, J. A. Sheikh, and A. Baran, “Pairing-induced speedup of nuclear spontaneous fission,” Phys. Rev. C 90, 061304 (2014).
  * Matheson  _et al._ (2019) Z. Matheson, S. A. Giuliani, W. Nazarewicz, J. Sadhukhan, and N. Schunck, “Cluster radioactivity of 294Og,” Phys. Rev. C 99, 041304 (2019).
  * Giuliani and Robledo (2013) S. A. Giuliani and L. M. Robledo, “Fission properties of the Barcelona-Catania-Paris-Madrid energy density functional,” Phys. Rev. C 88, 054325 (2013).
  * Rodríguez-Guzmán and Robledo (2014) R. R. Rodríguez-Guzmán and L. M. Robledo, “Microscopic description of fission in neutron-rich plutonium isotopes with the Gogny-D1M energy density functional,” Eur. Phys. J. A 50, 142 (2014).
  * Giuliani  _et al._ (2014) S. A. Giuliani, L. M. Robledo, and R. R. Rodríguez-Guzmán, “Dynamic versus static fission paths with realistic interactions,” Phys. Rev. C 90, 054311 (2014).
  * Zhao  _et al._ (2016) J. Zhao, B.-N. Lu, T. Nikšić, D. Vretenar, and S.-G. Zhou, “Multidimensionally-constrained relativistic mean-field study of spontaneous fission: Coupling between shape and pairing degrees of freedom,” Phys. Rev. C 93, 044315 (2016).
  * Sheikh  _et al._ (2009) J. A. Sheikh, W. Nazarewicz, and J. C. Pei, “Systematic study of fission barriers of excited superheavy nuclei,” Phys. Rev. C 80, 011302 (2009).
  * Stazczak  _et al._ (2011) A. Stazczak, A. Baran, and W. Nazarewicz, “Breaking of axial and reflection symmetries in spontaneous fission of fermium isotopes,” Int. J. Mod. Phys. E 20, 552–556 (2011).
  * Younes and Gogny (2009) W. Younes and D. Gogny, “Microscopic calculation of ${}^{240}\mathrm{Pu}$ scission with a finite-range effective force,” Phys. Rev. C 80, 054313 (2009).
  * Zdeb  _et al._ (2021) A. Zdeb, M. Warda, and L. M. Robledo, “Description of the multidimensional potential-energy surface in fission of ${}^{252}\mathrm{Cf}$ and ${}^{258}\mathrm{No}$,” Phys. Rev. C 104, 014610 (2021).
  * Bonilla  _et al._ (2022) E. Bonilla, P. Giuliani, K. Godbey, and D. Lee, “Training and projecting: A reduced basis method emulator for many-body physics,” Phys. Rev. C 106, 054322 (2022).
  * Giuliani  _et al._ (2023) P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, “Bayes goes fast: Uncertainty quantification for a covariant energy density functional emulated by the reduced basis method,” Front. in Phys. 10 (2023), 10.3389/fphy.2022.1054524.


# An efficient, concatenated, bosonic code for additive Gaussian noise

Kosuke Fukui  Department of Applied Physics, School of Engineering, The
University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan  Nicolas C. Menicucci  Centre
for Quantum Computation & Communication Technology, School of Science, RMIT
University, Melbourne, VIC 3000, Australia

###### Abstract

Bosonic codes offer noise resilience for quantum information processing. A
common type of noise in this setting is additive Gaussian noise, and a long-
standing open problem is to design a concatenated code that achieves the
hashing bound for this noise channel. Here we achieve this goal using a
Gottesman-Kitaev-Preskill (GKP) code to detect and discard error-prone qubits,
concatenated with a quantum parity code to handle the residual errors. Our
method employs a linear-time decoder and has applications in a wide range of
quantum computation and communication scenarios.

Introduction.—Quantum-error-correcting codes known as _bosonic codes_ [1]
protect discrete quantum information encoded in one or more bosonic modes. The
infinite-dimensional nature of the bosonic Hilbert space allows for more
sophisticated encodings than the more common single-photon encoding [2, 3, 4]
or that of a material qubit like a transmon [5, 6, 7]. The variety of codes
allows for multiple ways to protect the encoded quantum information against
decoherence at the physical level, with better than break-even performance
demonstrated recently [8, 9]. Based on the ubiquitous quantum harmonic
oscillator, these versatile codes find applications in optical, solid-state,
and vibrational systems. We call qubits encoded in a bosonic code _bosonic
qubits_.

While quantum supremacy has been demonstrated in both solid-state qubits [10]
and optics [11], the ultimate goal of a large-scale, fault-tolerant quantum
computer will require additional innovations, and its ultimate architecture
remains an open question. Such a device’s requirements [12] can be broadly
classified into scalability (many qubits) and fault tolerance (of good
quality), and architectures designed to use bosonic qubits as the information
carriers have recently demonstrated prominent advances in both areas.

Progress on scalability has been most significant in optics through
demonstrations of computationally universal continuous-variable (CV) cluster
states [13, 14] comprising ${\sim 10^{4}}$ modes [15, 16] and measurement-
based implementation of CV quantum gates [17, 18]. When used to process the
bosonic qubits proposed by Gottesman, Kitaev, and Preskill (GKP) [19]—and with
high enough squeezing—these architectures can be made fault tolerant [20, 21].

The GKP qubit [19] has emerged as a promising bosonic qubit for fault
tolerance due to its excellent performance against common types of noise when
compared with other known codes [1]. Experiments involving trapped ions [22]
and superconducting circuits [23] have demonstrated a GKP qubit, with the
latter boasting a squeezing level close to 10 dB. This level is sufficient for
fault tolerance in some proposed architectures [24, 25] and is approaching
that required by others [26, 27]. Numerous proposals exist to produce these
states in optics, as well (see Ref. [28] and references therein). Furthermore,
the GKP qubit performs well for quantum communication [29, 30, 31] thanks to
the robustness against photon loss [1]. In fact, recent results show that
using GKP qubits may greatly enhance long-distance quantum communication [32,
33].

A common type of noise in bosonic systems is modelled by the Gaussian quantum
channel (GQC) [34, 19], also known as additive Gaussian noise. While the
landscape of Gaussian operations includes other types of channels [35, 36],
the GQC is a particularly common one [37], and it is the one that we focus on
for this work. The GQC is a simple, canonical type of noise for analyzing
bosonic code performance [19, 34]. Buoyed by the fact that displacements form
an operator basis, protecting against the GQC allows some level of protection
against all types of bosonic noise [19]. In fact, the GKP encoding is
specifically designed to protect against the GQC (and thus against bosonic
noise in general), but despite this, its performance “out of the box” as a
single-mode code is suboptimal against the GQC [34, 19].

A long-standing open problem in CV quantum information is to design a simple
and efficient concatenated code that achieves the hashing bound (discussed
below) of the GQC [34]. This bound can be achieved by GKP-type codes based on
high-dimensional sphere packing [34], but the authors of that work were
unsatisfied with this because such a code is not concatenated and thus offers
no obvious structure that might be exploited [38, 39] to further improve its
performance.

Analog quantum error correction (QEC) [40] makes strides toward achieving the
goal of Ref. [34] by using the real-valued syndrome of a GKP qubit to improve
error recovery in a concatenated code. Likelihoods of the error patterns
obtained from the error syndrome are compared, and the most likely error
pattern is selected. In fact, when used with a suitable qubit code, analog QEC
can achieve the hashing bound of the GQC [40, 24]. This would seem to be the
end of the story, except for one major drawback: The decoder for analog QEC
employs a type of belief propagation [41] that becomes unwieldy as the qubit-
level code gets bigger. This is because one must model the entire probability
distribution of the multimode code to select the most likely CV error pattern.

What we would like instead is a simple CV-level decoder that generates
discrete outcomes that can be fed directly into a qubit-level code at the next
level of concatenation. This is the key innovation that makes further
improvements feasible since more complicated codes or additional layers of
concatenation do not require modifying the CV-level decoding scheme, thus
keeping the decoder simple and efficient.

In this work, we achieve this goal. Our innovation uses the CV-level
measurement outcome from GKP error correction to decide whether to keep the
qubit or discard it entirely and treat it as a located erasure error. This is
a simple, local decoding step and doesn’t require complicated modelling of CV-
level errors. The quantum parity code (QPC) [42] is well suited to dealing
with the missing qubits [43, 44, 45], and we numerically show that
concatenating the GKP code with a QPC achieves the hashing bound of the GQC
with a small code and straightforward decoding in linear time.

GKP qubit.—The GKP code encodes a qubit in an oscillator in a way that
protects against errors caused by small displacements in the $q$ (position)
and $p$ (momentum) quadratures [19]. (We use conventions $a=(q+ip)/\sqrt{2}$,
$[q,p]=i$, $\hbar=1$, vacuum variance $=1/2$.) The ideal code states of the
GKP code are Dirac combs in $q$ and in $p$. Physical states are finitely
squeezed approximations to these and are often modelled as a comb of Gaussian
peaks of width (i.e., standard deviation) $\sigma$ with separation
$\sqrt{\pi}$ modulated by a larger Gaussian envelope of width $1/\sigma$.
Since these approximate states are not orthogonal, there is a probability of
misidentifying $\ket{{0}}$ as $\ket{{1}}$ (and vice versa) in a measurement of
logical-$Z$, which is implemented by a $q$ measurement and binning to the
nearest integer multiple of $\sqrt{\pi}$. Similarly, $\ket{+}$ and $\ket{-}$
may be misidentified when measuring logical-$X$ with a $p$ measurement. A
qubit-level measurement error occurs when the measured outcome is more than
$\sqrt{\pi}/2$ away from the correct outcome.

Gaussian quantum channel.—The GKP code is tailored to combat the GQC, which
randomly displaces the state in phase space according to a Gaussian
distribution [19, 34]. The GQC is described by the superoperator
$\mathcal{G}_{\xi}$ acting on density operator $\rho$ as

$\displaystyle\rho\to\mathcal{G}_{\xi}(\rho)$
$\displaystyle=\frac{1}{\pi{\xi}^{2}}\int
d^{2}\alpha\,\mathrm{e}^{-|\alpha|^{2}/{{\xi}^{2}}}D(\alpha)\rho
D(\alpha)^{\dagger},$ (1)

where $D(\alpha)=e^{\alpha a^{\dagger}-\alpha^{*}a}$ is the phase-space
displacement operator. With $\alpha=\frac{1}{\sqrt{2}}(u+iv)$, the position
$q$ and momentum $p$ are displaced independently as $q\to q+u$, $p\to p+v$,
where $u$ and $v$ are real Gaussian random variables with mean zero and
variance $\xi^{2}$. Therefore, the GQC maintains the locations of the Gaussian
peaks in the probability for the measurement outcome, but it increases the
variance of each spike by $\xi^{2}$ in both quadratures.

Hashing bound.—The hashing bound, first introduced for qubit Pauli channels
[46, 47, 48, 49], may be generalised to the GQC by maximising the one-shot
coherent information over Gaussian states [50, 34]. The hashing bound is a
lower bound on the quantum communication rate through a channel having a set
level of noise. It is only known to be a tight lower bound for some channels
[51], and the GQC is not one of them [36]. Conversely, by fixing the rate, the
hashing bound gives a level of noise that should still allow for quantum
communication at the chosen rate using some quantum-error-correcting code.

When choosing a rate of zero, the (zero-rate) hashing bound represents a
threshold level of noise below which finite-rate quantum communication should
be possible—that is, it represents a target for a minimum error threshold when
designing codes for a specific class of channels. When a code displays an
error threshold that matches the hashing bound, we say the code “achieves the
hashing bound” for the channel. In fact, concatenated qubit codes exist that
_exceed_ the hashing bound for certain qubit channels [38, 39]—that is, they
have an error threshold that is strictly larger than that prescribed by the
hashing bound—leading to hopes that similar results may be true for
concatenated codes applied to the GQC [34]. The (zero-rate) hashing bound for
the GQC $\mathcal{G}_{\xi}$, is the standard deviation [50, 34]

$\displaystyle\xi_{\text{HB}}\coloneqq\frac{1}{\sqrt{e}}\approx 0.607\,.$ (2)

Figure 1: Effect of additive Gaussian noise [52] on meassuring a GKP qubit.
(a) Effect of shift by $u\pmod{2\sqrt{\pi}}$, distributed according to $p(u)$
[Eq. (3)], on an ordinary measurement of a GKP qubit [19]. (b) The highly
reliable measurement (HRM) [24] flags outcomes in the $2\delta$-wide “danger
zone” (yellow) as unreliable. (c) Postselected error probability of the HRM
for several values of $\delta$. (d) Corresponding success probability. Note:
$(\text{Squeezing level in
dB})=-10\log_{10}(\sigma^{2}/\sigma_{\text{vac}}^{2})$, where the vacuum
variance ${\sigma_{\text{vac}}^{2}=\tfrac{1}{2}}$.

Noise model.—GKP error correction, in both its original [19] (Steane-style
[53]) form and in its teleportation-based [54] (Knill-style [55]) form,
involve measuring the deviation of the state’s support in each quadrature
($q,p$) away from an integer multiple of $\sqrt{\pi}$. These measurement
outcomes—each of the form $s_{\rm m}=n\sqrt{\pi}+\Delta_{\rm m}$, where $n$ is
an integer and $|\Delta_{\rm m}|\leq\sqrt{\pi}/2$—together form the syndrome.
Normally, each value of $\Delta_{\rm m}$ locally determines the displacement
to apply in order to correct the error—either snapping back to grid in the
original method [19] or applying a logical Pauli in the teleportation-based
method [54]. Analog QEC [40, 24] instead feeds all these real-valued syndrome
data $s_{\rm m}$ directly to a higher-level decoder, which makes a global
decision. Our proposal keeps aspects of both approaches. We use $\Delta_{\rm
m}$ to locally decide whether to try to correct the error or to give up and
report the qubit as lost to the next-level decoder.

We model a damaged GKP codeword as an ideal one [19] that has been displaced
by a definite (but unknown) amount in each quadrature. This approximately
models the errors due to both coherent and incoherent noise [52, 20, 40, 24]
and simplifies the analysis. Due to the $2\sqrt{\pi}$-periodicity of all GKP
codewords, displacements by $u$ are equivalent to those by $u+2k\sqrt{\pi}$
for any integer $k$. Thus, given any initial distribution $p_{0}(u)$ of the
unknown displacement $u$ in a single quadrature, its effect on a GKP codeword
is captured by folding $p_{0}(u)$ into the wrapped distribution
$p(u)=\sum_{k\in\mathbb{Z}}p_{0}(u+2k\sqrt{\pi})$, whose domain is
$[-\sqrt{\pi},\sqrt{\pi})$. When $p_{0}$ is a zero-mean Gaussian of variance
$\sigma^{2}$, it wraps into

$\displaystyle
p(u)=\frac{1}{2\sqrt{\pi}}\vartheta\left(-\frac{u}{2\sqrt{\pi}},\frac{i\sigma^{2}}{2}\right),$
(3)

where ${\vartheta(z,\tau)=\sum_{m\in\mathbb{Z}}\exp\bigl{[}2\pi
i\bigl{(}\tfrac{1}{2}m^{2}\tau+mz\bigr{)}\bigr{]}}$ is a Jacobi theta function
of the third kind. Figure 1(a) shows this distribution and the logical effect
of a shift by ${u\pmod{2\sqrt{\pi}}}$ on measuring a GKP codeword.

Highly-reliable measurement.—Logical errors result when the GKP syndrome value
$s_{\rm m}$, which is wrapped$\mod\sqrt{\pi}$, misidentifies $u$ as
$u\pm\sqrt{\pi}$ [19]. The highly-reliable measurement (HRM) [24] buffers
against this possibility by introducing a “danger zone” of outcomes
$0\leq\sqrt{\pi}/2-|\Delta_{\rm m}|<\delta$ for some $\delta>0$. Outcomes in
this zone are flagged as unreliable, with $\delta\to 0$ recovering the usual
case [19]. This corresponds to flagging as unreliable any displacement
$u\pmod{2\sqrt{\pi}}$ that falls within $\delta$ of a crossover point
$\pm\sqrt{\pi}/2$, as shown in Fig. 1(b). When the HRM flags a result $s_{\rm
m}$ as unreliable, the corresponding qubit is discarded and treated as a
located erasure error ($s_{\rm m}\to E$), while otherwise the result is kept
and binned as usual ($s_{\rm m}\to\pm 1$) [19]. The HRM is thus a ternary
(3-outcome) decoder for GKP qubits that maps each raw CV outcome $s_{\rm m}$
from $\mathbb{R}\to\\{\pm 1,E\\}$.

Given a definite displacement $u\in[-\sqrt{\pi},\sqrt{\pi})$, we define
probabilities for three cases: the measurement result is correct,
$P^{(\text{c})}={\Pr(|u|<\sqrt{\pi}/2-\delta)}$; the result is incorrect,
$P^{(\text{i})}={\Pr(|u|>\sqrt{\pi}/2+\delta)}$; or the result is unreliable
and the qubit discarded,
$P^{(\text{d})}={\Pr(-\delta<|u|-\sqrt{\pi}/2<\delta)}$. We further define the
_success probability_ , $1-P^{(\text{d})}$, as the probability the qubit was
not discarded and the _postselected error probability_ , $P_{\rm
post}^{(\text{i})}=P^{(\text{i})}/(1-P^{(\text{d})})$, as the probability of
getting an incorrect outcome within the sample of qubits that are not
discarded. Decreasing the postselected error probability (by increasing
$\delta$) reduces the success probability [24], as shown in Fig. 1(c,d).

Trading errors for loss—The HRM is the key to achieving the hashing bound of
the GQC without the computational overhead of analog QEC. Analog QEC requires
modelling the joint likelihood of real-valued outcomes over multi-mode
codewords, while the HRM maps locally detected unreliable results to lost
qubits at known locations.

Loss-tolerant QEC codes were originally proposed to overcome loss of
individual photons—the main hurdle in QC based on a photonic qubit. Here,
concatenating GKP qubits with one of these codes compensates for the discarded
(“lost”) qubits due to using the HRM. This tradeoff of (unlocated) errors for
(located) erasure makes the logical qubit more robust. In the following, we
concatenate GKP qubits with the QPC proposed by Ralph et al. [42] and
implement teleportation-based QEC as proposed by Muralidharan et al. [43].

The $(n,m)$-QPC [43] is an $nm$-qubit code built from $n$ blocks of $m$
qubits. Logical basis states are $\ket{\pm}_{\rm
L}={2^{-n/2}\bigl{(}\ket{0}^{\otimes m}\pm\ket{1}^{\otimes m}\bigr{)}^{\otimes
n}}$. In our code, the physical qubit states are square-lattice GKP states
[19] of a single bosonic mode—i.e., $\ket{0}=\ket{0_{\text{GKP}}}$, and
$\ket{1}=\ket{1_{\text{GKP}}}$.

Figure 2: Trading errors for loss to improve error correction. (a)
Teleportation-based QEC using GKP qubits [19] concatenated with the quantum
parity code (QPC) [42, 43]. All gates are shown at the logical level. The
feedforward operations $D_{p}$ and $D_{q}$ implement logical-level $X$ and $Z$
operations via physical-level displacements [19]. The outcomes of the logical
Bell measurement decide what feedforward is required [43]. (b) Encoded
measurement with the $(n,m)$-QPC in the $X$ basis ($n=5$, $m=4$ shown).
Performing the HRM in $p$ at the physical level [Fig. 1(b)] gives, for each
qubit, either a binary outcome ($\pm 1$) or a located erasure error ($E$). The
latter is indicated in red (with the unreliable binary outcome underneath) and
occurs when the CV-level outcome is in the “danger zone” (see Fig. 1) and
therefore unreliable. A horizontal block is ignored if it has any discarded
outcomes. All remaining blocks have their parity (product) taken, after which
a majority vote of those parities determines the logical outcome—with heralded
failure if there is no majority. (c) Encoded measurement with the $(5,4)$-QPC
in the $Z$ basis. This is similar to (b), but this time the HRM is done in
$q$, and majority voting within a block precedes taking the parity of the
blocks’ voting outcomes—with heralded failure if any block has no majority.
See Ref. [43] for more details. Without postselection, both logical outcomes
would have been flipped (likely due to uncorrected errors in the discarded
values).

We implement syndrome measurements by adapting the teleportation-based
protocol of Ref. [43] to this setting, as shown in Fig. 2(a). Using
teleportation to do the correction guarantees that the output state is already
in the logical subspace, and only logical corrections are required at the end
[55]. We use HRMs to implement the base-level qubit measurements, so we get
(after postselection) either a qubit value ($\pm 1$) or a located erasure
($E$) for every measured qubit. These outcomes are processed into logical
outcomes using the method of Ref. [43], which is described briefly in Fig.
2(b,c). These logical outcomes determine the logical correction to apply at
the end. In the example shown, $X_{\rm L}=-1$ and $Z_{\rm L}=-1$, so we
perform logical bit- and phase-flip operations on the output qubit to correct
it. At the physical level, these occur via displacements by integer multiples
of $\sqrt{\pi}$ [19].

Numerical simulation.—We evaluate our proposed QEC method using a Monte Carlo
simulation of the circuit in Fig. 2(a). The input state passes through a GQC
[Eq. (1)]. Our model simulates code-capacity noise (i.e., assuming no errors
aside from the channel noise itself [56, 57]) in order to evaluate the best
possible performance of our code under the GQC and to compare with previous
results [19, 34, 40, 24]. We consider two cases: (a) the conventional case,
which corresponds to choosing $\delta=0$, and (b) an optimized choice of
$\delta=0.223\sqrt{\pi}$ (explained below).

Figure 3: failure probabilities using the $(n,m)$-QPC, as shown in Fig. 2, for
(a) $\delta=0$ (conventional GKP error correction [19]) and (b) an optimized
value of $\delta=0.223\sqrt{\pi}$. The former (a) has a threshold $\xi\approx
0.555$ (identified by the crossover point), which matches previous work [19,
34]. The latter (b) has a threshold $\xi\approx 0.607$, which numerically
achieves the hashing bound of the GQC [Eq. (2)].

A logical-$X$ (-$Z$) error occurs when the decoding procedure leads to either
(1) the wrong value for $X_{\rm L}$ ($Z_{\rm L}$) or (2) heralded failure of
that measurement—see Fig. 2(b,c), caption. While one could treat the heralded
failures in a more sophisticated fashion (concatenating with a higher-level
code, for instance), in our simulations we simply treat these as logical
errors on that outcome. The probability of a logical-$X$ (-$Z$) error is
denoted $E_{\rm X}$ ($E_{\rm Z}$), and the overall logical error probability
$p_{E}$ is $1-(1-E_{\rm X})(1-E_{\rm Z})$, which we call the _failure
probability_ of the QEC.

Figure 3 shows the performance of both methods as a function of the standard
deviation of the GKP qubit for several sizes $(n,m)$ of the QPC. These sizes
are chosen to ensure approximately symmetric noise at the output (i.e.,
$E_{\rm X}\sim E_{\rm Z}$). The conventional method (a), which uses
$\delta=0$, gives a threshold of $\xi\approx 0.555$, matching previous work
with concatenated codes and simple decoding [19, 34]. For our new method (b),
we optimize the value of $\delta$ to maximize the threshold value of the
standard deviation, giving $\delta=0.223\sqrt{\pi}$, which gives a loss
probability $\sim 38\%$ at threshold. Importantly, the corresponding threshold
numerically achieves the hashing bound of the GQC: $\xi\approx 0.607$ [Cf. Eq.
(2)].

Conclusion.—The key insight of this work is that one does not need to model
the full likelihood function [40, 24] to correctly interpret GKP syndrome
information [19] within a concatenated code. Instead, the real-valued outcomes
can be coarse grained to one of three qubit-level outcomes through the HRM
mapping $\mathbb{R}\to\\{\pm 1,E\\}$, where $E$ represents an untrustworthy
value. These ternary outcomes suffice to achieve the hashing bound of the GQC
by treating $E$ outcomes as erasure errors and concatenating with a qubit-
level code designed to handle such errors [42, 43, 44, 45].

The innovation of this work over previous hashing-bound-achieving methods [40,
24] lies in the efficiency and versatility of the decoder. Respectively, (1)
decoding happens in linear time since the CV-level decoding is entirely local;
and (2) the HRM wraps each GKP qubit in a simple _error-detecting_ code, so
concatenating with any qubit-level code designed to handle erasures [58, 59,
60] should benefit from this type of outcome mapping. Further applications and
extensions include improved decoding in GKP-based architectures (e.g., [26,
27]) and in codes that exploit biased noise (e.g., [57, 61, 62]).

Acknowledgments.—We thank Joe Fitzsimons, Tim Ralph, Ben Baragiola, and
Giacomo Pantaleoni for discussions. We acknowledge the organizers of the BBQ
2019 workshop, where early results of this work were presented. This work is
supported by the Australian Research Council (ARC) Centre of Excellence for
Quantum Computation and Communication Technology (Project No. CE170100012).
K.F. acknowledges financial support from donations from Nichia Corporation.

## References

  * Albert  _et al._ [2018] V. V. Albert, K. Noh, K. Duivenvoorden, D. J. Young, R. Brierley, P. Reinhold, C. Vuillot, L. Li, C. Shen, S. Girvin, _et al._ , Performance and structure of single-mode bosonic codes, Physical Review A 97, 032346 (2018).
  * Knill  _et al._ [2001] E. Knill, R. Laflamme, and G. J. Milburn, A scheme for efficient quantum computation with linear optics, nature 409, 46 (2001).
  * Kok  _et al._ [2007] P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P. Dowling, and G. J. Milburn, Linear optical quantum computing with photonic qubits, Reviews of modern physics 79, 135 (2007).
  * Rudolph [2017] T. Rudolph, Why i am optimistic about the silicon-photonic route to quantum computing, APL Photonics 2, 030901 (2017).
  * Kjaergaard  _et al._ [2020] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, Superconducting qubits: Current state of play, Annual Review of Condensed Matter Physics 11, 369 (2020).
  * Koch  _et al._ [2007] J. Koch, M. Y. Terri, J. Gambetta, A. A. Houck, D. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Charge-insensitive qubit design derived from the cooper pair box, Physical Review A 76, 042319 (2007).
  * Schreier  _et al._ [2008] J. Schreier, A. A. Houck, J. Koch, D. I. Schuster, B. Johnson, J. Chow, J. M. Gambetta, J. Majer, L. Frunzio, M. H. Devoret, _et al._ , Suppressing charge noise decoherence in superconducting charge qubits, Physical Review B 77, 180502 (2008).
  * Ofek  _et al._ [2016] N. Ofek, A. Petrenko, R. Heeres, P. Reinhold, Z. Leghtas, B. Vlastakis, Y. Liu, L. Frunzio, S. Girvin, L. Jiang, _et al._ , Extending the lifetime of a quantum bit with error correction in superconducting circuits, Nature 536, 441 (2016).
  * Hu  _et al._ [2019] L. Hu, Y. Ma, W. Cai, X. Mu, Y. Xu, W. Wang, Y. Wu, H. Wang, Y. Song, C.-L. Zou, _et al._ , Quantum error correction and universal gate set operation on a binomial bosonic logical qubit, Nature Physics 15, 503 (2019).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Zhong  _et al._ [2020] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, _et al._ , Quantum computational advantage using photons, Science 370, 1460 (2020).
  * DiVincenzo [2000] D. P. DiVincenzo, The physical implementation of quantum computation, Fortschritte der Physik: Progress of Physics 48, 771 (2000).
  * Menicucci  _et al._ [2006] N. C. Menicucci, P. Van Loock, M. Gu, C. Weedbrook, T. C. Ralph, and M. A. Nielsen, Universal quantum computation with continuous-variable cluster states, Physical review letters 97, 110501 (2006).
  * Menicucci [2011] N. C. Menicucci, Temporal-mode continuous-variable cluster states using linear optics, Physical Review A 83, 062314 (2011).
  * Asavanant  _et al._ [2019] W. Asavanant, Y. Shiozawa, S. Yokoyama, B. Charoensombutamon, H. Emura, R. N. Alexander, S. Takeda, J.-i. Yoshikawa, N. C. Menicucci, H. Yonezawa, _et al._ , Generation of time-domain-multiplexed two-dimensional cluster state, Science 366, 373 (2019).
  * Larsen  _et al._ [2019] M. V. Larsen, X. Guo, C. R. Breum, J. S. Neergaard-Nielsen, and U. L. Andersen, Deterministic generation of a two-dimensional cluster state, Science 366, 369 (2019).
  * Asavanant  _et al._ [2020] W. Asavanant, B. Charoensombutamon, S. Yokoyama, T. Ebihara, T. Nakamura, R. N. Alexander, M. Endo, J.-i. Yoshikawa, N. C. Menicucci, H. Yonezawa, _et al._ , One-hundred step measurement-based quantum computation multiplexed in the time domain with 25 mhz clock frequency, arXiv preprint arXiv:2006.11537 (2020).
  * Larsen  _et al._ [2020] M. V. Larsen, X. Guo, C. R. Breum, J. S. Neergaard-Nielsen, and U. L. Andersen, Deterministic multi-mode gates on a scalable photonic quantum computing platform, arXiv preprint arXiv:2010.14422 (2020).
  * Gottesman  _et al._ [2001] D. Gottesman, A. Kitaev, and J. Preskill, Encoding a qubit in an oscillator, Physical Review A 64, 012310 (2001).
  * Menicucci [2014] N. C. Menicucci, Fault-tolerant measurement-based quantum computing with continuous-variable cluster states, Physical review letters 112, 120504 (2014).
  * Walshe  _et al._ [2019] B. W. Walshe, L. J. Mensen, B. Q. Baragiola, and N. C. Menicucci, Robust fault tolerance for continuous-variable cluster states with excess antisqueezing, Physical Review A 100, 010301 (2019).
  * Flühmann  _et al._ [2019] C. Flühmann, T. L. Nguyen, M. Marinelli, V. Negnevitsky, K. Mehta, and J. P. Home, Encoding a qubit in a trapped-ion mechanical oscillator, Nature 566, 513 (2019).
  * Campagne-Ibarcq  _et al._ [2020] P. Campagne-Ibarcq, A. Eickbusch, S. Touzard, E. Zalys-Geller, N. E. Frattini, V. V. Sivak, P. Reinhold, S. Puri, S. Shankar, R. J. Schoelkopf, _et al._ , Quantum error correction of a qubit encoded in grid states of an oscillator, Nature 584, 368 (2020).
  * Fukui  _et al._ [2018] K. Fukui, A. Tomita, A. Okamoto, and K. Fujii, High-threshold fault-tolerant quantum computation with analog quantum error correction, Physical review X 8, 021054 (2018).
  * Fukui [2019] K. Fukui, High-threshold fault-tolerant quantum computation with the gkp qubit and realistically noisy devices, arXiv preprint arXiv:1906.09767 (2019).
  * Bourassa  _et al._ [2020] J. E. Bourassa, R. N. Alexander, M. Vasmer, A. Patil, I. Tzitrin, T. Matsuura, D. Su, B. Q. Baragiola, S. Guha, G. Dauphinais, _et al._ , Blueprint for a scalable photonic fault-tolerant quantum computer, arXiv preprint arXiv:2010.02905 (2020).
  * Larsen  _et al._ [2021] M. V. Larsen, C. Chamberland, K. Noh, J. S. Neergaard-Nielsen, and U. L. Andersen, A fault-tolerant continuous-variable measurement-based quantum computation architecture, arXiv preprint arXiv:2101.03014 (2021).
  * Tzitrin  _et al._ [2020] I. Tzitrin, J. E. Bourassa, N. C. Menicucci, and K. K. Sabapathy, Progress towards practical qubit computation using approximate gottesman-kitaev-preskill codes, Physical Review A 101, 032315 (2020).
  * Bennett and Brassard [1984] C. H. Bennett and G. Brassard, Quantum cryptography: Public key distribution and coin tossing, in  _Proceedings of IEEE International Conference on Computers, Systems and Signal Processing_ , Vol. 175 (New York, 1984) p. 8.
  * Briegel  _et al._ [1998] H.-J. Briegel, W. Dür, J. I. Cirac, and P. Zoller, Quantum repeaters: the role of imperfect local operations in quantum communication, Physical Review Letters 81, 5932 (1998).
  * Kimble [2008] H. J. Kimble, The quantum internet, Nature 453, 1023 (2008).
  * Fukui  _et al._ [2020] K. Fukui, R. N. Alexander, and P. van Loock, All-optical long-distance quantum communication with gottesman-kitaev-preskill qubits, arXiv preprint arXiv:2011.14876 (2020).
  * Rozpędek  _et al._ [2020] F. Rozpędek, K. Noh, Q. Xu, S. Guha, and L. Jiang, Quantum repeaters based on concatenated bosonic and discrete-variable quantum codes, arXiv preprint arXiv:2011.15076 (2020).
  * Harrington and Preskill [2001] J. Harrington and J. Preskill, Achievable rates for the gaussian quantum channel, Physical Review A 64, 062301 (2001).
  * Caruso  _et al._ [2006] F. Caruso, V. Giovannetti, and A. S. Holevo, One-mode bosonic gaussian channels: a full weak-degradability classification, New Journal of Physics 8, 310 (2006).
  * Holevo [2007] A. S. Holevo, One-mode quantum gaussian channels: Structure and quantum capacity, Problems of Information Transmission 43, 1 (2007).
  * Eisert and Wolf [2007] J. Eisert and M. M. Wolf, Gaussian Quantum Channels, in  _Quantum Information with Continuous Variables of Atoms and Light_ (World Scientific, London, 2007) pp. 23–42.
  * Shor and Smolin [1996] P. W. Shor and J. A. Smolin, Quantum error-correcting codes need not completely reveal the error syndrome, arXiv preprint quant-ph/9604006 (1996).
  * DiVincenzo  _et al._ [1998] D. P. DiVincenzo, P. W. Shor, and J. A. Smolin, Quantum-channel capacity of very noisy channels, Physical Review A 57, 830 (1998).
  * Fukui  _et al._ [2017] K. Fukui, A. Tomita, and A. Okamoto, Analog quantum error correction with encoding a qubit into an oscillator, Physical review letters 119, 180507 (2017).
  * Poulin [2006] D. Poulin, Optimal and efficient decoding of concatenated quantum block codes, Physical Review A 74, 052333 (2006).
  * Ralph  _et al._ [2005] T. C. Ralph, A. Hayes, and A. Gilchrist, Loss-tolerant optical qubits, Physical review letters 95, 100501 (2005).
  * Muralidharan  _et al._ [2014] S. Muralidharan, J. Kim, N. Lütkenhaus, M. D. Lukin, and L. Jiang, Ultrafast and fault-tolerant quantum communication across long distances, Physical review letters 112, 250501 (2014).
  * Munro  _et al._ [2012] W. J. Munro, A. M. Stephens, S. J. Devitt, K. A. Harrison, and K. Nemoto, Quantum communication without the necessity of quantum memories, Nature Photonics 6, 777 (2012).
  * Ewert  _et al._ [2016] F. Ewert, M. Bergmann, and P. van Loock, Ultrafast long-distance quantum communication with static linear optics, Physical review letters 117, 210501 (2016).
  * Wilde [2017] M. Wilde,  _Quantum Information Theory_ , second edition ed. (Cambridge University Press, Cambridge, UK ; New York, 2017).
  * Lloyd [1997] S. Lloyd, Capacity of the noisy quantum channel, Physical Review A 55, 1613 (1997).
  * Shor [2002] P. Shor, The quantum channel capacity and coherent information (2002).
  * Devetak [2005] I. Devetak, The private classical capacity and quantum capacity of a quantum channel, IEEE Transactions on Information Theory 51, 44 (2005).
  * Holevo and Werner [2001] A. S. Holevo and R. F. Werner, Evaluating capacities of bosonic gaussian channels, Physical Review A 63, 032312 (2001).
  * Devetak and Shor [2005] I. Devetak and P. W. Shor, The capacity of a quantum channel for simultaneous transmission of classical and quantum information, Communications in Mathematical Physics 256, 287 (2005).
  * Glancy and Knill [2006] S. Glancy and E. Knill, Error analysis for encoding a qubit in an oscillator, Physical Review A 73, 012325 (2006).
  * Steane [1997] A. M. Steane, Active stabilization, quantum computation, and quantum state synthesis, Physical Review Letters 78, 2252 (1997).
  * Walshe  _et al._ [2020] B. W. Walshe, B. Q. Baragiola, R. N. Alexander, and N. C. Menicucci, Continuous-variable gate teleportation and bosonic-code error correction, Physical Review A 102, 062411 (2020).
  * Knill [2005] E. Knill, Scalable quantum computing in the presence of large detected-error rates, Physical Review A 71, 042322 (2005).
  * Chamberland and Jochym-O’Connor [2017] C. Chamberland and T. Jochym-O’Connor, Error suppression via complementary gauge choices in reed-muller codes, Quantum Science and Technology 2, 035008 (2017).
  * Tuckett [2020] D. K. Tuckett, Tailoring surface codes: Improvements in quantum error correction with biased noise, (2020).
  * Varnava  _et al._ [2006] M. Varnava, D. E. Browne, and T. Rudolph, Loss tolerance in one-way quantum computation via counterfactual error correction, Physical review letters 97, 120501 (2006).
  * Barrett and Stace [2010] S. D. Barrett and T. M. Stace, Fault tolerant quantum computation with very high threshold for loss errors, Physical review letters 105, 200502 (2010).
  * Tillich and Zémor [2013] J.-P. Tillich and G. Zémor, Quantum ldpc codes with positive rate and minimum distance proportional to the square root of the blocklength, IEEE Transactions on Information Theory 60, 1193 (2013).
  * Bonilla-Ataides  _et al._ [2020] J. P. Bonilla-Ataides, D. K. Tuckett, S. D. Bartlett, S. T. Flammia, and B. J. Brown, The xzzx surface code, arXiv preprint arXiv:2009.07851 (2020).
  * Hänggli  _et al._ [2020] L. Hänggli, M. Heinze, and R. König, Enhanced noise resilience of the surface–gottesman-kitaev-preskill code via designed bias, Physical Review A 102, 052408 (2020).