# Covariant density functional theory input for r-process simulations in
actinides and superheavy nuclei: the ground state and fission properties

A. Taninah  Department of Physics and Astronomy, Mississippi State University,
MS 39762  S. E. Agbemava  Department of Physics and Astronomy, Mississippi
State University, MS 39762  Ghana Atomic Energy Commission, National Nuclear
Research Institute, P.O.Box LG80, Legon, Ghana.  A. V. Afanasjev  Department
of Physics and Astronomy, Mississippi State University, MS 39762

###### Abstract

The systematic investigation of the ground state and fission properties of
even-even actinides and superheavy nuclei with $Z=90-120$ from the two-proton
up to two-neutron drip lines with proper assessment of systematic theoretical
uncertainties has been performed for the first time in the framework of
covariant density functional theory (CDFT). These results provide a necessary
theoretical input for the r-process modeling in heavy nuclei and, in
particular, for the study of fission cycling. Four state-of-the-art globally
tested covariant energy density functionals (CEDFs), namely, DD-PC1, DD-ME2,
NL3* and PC-PK1, representing the major classes of the CDFT models are
employed in the present study. Ground state deformations, binding energies,
two neutron separation energies, $\alpha$-decay $Q_{\alpha}$ values and half-
lives and the heights of fission barriers have been calculated for all these
nuclei. Theoretical uncertainties in these physical observables and their
evolution as a function of proton and neutron numbers have been quantified and
their major sources have been identified. Spherical shell closures at $Z=120$,
$N=184$ and $N=258$ and the structure of the single-particle (especially,
high-$j$) states in their vicinities as well as nuclear matter properties of
employed CEDFs are two major factors contributing into theoretical
uncertainties. However, different physical observables are affected in a
different way by these two factors. For example, theoretical uncertainties in
calculated ground state deformations are affected mostly by former factor,
while theoretical uncertainties in fission barriers depend on both of these
factors.

##  I Introduction

The majority of the nuclei found in nature are formed in the astrophysical
rapid neutron-capture process (r-process). Indeed, the r-process is
responsible for the synthesis of approximately half of the nuclei in nature
beyond Fe [1] and it is the only process which leads to the creation of nuclei
heavier than Bi [2]. It takes place at extremely high neutron densities
($N_{n}\geq 10^{20}$ cm-3) which are high enough to make neutron capture
faster than $\beta$-decay even for the nuclei with neutron excess between 15
to 30 units from the stability line. The production of neutron-rich nuclei
located in the vicinity of the neutron dripline is enabled under these
conditions via neutron capture and $(\gamma,n)$ photodisintegration during the
r-process. Once the neutron source ceases, the progenitor nuclei decay either
via $\beta^{-}$ decay or $\alpha$ emission or by fission processes (such as
neutron-induced, $\beta$-delayed and spontaneous fissions) towards stability
and form the stable isotopes of elements up to the heaviest species Th, U and
Pu. The typical timescale of the r-process is in the seconds range [3, 2, 4].

Over the years different possible astrophysical sites have been and still are
considered as possible candidates for the r-process. These include core-
collapse supernovae, magneto-rotational core-collapse supernovae, accretion
disk outflows from collapsars, neutron star (NS) mergers and neutron star -
black hole mergers etc [3, 2, 4, 5]. So far only the NS merger is
experimentally confirmed as a site of the r-process via the observation of
gravitational waves from the GW170817 neutron star merger [6] with
simultaneous observation of the AT 2017gfo macronova/kilonova afterglow [7].
In NS mergers, the r-process material originates in the NS crust, and the
composition of the crust and how it responds to stress caused by the merger
dictates the amount of the r-process material which is ejected. NS merger
produces approximately $10^{-2}$ M⊙ of ejected r-process matter in the dynamic
ejecta and similar amount in the accretion disk outflows [8, 9, 4]. Although
some uncertainties still exist, at present the NS merger is considered as the
major astrophysical site of the r-process providing the dominant source of
heavy nuclei [5, 2, 4].

The modeling of the r-process in such neutron-rich environments depends
sensitively on nuclear masses, $\alpha$\- and $\beta$-decay half-lives,
neutron capture and fission properties of the nuclei the majority of which
will never be measured in laboratory conditions [10, 5, 2, 4]. Nuclear masses
determine the flow path of the r-process, $\beta$-decay rates are responsible
for the speed with which the r-process moves matter to heavier nuclei,
$\alpha$-decays become important in heavy nuclei as a competing decay channel
and neutron captures drive the nuclei towards neutron-rich side of nuclear
landscape. Of special interest in the context of the present manuscript are
fission properties. Fission needs to be considered in the r-process
simulations if the neutron-to-seed ratio is large enough to produce fissioning
nuclei [11, 12, 13, 14]. If the initial neutron-to-seed ratio is large ($\geq
100$) the r-process can reach the region near and beyond neutron shell closure
at $N=184$, where fission plays a dominant role (the examples of the
distribution of abundances of actinides and superheavy elements as obtained in
a pair of the r-process simulations are shown in Fig. 2). This is exactly the
case for the NS mergers [12]. In this case, all fission channels (neutron
induced, beta-delayed, neutrino induced and spontaneous fissions) need to be
considered. Fission leads to the termination of hot r-process by means of
fission-cycling which returns matter to lighter nuclei [10]. It also
determines the strength of fission cycling, the ratio of the actinides to
light and medium mass r-process nuclei and thus the shape of final element
abundance pattern. In addition, it defines the possibility of the formation of
neutron-rich superheavy nuclei in the r-process [15].

Figure 1: The part of nuclear chart under study. Black, green, blue and orange
lines are used for the results obtained with the DD-PC1, DD-ME2, NL3* and PC-
PK1 CEDFs, respectively. Two-proton and two-neutron drip lines predicted by
four CEDFs are shown by dashed lines. Two samples of the distribution of the
abundances of heavy and superheavy elements in the r-process simulations are
shown by open squares (based on bottom panel of Fig. 8 of Ref. [15]) and red
circles (based on Fig. 3 of Ref. [16]). The former results correspond to hot
r-process conditions and are based on the ETFS/ETFSI combination of fission
barriers and mass predictions (see Ref. [15]). The latter results have been
obtained in Ref. [16] based on fission properties obtained with the BCPM
energy density functional. Note that these r-process calculations are
restricted to the $Z\leq 110$ nuclei. The r-process path is shown here
approximately by solid lines corresponding to two-neutron separation energy
$S_{2n}=4.0$ MeV222The r-process proceeds along the lines of constant neutron
separation energies towards heavy nuclei that for typical conditions during
the r-process corresponds to $S^{0}_{n}$ approximately located between 2 and 3
MeV [17, 1]. However, due to neutron pairing being stronger at even neutron
numbers $N$ the most abundant isotopes always have even $N$ values. For this
reason, we follow Refs. [17, 1] and characterize the r-process path [the path
in the $(Z,N)$ plane corresponding to an isotope with highest abundance in
each isotopic chain] as the path which satisfies the condition that two-
neutron separation energy $S_{2n}$ has the value $S_{2n}=2S^{0}_{n}$.. Dot-
dashed lines show the beta-stability lines for four functionals.

Thus, in the situation when experimental data are not known the outcome of the
r-process modeling sensitively depends on the quality of employed theoretical
frameworks and associated theoretical uncertainties and their propagation on
going to neutron-rich nuclei. By tradition, the output of different
theoretical frameworks is used for different physical observables (such as
masses, the rates and half-lives of different decay channels and reactions
etc.) in the the r-process modeling. Existing r-process calculations, which
include information on fission properties, are based on the fission barrier
heights obtained in non-relativistic models [14, 10, 18]. So far fission
barrier heights obtained in finite range droplet model (FRDM), Thomas-Fermi
(TF) model and extended Thomas-Fermi model with Strutinsky integral approach
(ETFSI-Q), Hartree-Fock-Bogoliubov (HFB) model with Skyrme HFB-14 energy
density functional (EDF) and BCPM EDF have been used in these calculations.
Moreover the sets of fission barriers relevant for the r-process simulations
have been generated in the FRDM in Ref. [18], in the ETFSI-Q approach with the
SkSC4 functional in Ref. [19], in the HFB models with Skyrme HFB-14, SV-min,
SLy6, SkI3, SV-bas EDFs [20, 21, 22], Gogny D1M* [23] and BCPM [18]
functionals. Note that all these calculations assume axial symmetry of the
nuclei.

Covariant density functional theory (CDFT) [24] is an approach alternative to
above mentioned non-relativistic methods and so far it has not been applied
for a systematic study of fission properties of the nuclei relevant for the
r-process modeling. However, this theory has a number of advantages over non-
relativistic methods which are discussed below. Covariant energy density
functionals (CEDF) exploit basic properties of QCD at low energies, in
particular symmetries and the separation of scales [25]. They provide a
consistent treatment of the spin degrees of freedom and spin-orbit splittings
([26, 27]); the latter has an essential influence on the underlying shell
structure. In addition, these functionals include nuclear magnetism [28], i.e.
a consistent description of currents and time-odd mean fields important for
odd-mass nuclei [29], the excitations with unsaturated spins, magnetic moments
[30] and nuclear rotations [31, 32]. Because of Lorentz invariance no new
adjustable parameters are required for the time-odd parts of the mean fields
[29]. This is contrary to the case of non-relativistic Skyrme DFTs in which
several prescriptions for fixing time-odd mean fields exist [33, 34]333
Unfortunately, the role of time-odd mean fields in Gogny DFTs has not been
studied so far and it is unknown whether they are uniquely defined.. This fact
could be extremely important in the applications to fission processes
including dynamical correlations since time-odd mean fields have a significant
impact on collective masses [35, 36].

The goal of this study is to close this gap in our knowledge and to perform
first systematic investigation within the CDFT framework of the ground state
and fission properties of the nuclei with proton numbers $Z=90-120$ located
between two-proton and two-neutron drip lines (see Fig. 2). This study will
not only provide an input for the r-process modeling but also evaluate the
extension of nuclear landscape up to two-neutron drip line as well as estimate
relevant theoretical uncertainties and their sources in the description of
physical observables of interest. In addition, it will allow for the first
time to compare the predictions for fission barriers in the nuclei relevant
for r-process modeling obtained in relativistic and non-relativistic
frameworks.

Considering the region of the nuclear chart in which the r-process is expected
to take place and the fact that there are no experimental data to benchmark
theoretical results, it is important to estimate theoretical uncertainties in
the predictions of physical observables of interest [37, 38, 39]. Theoretical
uncertainties emerge from the underlying theoretical approximations. In the
DFT framework, there are two major sources of these approximations, namely,
the range of interaction and the form of the density dependence of the
effective interaction [40, 41]. In the non-relativistic case one has zero
range Skyrme and finite range Gogny forces and different density dependencies
[40]. A similar situation exists also in the relativistic case: point coupling
and meson exchange models have an interaction of zero and of finite range,
respectively [24, 42, 43, 44]. The density dependence is introduced either
through an explicit dependence of the coupling constants [45, 42, 44] or via
non-linear meson couplings [41, 43]. This ambiguity in the definition of the
range of the interaction and its density dependence leads to several major
classes of the covariant energy density functionals (CEDF) which were
discussed in Ref. [39].

Since statistical uncertainties in the physical observables are smaller than
systematic ones (see Ref. [46]), we focus here on the latter ones. They are
related to the choice of EDF. We follow our previous publications on this
topic [39, 47, 48, 49] and define systematic theoretical uncertainty for a
given physical observable (which we call in the following “spreads”) via the
spread of theoretical predictions as [39]

$\Delta O(Z,N)=|O_{max}(Z,N)-O_{min}(Z,N)|,$ (1)

where $O_{max}(Z,N)$ and $O_{min}(Z,N)$ are the largest and smallest values of
the physical observable $O(Z,N)$ obtained within the set of CEDFs under
investigation for the $(Z,N)$ nucleus. Note that these spreads are only a
crude approximation to the systematic theoretical errors discussed in Ref.
[38] since they are obtained with a very small number of functionals which do
not form an independent statistical ensemble. Note also that these systematic
errors are not well defined in unknown regions of nuclear chart or deformation
since systematic biases of theoretical models could not be established in
these regions in the absence of experimental data and/or an exact theory.

In order to consider several possible scenarios in the evolution of physical
observables as a function of proton and neutron numbers and to evaluate
systematic theoretical uncertainties, the CEDFs NL3* [43], DD-ME2 [42], DD-PC1
[44] and PC-PK1 [50] are used here444The compilation of Ref. [51] published in
2014 indicates the existence of 263 CEDFs ranging from simplest ansatz non-
linear meson exchange functionals such as NL1 [52] and NL3* [43] to more
microscopically motivated CEDFs such as G1, G2 [53] and DD-ME$\delta$ [54]. In
addition, a number of new functionals were fitted in the time period between
2014 and 2020 (see, for example, Refs. [55, 56, 57, 46, 58]) so at present the
total number of available CEDFs is likely to be in the vicinity of 300.
Because of extremely time-consuming nature of numerical calculations in this
project, we use only the indicated last generation functionals. They
outperform previous generation functionals in terms of the accuracy of global
description of ground state observables such as binding energies and charge
radii [39, 49, 59, 60], properly describe the regions of octupole deformation
[48, 61] and are able to reproduce experimentally known fission barriers in
actinides [62, 63, 64]. for all $Z=90-120$ even-even nuclei located between
two-proton and two-neutron drip lines555Present study partially builds on
previous results obtained by us. These are ground state properties of the
$Z=90-104$ nuclei located between two-proton and two-neutron drip lines
obtained in reflection symmetric RHB calculations with the NL3*, DD-ME2 and
DD-PC1 CEDFs in Ref. [39], ground state properties of octupole deformed nuclei
with $N<210$ obtained in reflection asymmetric RHB calculations with NL3*, DD-
ME2, DD-PC1 and PC-PK1 functionals in Ref. [48, 61] and ground state
properties and the heights of inner fission barriers of superheavy nuclei with
$Z=100-120,N\leq 196$ obtained in reflection symmetric RHB calculations with
NL3*, DD-ME2, DD-PC1 and PC-PK1 functionals in Refs. [59, 47]. Additional
information on the extension of nuclear landscape to $Z>120$ obtained with DD-
PC1 CEDF can be found in Refs. [65, 66, 67].. These are state-of-the-art
functionals representing the major classes of CDFTs (for more details see the
discussion in Sect. II of Ref. [39] and the introduction to Ref. [59]). Their
performance and related theoretical uncertainties have recently been analyzed
globally in Refs. [39, 68, 49, 48] and in particular in superheavy nuclei in
Refs. [59, 47]. They are characterized by an improved accuracy of the
description of experimental data as compared with the previous generation of
CEDFs. The fact that the NL3*, DD-PC1 and PC-PK1 functionals reproduce
empirical data on fission barrier heights in actinides [62, 63, 69, 64] is
especially important in the context of the present study.

The manuscript is organized as follows. Theoretical framework and the details
related to the calculations of the ground states and fission barriers are
discussed in Sec. II. Sec. III is devoted to the analysis of the results of
the calculations for ground state properties. Theoretical results for
$\alpha$-decay properties and related theoretical uncertainties are presented
in Sec. IV. The heights of primary fission barriers, their distribution in the
$(Z,N)$ plane, related theoretical uncertainties and the comparison with non-
relativistic results are considered in Sec. V. Finally, Sec. VI summarizes the
results of our work.

Figure 2: Schematic illustration of different types of deformation energy
curves and the selection of respective ground states (see text for details).
Local minima are labelled by the letters A and B and the saddle points of
respective fission barriers by $\rm S_{A}$ and $\rm S_{B}$. Solid circles
indicate the minima selected as the ground states. The heights of inner and
outer fission barriers with respect of corresponding minima (shown by dashed
lines) are indicated by $\rm B_{in}$ and $\rm B_{out}$. $\rm B_{out-iso}$ is
the height of outer fission barrier with respect of fission isomer.

##  II Theoretical framework

Figure 3: Deformation energy curves obtained in axial RHB calculations with
DDPC1 functional for the Ds isotopes. The isotopes are indicated by respective
neutron numbers. Thick and thin solid lines are used for the RS-RHB and RA-RHB
results, respectively. The results of the RA-RHB calculations are shown only
in the deformation range in which they are lower in energy than the RS-RHB
ones. Solid circles indicate the ground states and the asterisks denote the
saddles of outer fission barriers which are not affected by octupole
deformation. Blue, orange, black, red and green lines are used to indicate
neutron numbers the last digits of which are 0, 2, 4, 6, and 8, respectively.
Note that the energy on the vertical axis spans different ranges in different
panels.

Figure 4: The same as Fig. 3 but for the Th isotopes. Note that the
deformation range has been extended on horizontal axis as compared with Fig.
3. In order to save computational time the RA-RHB calculations have been
carried out only up to $\beta_{2}=2.0$.

In the present manuscript, the RHB framework with finite range pairing and its
separable limit are used for a systematic study of ground state properties of
all even-even actinides ($Z=90-102$) and superheavy ($Z=104-120$) nuclei from
the proton to neutron drip line. It has the proper coupling to the continuum
at the neutron drip line and, therefore, it allows a correct description of
weakly bound nuclei close to the neutron drip line.

Figure 5: Proton quadrupole deformations $\beta_{2}$ obtained in the RS-RHB and RA-RHB calculations with indicated CEDFs (see Refs. [48, 61] for the details of the regions of octupole deformation). Note that last bound $Z=120$ nucleus appears at $N=324$ in the calculations with the CEDF PC-PK1 (see Table 3 below); for simplicity it is not shown in panel (d).  Table 1: The nuclei in which extremely superdeformed minimum is the lowest in energy in the calculations with CEDF DD-PC1. The columns 1 and 2 show the proton and neutron numbers of the nuclei. Third column displays the energy $E_{N-S}$ [in MeV] by which the ESD minimum is lower than the normal-deformed minimum. The deformations of the ESD minimum $\beta_{2}^{min}$ and the saddle of second fission barrier $\beta_{2}^{saddle}$, $\beta_{3}^{saddle}$ are presented in the columns 4 and 5, respectively. The energies [in MeV] of the second fission barrier with respect of the ESD minimum, obtained in the RS-RHB and RA-RHB calculations, are shown in the columns 6 and 7. Note that the values presented in the columns 5-7 are obtained in the calculations with $N_{F}=26$; this is done in order to have a comparable numerical accuracy with the one obtained at normal deformed minimum.  Z | N | $E_{N-S}$ | $\beta_{2}^{min}$ | $\beta_{2}^{saddle},\beta_{3}^{saddle}$ |  $E_{B}^{II}$[MeV] |  $E_{B}^{II}$ [MeV]
---|---|---|---|---|---|---
|  | [MeV] |  |  | (RS-RHB) | (RA-RHB)
1 | 2 | 3 | 4 | 5 | 6 | 7
98 | 142 | 0.603 | 0.88 | 1.21, 0.33 | 7.207 | 2.625
98 | 144 | 0.306 | 0.90 | 1.34, 0.46 | 7.373 | 3.503
98 | 228 | 2.341 | 1.01 | 1.29, 0.34 | 6.622 | 2.906
98 | 230 | 2.083 | 1.01 | 1.30, 0.35 | 5.072 | 3.203
100 | 146 | 0.876 | 0.97 | 1.32, 0.37 | 5.890 | 2.985
100 | 230 | 2.431 | 1.01 | 1.30, 0.35 | 5.750 | 3.270
100 | 232 | 2.336 | 1.01 | 1.30, 0.37 | 4.060 | 2.739
102 | 146 | 2.038 | 0.99 | 1.32, 0.34 | 4.813 | 2.749
102 | 148 | 0.629 | 0.98 | 1.28, 0.31 | 4.093 | 2.250
102 | 232 | 2.591 | 1.02 | 1.29, 0.36 | 3.476 | 2.840
102 | 234 | 3.567 | 1.03 | 1.31, 0.31 | 2.674 | 2.304
104 | 146 | 3.435 | 0.99 | 1.33, 0.28 | 6.042 | 2.410
104 | 148 | 1.921 | 1.00 | 1.31, 0.27 | 5.579 | 2.204
104 | 150 | 0.924 | 1.00 | 1.29, 0.27 | 5.370 | 2.505
106 | 148 | 3.638 | 1.08 | 1.32, 0.26 | 4.112 | 2.340
106 | 150 | 2.374 | 1.10 | 1.30, 0.23 | 3.753 | 2.345

Figure 6: (a) Deformation energy curves obtained in axial RS-RHB calculations
with the CEDF DD-PC1 for the 240Cf nucleus. (b) Potential energy surface in
the $(\beta_{2},\beta_{3}$) plane obtained in the RA-RHB calculations.
Extremely superdeformed minimum is indicated by open white circle.

Figure 7: Proton quadrupole deformation spreads $\Delta\beta_{2}(Z,N)$ as a
function of proton and neutron number. $\Delta\beta_{2}(Z,N)=|\beta_{2}^{\rm
max}(Z,N)-\beta_{2}^{\rm min}(Z,N)|$, where $\beta_{2}^{\rm max}(Z,N)$ and
$\beta_{2}^{\rm min}(Z,N)$ are the largest and smallest proton quadrupole
deformations obtained with four employed CEDF for the $(Z,N)$ nucleus.

The RHB equations for the fermions are given by [70]

$\displaystyle\begin{pmatrix}\hat{h}_{D}-\lambda&\hat{\Delta}\\\
-\hat{\Delta}^{*}&-\hat{h}_{D}^{\,*}+\lambda\end{pmatrix}\begin{pmatrix}U({\bm{r}})\\\
V({\bm{r}})\end{pmatrix}_{k}=E_{k}\begin{pmatrix}U({\bm{r}})\\\
V({\bm{r}})\end{pmatrix}_{k},$ (2)

Here, $\hat{h}_{D}$ is the Dirac Hamiltonian for the nucleons with mass $m$;
$\lambda$ is the chemical potential defined by the constraints on the average
particle number for protons and neutrons; $U_{k}({\bm{r}})$ and
$V_{k}({\bm{r}})$ are quasiparticle Dirac spinors [71, 72, 70] and $E_{k}$
denotes the quasiparticle energies. The Dirac Hamiltonian

$\hat{h}_{D}=\bm{\alpha}(\bm{p}-\bm{V})+V_{0}+\beta(m+S).$ (3)

contains an attractive scalar potential

$\displaystyle S(\bm{r})=g_{\sigma}\sigma(\bm{r}),$ (4)

a repulsive vector potential

$\displaystyle
V_{0}(\bm{r})~{}=~{}g_{\omega}\omega_{0}(\bm{r})+g_{\rho}\tau_{3}\rho_{0}(\bm{r})+eA_{0}(\bm{r}),$
(5)

and a magnetic potential

$\displaystyle\bm{V}(\bm{r})~{}=~{}g_{\omega}\bm{\omega}(\bm{r})+g_{\rho}\tau_{3}\bm{\rho}(\bm{r})+e\bm{A}(\bm{r}).$
(6)

The last term breaks time-reversal symmetry and induces currents. For example,
time-reversal symmetry is broken when the time-reversed orbitals are not
occupied pairwise; this takes place in odd-mass nuclei [29]. However, nuclear
magnetism [28], i.e. currents and time-odd mean fields, plays no role in the
studies of ground states and fission barriers in even-even nuclei. Thus,
magnetic potential is neglected in the present RHB calculations.

In order to avoid the uncertainties connected with the definition of the size
of the pairing window [73], we use the separable form of the finite range
Gogny pairing interaction introduced by Tian et al [74]. Its matrix elements
in $r$-space have the form

$\displaystyle
V({\bm{r}}_{1},{\bm{r}}_{2},{\bm{r}}_{1}^{\prime},{\bm{r}}_{2}^{\prime})$
$\displaystyle=$ $\displaystyle=-G\delta({\bm{R}}-$
$\displaystyle\bm{R^{\prime}}$
$\displaystyle)P(r)P(r^{\prime})\frac{1}{2}(1-P^{\sigma})$ (7)

with ${\bm{R}}=({\bm{r}}_{1}+{\bm{r}}_{2})/2$ and
${\bm{r}}={\bm{r}}_{1}-{\bm{r}}_{2}$ being the center of mass and relative
coordinates. The form factor $P(r)$ is of Gaussian shape

$\displaystyle P(r)=\frac{1}{(4\pi a^{2})^{3/2}}e^{-r^{2}/4a^{2}}$ (8)

The two parameters $G=728$ MeV fm3 and $a=0.644$ fm of this interaction are
the same for protons and neutrons and have been derived in Ref. [74] by a
mapping of the 1S0 pairing gap of infinite nuclear matter to that of the Gogny
force D1S [75]. This pairing provides a reasonable description of pairing
properties in the actinides (see Refs. [76, 39, 77]) and has been used in our
previous studies of different phenomena in actinides, super- and hyperheavy
nuclei in Refs. [39, 59, 48, 61, 65, 66].

The truncation of the basis is performed in such a way that all states
belonging to the major shells up to $N_{F}=20$ fermionic shells for the Dirac
spinors and up to $N_{B}=20$ bosonic shells for the meson fields are taken
into account. Note that the latter applies only to the NL3* and DD-ME2
functionals which contain meson exchange. As follows from investigation of
Refs. [39, 78] this truncation of basis provides sufficient numerical
accuracy.

The calculations are performed in the following way:

  * •

Reflection-symmetric constrained axial RHB calculations (further RS-RHB) are
performed for each nucleus and the potential energy curve is defined in a
large deformation range from $\beta_{2}=-1.0$ up to $\beta_{2}=1.6$ in step of
$\Delta\beta_{2}=0.05$ by means of the constraint on the quadrupole moment
$q_{20}$. The calculations are performed by successive diagonalizations using
the method of quadratic constraints [79]. The parallel version of computer
code allows simultaneous calculations for a significant number of nuclei and
deformation points in each nucleus. For each nucleus, we minimize

$E_{quad}=E_{RHB}+C_{20}(\langle\hat{Q}_{20}\rangle-q_{20})^{2}$ (9)

where $E_{RHB}$ is the total energy and $\langle\hat{Q}_{20}\rangle$ denotes
the expectation value of the mass quadrupole operator,

$\hat{Q}_{20}=2z^{2}-x^{2}-y^{2}$ (10)

$q_{20}$ is the constrained value of the multipole moment, and $C_{20}$ the
corresponding stiffness constant [79]. In order to provide the convergence to
the exact value of the desired multipole moment we use the method suggested in
Ref. [80]. Here the quantity $q_{20}$ is replaced by the parameter
$q_{20}^{eff}$, which is automatically modified during the iteration in such a
way that we obtain $\langle\hat{Q}_{20}\rangle=q_{20}$ for the converged
solution. This method works well in our constrained calculations.

  * •

In addition, reflection-asymmetric (octupole deformed) constrained axial RHB
calculations (further RA-RHB) are performed in discussed below cases using
parallel version of the code developed in Ref. [48]. In these calculations,
the constraint

$E_{quad}+C_{30}(\langle\hat{Q}_{30}\rangle-q_{30})^{2}$ (11)

is employed in addition to the constraint on quadrupole moment (see Eq. (9)).
Here $\langle\hat{Q}_{30}\rangle$ denotes the expectation value of the mass
octupole operator

$\hat{Q}_{30}=z(2z^{2}-3x^{2}-3y^{2}).$ (12)

Note that we also fix the (average) center of mass of the nucleus at the
origin with the constraint

$\langle\hat{Q}_{10}\rangle=0$ (13)

on the center-of-mass operator ${\hat{Q}}_{10}$ in order to avoid a spurious
motion of the center of mass. In the present paper, reflection asymmetric RHB
calculations have been performed for the ground states of the nuclei not
covered in previous systematic studies of octupole deformation in CDFT (see
Refs. [48, 61]). We have not found any additional (as compared with those
given in Refs. [48, 61]) nuclei which possess octupole deformation in the
ground state. So full information on the octupole deformation of the ground
states can be found in these references. The information (which follows from
Refs. [48, 61]) about the gain in binding energy due to octupole deformation
at the ground state and its impact on ground state quadrupole deformation and
fission barrier heights is fully taken into account in the present paper. In
addition, RA-RHB calculations have been performed in some nuclei in order to
define the heights of outer fission barriers (see the discussion below for
more details).

Figure 8: The binding energy spreads $\Delta E(Z,N)$ as a function of proton
and neutron number. $\Delta E(Z,N)=|E_{\rm max}(Z,N)-E_{\rm min}(Z,N)|$, where
$E_{\rm max}(Z,N)$ and $E_{\rm min}(Z,N)$ are the largest and smallest binding
energies obtained with employed set of CEDFs for the $(Z,N)$ nucleus. Upper
panel shows the binding energy spreads for four employed functionals, while
the bottom one the spreads for a set of functionals in which PC-PK1 is
excluded.

The charge quadrupole and octupole moments are defined as

$\displaystyle Q_{20}$ $\displaystyle=$ $\displaystyle\int
d^{3}r\rho({\bm{r}})\,(2z^{2}-r^{2}_{\perp}),$ (14) $\displaystyle Q_{30}$
$\displaystyle=$ $\displaystyle\int
d^{3}r\rho({\bm{r}})\,z(2z^{2}-3r^{2}_{\perp})$ (15)

with $r^{2}_{\perp}=x^{2}+y^{2}$. In principle these values can be directly
compared with experimental data. However, it is more convenient to transform
these quantities into dimensionless deformation parameters $\beta_{2}$ and
$\beta_{3}$ using the relations

$\displaystyle Q_{20}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{16\pi}{5}}\frac{3}{4\pi}ZR_{0}^{2}\beta_{2},$ (16)
$\displaystyle Q_{30}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{16\pi}{7}}\frac{3}{4\pi}ZR_{0}^{3}\beta_{3}$ (17)

where $R_{0}=1.2A^{1/3}$. These deformation parameters are more frequently
used in experimental works than quadrupole and octupole moments. In addition,
the potential energy surfaces (PES) are plotted in this manuscript in the
($\beta_{2},\beta_{3}$) deformation plane.

Because of different patterns of deformation energy curves (see Figs. 3 and
4), a special care is used when assigning a specific minimum to the ground
state. A basic rule in this process is the assumption that local minimum
surrounded by the barrier, the height of which is less than 2 MeV, is
considered as extremely unstable666This low fission barrier of 2 MeV or less
would translate into a high penetration probability for spontaneous fission so
that such minima are expected to be extremely unstable. In addition, the
inclusion of octupole deformation (or triaxial deformation in some nuclei
[78]) in the case of superdeformed minima surrounded by such low fission
barriers could either completely eliminate or substantially reduce them (see
Refs. [81, 78, 66]).. The procedure of the selection of the ground state is
discussed below. The situation shown in Fig. 2a is the simplest one: single-
humped (inner) fission barrier acts against the fission into two fragments and
when $\rm B_{in}>2$ MeV the assignment of the normal-deformed prolate minimum
to the ground state is straightforward. It changes if the height of this
fission barrier decreases and becomes smaller than 2 MeV (see Fig. 2b). Then
highly-deformed oblate minimum B becomes a ground state; it has larger and
broader fission barrier as compared with minimum A. Thus, it is expected that
this ground state will live significantly longer than the state associated
with minimum A.

Figure 9: Two-neutron separation energies $S_{2n}(Z,N)$ obtained in the RHB
calculations with indicated CEDFs.

More complex situation involving two humped fission barrier is shown in Fig.
2(c). If $B_{out-iso}<2$ MeV, then the minimum B corresponding to fission
isomer is considered extremely unstable and the minimum A is associated with
the ground state. Note that in some cases the minimum B can be lower in energy
than minimum A. If that is the case and if $B_{out-iso}>2$ MeV then the
minimum B is considered as the ground state. Note that the heights of outer
fission barriers are frequently lowered when octupole deformation is included
in the calculations [81, 82, 83, 69, 78, 63]. Thus, if $B_{out-iso}>2$ MeV in
RS-RHB calculations, we perform RA-RHB calculations in the region of the
$(\beta_{2},\beta_{3})$ plane covering the minimum B and the saddle of outer
fission barrier on the grid with the steps of
$\Delta\beta_{2}=\Delta\beta_{3}=0.05$. This allows to establish whether
minimum B could be considered as relatively stable or unstable. Similar
calculations are also performed in the cases when $\rm B_{in}<B_{out}$ in the
RS-RHB calculations. This is because we consider only the height of the
primary (highest) fission barrier (PFB) in Sec. V in the case of doubly-humped
structure of the barrier to minimize the computational cost. Note that the
calculations leading to the definition of the fission path and the saddle
point in the RA-RHB code are by roughly two orders of magnitude more time
consuming than those in the RS-RHB code.

The procedure outlined above takes into account potential stability of the
nuclei in respective energy minima with respect of fission and it is
especially important in superheavy nuclei some local minima of which are
characterized by small fission barriers (see Fig. 3). Note that after defining
the minimum corresponding to the ground state, the RS-RHB and RA-RHB
calculations without constraint are performed in it and precise binding energy
and equilibrium of the ground state is determined. In addition, the height(s)
of the fission barrier(s) is(are) defined.

In the calculations with the PC-PK1 and NL3* functionals there are two small
islands of the nuclei located in the $Z\approx 114-118,N\approx 238-240$ and
$Z\approx 106-110,N\approx 190-194$ regions in which calculated deformation
energy curves reveal several local minima (somewhat similar to the deformation
energy curves shown at the bottom of Fig. 3b). However, all these minima are
surrounded by very low fission barriers with the heights smaller than 2.0 MeV.
Moreover, many of these minima have even lower heights of respective fission
barriers (on the level of 1.0 MeV or smaller). These nuclei are expected to be
unstable against fission and in principle it does not matter which of the
calculated fission barriers is used. For these nuclei, we select the ground
state guided by the flow of the $\beta$-decays in the r-process: the selected
local minimum (and thus the corresponding ground state deformation and fission
barrier height) in the $(Z,N)$ nucleus has the deformation closest to the one
of the well established ground state in the $(Z-2,N+2)$ nucleus.

Figure 10: The $S_{2n}$ spreads $\Delta S_{2n}(Z,N)$ as a function of proton
and neutron number. $\Delta S_{2n}(Z,N)=|S_{2n}^{\rm max}(Z,N)-S_{2n}^{\rm
min}(Z,N)|$, where $S_{2n}^{\rm max}(Z,N)$ and $S_{2n}^{\rm min}(Z,N)$ are the
largest and smallest $S_{2n}$ values obtained with four employed CEDFs for the
$(Z,N)$ nucleus.

So far all existing global calculations of the fission barriers for the
r-process simulations have been performed in non-relativistic models. These
include the calculations within the FRDM [84, 18], the TF and ETFSI-Q
approaches [19, 15], the HFB approaches with different Skyrme functionals [21,
22, 85], Gogny D1M* [23] and BCPM [18] EDFs. Because of their global
character, all these calculations are restricted to axial symmetry. We also
assume axial symmetry in our calculations because triaxial RMF+BCS and RHB
calculations (see Refs. [62, 78, 47, 66, 86]) are too time consuming to be
performed on a global scale. Note also that dynamical correlations are not
taken into account in our calculations of fission barriers because of the
reasons discussed in Appendix A.

##  III Ground state properties

The distributions of calculated proton deformations $\beta_{2}$ in the $(Z,N)$
plane obtained with four employed CEDFs are shown in Fig. 5. The width of the
gray region (the gray color corresponds to spherical and near-spherical
shapes) along a specific magic number corresponding to a shell closure
indicates the impact of this shell closure on the structure of neighboring
nuclei. Note that proton and neutron shell gaps act simultaneously in the
vicinity of doubly magic spherical nuclei. Thus, the effect of a single gap is
more quantifiable away from these nuclei. One can see that neutron $N=184$ and
$N=258$777 Note that appreciable $N=258$ spherical shell gap appears also in
the calculations of some superheavy nuclei with other CEDFs such as G1, G2
[87], NL3 [88, 87, 89], NL-Z2 [88], NL1, NLSH, TM1, TW99, DD-ME1, PK1, and
PK1R [89]. However, these calculations are restricted to spherical shapes and
thus it is not clear how large is the impact of this gap on ground state
deformations in the region near the $N=258$ line.  spherical shell gaps have a
pronounced impact on calculated deformations while the impact of the proton
$Z=120$ spherical shell closure is limited to the $N\sim 170-184$ nuclei (see
Refs. [90, 59] for examples of their size dependence on the functional). In
addition, as illustrated in Ref. [86] on the example of the PC-PK1 functional,
the inclusion of the correlations beyond mean field washes out the impact of
the $Z=120$ shell closure and leads to oblate deformed ground states in the
majority of the $Z=120$ nuclei with $N=172-186$. The predictions of the DD-PC1
functional differs substantially from other CEDFs: the impact of above
mentioned shell closures are substantially reduced in it and as a consequence
the regions with $Z\sim 120,N\sim 184$ and $Z\sim 120,N\sim 258$ are dominated
by oblate ground states contrary to spherical states in other functionals.
Note that this functional provides the best global description of experimental
binding energies (see Ref. [39]). This, however, does not guarantee that it
will outperform other functionals in the description of physical observables
of interest in the region of superheavy nuclei (see Table I in Ref. [59]).

The calculations reveal a number of nuclei scattered across the part of
nuclear chart under study which have extremely superdeformed (ESD) minimum
with $\beta_{2}\sim 1.0$ located at lower energy than normal-deformed minimum
(see Fig. 6a and Table 1). In these nuclei the ESD minimum is surrounded by
outer fission barrier the height of which exceeds 2.0 MeV (being typically in
the range of $2.0-3.0$ MeV) in the RA-RHB calculations (see Fig. 6b and Table
1). Although the fission barrier is low, some of these ESD minima could be
potentially stabilized against fission for physically sufficient time because
of broad fission barrier in the $(\beta_{2},\beta_{3})$ plane. If that would
be a case, they would represent the ground states. However, they are not
included into Fig. 5 because of the following reasons. First, there are
significant theoretical uncertainties in the predictions of fission barriers
(see present paper and Refs. [78, 47]) as well as in relative energies of the
minima with different deformations [73]. Second, only few nuclei in the
$Z\approx 100,N\approx 230$ region could be potentially extremely supedeformed
in the ground states (see Table 1). However, the flow of matter in the
r-process between two nuclei with drastically different deformations of the
ground states will be most likely significantly suppressed because of
considerable differences in the wave functions of these ground states. Thus,
it will proceed mostly along the dominant deformation of the ground states in
the region, namely, normal deformation, even if such states are excited in
energy in a few nuclei. Third, the majority of the nuclei in Table 1 are
neutron poor $Z\approx 102$ nuclei which do not play a role in the r-process.
There are experimental data on the 240,242Cf, 246Fm and 254Rf nuclei but only
for their ground states [91]. At present, these data do not allow to define
the deformations of the ground states. However, since it has been obtained in
the reactions (such as $\alpha$-decay, $\beta$-decay, electron capture and the
reactions on spherical Pb isotopes) which do not favor significant shape
changes, these ground states are most likely normal-deformed. More detailed
and focused experimental studies are needed in order to see whether ESD states
exist in such nuclei.

Table 2: Selected properties of symmetric nuclear matter at saturation: the incompressibility $K_{0}$, the symmetry energy $J$ and its slope $L_{0}$. Top four lines show the values for indicated CEDFs, while bottom two lines show two sets (SET2a and SET2b) of the constraints on the experimental/empirical ranges for the quantities of interest defined in Ref. [51]. The CEDF values which are located outside the limits of the SET2b constraint set are shown in bold.  CEDF |  $K_{0}$ [MeV] |  $J$ [MeV] |  $L_{0}$ [MeV]
---|---|---|---
1 | 2 | 3 | 4
NL3* [43] | 258 | 38.68 | 122.6
DD-ME2 [42] | 251 | 32.40 | 49.4
DD-PC1 [44] | 230 | 33.00 | 68.4
PC-PK1 [50] | 238 | 35.6 | 113
SET2a | 190-270 | 25-35 | 25-115
SET2b | 190-270 | 30-35 | 30-80

The spreads of theoretical predictions in quadrupole deformations $\beta_{2}$
obtained with four employed functionals are summarized in Fig. 7(a). The
largest spread of $\Delta\beta_{2}\approx 0.7$ is visible along the line of
$N/Z\approx 1.81$ which starts at $Z=104$. This corresponds to the boundary of
the transition from oblate to prolate shapes the exact position of which in
the $(Z,N)$ plane is functional dependent (see Fig. 5). It is defined by the
underlying single-particle structure at prolate and oblate shapes as well as
to a degree by the heights of outer fission barriers (see the rules for the
definition of the ground states described in Sec. II). Second region of the
largest spreads in $\Delta\beta_{2}$ is located along the $N\approx 184$ line
starting from $Z\approx 100$ and extending up to $Z=120$. Third region is
located along the $Z=120$ line from proton-drip line up to $N\approx 188$.
These two regions of large spreads in calculated quadrupole deformation emerge
from the differences in the predictions of ground state deformations (see Fig.
5) which in turn can be traced back to the sizes of the $Z=120$ and $N=184$
spherical shell closures and the densities of the single-particle states in
their vicinities (see Ref. [59]). The last region of large theoretical
uncertainties is located between $N\approx 236$ and $N=258$. In the region
around $N\approx 236$ these theoretical uncertainties are mostly due to the
uncertainties in the predictions of the boundary of the transition from
prolate to oblate shapes. For higher $N$ values, large $\Delta\beta_{2}$
values emerge from the transition from prolate or oblate shapes to spherical
ones and to a large degree are defined by the uncertainties in the prediction
of the size of the $N=258$ spherical shell closure (see Fig. 6d in Ref. [90])
and single-particle densities in its vicinity. With few exceptions theoretical
uncertainties in the predictions of ground state deformations in the part of
nuclear chart outside of above discussed regions are very small (see Fig.
7(a)).

It is important to understand to what extent the predictions of the ground
state deformations and related theoretical uncertainties in these predictions
are dependent on nuclear matter properties of employed CEDFs. All employed
CEDFs have the density $\rho_{0}$ and the energy per particle $E/A$ at the
saturation of symmetric nuclear matter (SNM) which are very close to each
other and to empirical estimates (see Table III in Ref. [49]). Thus, the
impact of only selected SNM properties listed in Table 2 on the ground state
deformations are discussed below. Let start from the consideration of the
predictions by the pair of the functionals DD-PC1 and DD-ME2. Their SNM
properties such as incompressibility $K_{0}$, the symmetry energy $J$ and its
slope $L_{0}$ are close to each other and are located within the SET2b
constraints on experimental/empirical ranges for physical observables of
interest (see Table 2). Despite that this pair of the functionals gives the
largest contribution into the spreads $\Delta\beta_{2}$ (compare panels (b)
and (a) of Fig. 7). On the contrary, the pairs of the functionals PC-PK1 and
NL3* (which have $J$ and $L_{0}$ values located outside the SET2b constraint
range [see Table 2]) as well as NL3* and DD-ME2 (which have drastically
different values of the $J$ and $L_{0}$ parameters [see Table 2]) have (with
very few exceptions) very similar predictions for the ground state
deformations across the part of nuclear chart under study. These exceptions
are related to some differences in the predictions of the boundaries between
oblate and prolate shapes as well as between prolate and spherical shapes.

These results for ground state deformations together with the analysis of the
results for binding energies and charge radii of the $Z\leq 104$ nuclei
presented in Ref. [49] strongly indicate that

  * •

the major source of the uncertainties in the predictions of ground state
deformations is related to local differences in underlying single-particle
structure and, in particular, to the size of spherical $Z=120$ and $N=184$ and
258 shell closures and the densities of the single-particle states in their
vicinities,

  * •

strict enforcement of the limits on the nuclear matter properties defined in
Ref. [51] will not necessary lead to the functionals with good description of
ground and excited state properties and will not reduce theoretical
uncertainties in the description of physical observables of interest in
high-$Z$ and/or neutron-rich nuclei.

Table 3: Two-proton and two-neutron drip lines predicted by the NL3* and PC-PK1999The analysis of Ref. [92] performed within the relativistic continuum Hartree-Bogoliubov theory and PC-PK1 functional leads to somewhat different predictions for the position of the two-neutron drip line for some isotopic chains as compared with our results. This is a consequence of the neglect of deformation effects in Ref. [92]. functionals (see Fig. 2 for graphical representation of drip lines). Neutron numbers (columns 2-5) corresponding to these drip lines are given for each even proton number $Z$ (column 1). An asterisk at a neutron number at the two-neutron drip line indicates isotope chains with additional two-neutron binding at higher $N$ values (peninsulas). Proton | Two-proton drip line | Two-neutron drip line
---|---|---
number $Z$ | NL3* | PCPK1 | NL3* | PCPK1
1 | 2 | 3 | 4 | 5
90 | 112 | 114 | 218 | 256
92 | 118 | 122 | 224 | 258
94 | 122 | 126 | 232 | 258
96 | 126 | 128 | 252 | 258
98 | 130 | 130 | 256 | 258
100 | 132 | 132 | 258 | 258
102 | 134 | 136 | 258 | 258
104 | 138 | 140 | 258 | 258
106 | 142 | 144 | 258 | 258*
108 | 146 | 148 | 258 | 288
110 | 150 | 154 | 258 | 292
112 | 154 | 158 | 258 | 298
114 | 158 | 162 | 262* | 302
116 | 162 | 166 | 286 | 312
118 | 166 | 170 | 294 | 318
120 | 170 | 174 | 298 | 324

Figure 11: The $Q_{\alpha}$ values for even-even actinides and superheavy
nuclei calculated with indicated CEDFs.

Table 3 compares the positions of the two-proton and two-neutron drip lines
obtained in the RHB calculations with the NL3* and PC-PK1 functionals; the
results for the DD-PC1 and DD-ME2 functionals could be found in Table IV of
Ref. [39]. One can see that the two-proton drip lines are very similar in
these two functionals; they differ by no more than four neutrons. This is in
line with earlier observations that theoretical uncertainties in the
predictions of the position of two-proton drip line are relatively small (see
Ref. [93] and Sec. VII in Ref. [39]). Note that two-proton drip lines in the
isotopic chains of interest obtained with PC-PK1 are very close to those
obtained with DD-PC1 (compare Table 3 in the present article with Table IV of
Ref. [39]).

Among the considered CEDFs the PC-PK1 functional provides the most neutron-
rich two-neutron drip line and the NL3* provides the second most neutron-rich
two-neutron drip lines (compare Table 3 in the present article with Table IV
of Ref. [39] and see Sec. VIII in Ref. [39]). All employed functionals reveal
the presence of the shell closure at $N=258$ (see Fig. 6d in Ref. [90]). The
size of this gap is almost the same in the NL3* and PC-PK1 functionals, but
(i) it is shifted down in energy by $\approx 400$ keV for PC-PK1 as compared
with NL3* and (ii) high-$j$ intruder orbitals $1k_{15/2}$ and $2i_{13/2}$,
which have a significant impact on the position of neutron drip line (see
discussion in Ref. [90]), are located at lower energies in the PC-PK1
functional as compared with the NL3* one. These features lead to the shift of
the two-neutron drip line to substantially higher neutron numbers in the PC-
PK1 CEDF as compared with NL3*. The sizes of the $N=258$ shell gaps are
smaller by $\approx 10$% and 20% in the DD-ME2 and DD-PC1 functionals as
compared with the ones in the PC-PK1 and NL3*. In addition, above mentioned
high-$j$ intruder orbitals in the calculations with the DD-ME2 and DD-PC1
functionals are located at higher energies as compared with the ones in NL3*.
As a consequence, their two-neutron drip lines are located at lower neutron
numbers as compared with the NL3*. These features are clearly seen in Fig. 5.

Fig. 8 shows the map of theoretical uncertainties $\Delta E(Z,N)$ in binding
energies. These uncertainties increase drastically when approaching the
neutron-drip line and in some nuclei they reach 50 MeV. Poorly defined
isovector properties of CEDFs is the major reason for that (see Ref. [39]).
Note that the $\Delta E(Z,N)$ spreads for the NL3*, DD-PC1 and DD-ME2
functionals are relatively modest [see Fig. 8(b)]101010The addition of the DD-
ME$\delta$ functional to this set of three functionals is not expected to
modify significantly $\Delta E(Z,N)$ (see Fig. 9 in Ref. [39]). and the major
contribution to $\Delta E(Z,N)$ is coming from the PC-PK1 functional (compare
panels (a) and (b) in Fig. 8). The fact that isovector properties of the PC-
PK1 functional are significanly different from those of NL3*, DD-PC1, DD-ME2
and DD-ME$\delta$ is also confirmed by the analysis of binding energies in the
Yb ($Z=70$) isotopic chain (see Fig. 3 in Ref. [49]). As follows from the
analysis of parametric correlations in different classes of CEDFs performed in
Ref. [58], a possible reason for that could be related to over-parametrization
of the isoscalar channel in this class of CEDFs111111The analysis of Ref. [58]
suggests that the number of parameters in the isoscalar channel of PC-PK1 CEDF
can be reduced from 4 to 1.. This, in turn, may lead to a somewhat wrong
balance of the isoscalar and isovector channels in known nuclei which reveals
itself in a more pronounced way via different (as compared with other
functionals) isovector dependence of binding energies in neutron-rich one.

Figure 9 presents the summary of two-neutron separation energies $S_{2n}(Z,N)$
obtained with four employed CEDFs. Note that some discontinuities in smooth
trends of the $S_{2n}(Z,N)$ distributions as a function of neutron number are
either due to the presence of substantial spherical shell gaps at $N=184$ or
$N=258$ or due to the crossing of the boundaries between prolate and oblate
shapes. For example, the impact of the $N=184$ spherical shell gap on the
$S_{2n}(Z,N)$ distributions is clearly visible in Figs. 9(b), (c) and (d) (see
also Figs. 5 (b), (c) and (d) for deformation distributions). On the contrary,
its impact is substantially suppressed in superheavy nuclei in the
calculations with CEDF DD-PC1 [see Fig. 9(a)] because of the reduced role of
the $N=184$ spherical shell gap in this functional.

Finally, the spreads $\Delta S_{2n}(Z,N)$ in two-neutron separation energies
are presented in Fig. 10. They are the lowest in known nuclei but in general
increase with increasing neutron number. The $\Delta S_{2n}(Z,N)$ values are
quite large ($\Delta S_{2n}(Z,N)\approx 2.2$ MeV) in the vicinity of two-
neutron drip lines and the $N=184$ and $N=258$ spherical shell gaps. However,
they become extremely large ($\Delta S_{2n}(Z,N)\approx 4.0$ MeV) at the
boundaries between prolate and oblate shapes. Similar to the spreads in
binding energies (see discussion of Fig. 8 above), the largest contribution to
the $\Delta S_{2n}(Z,N)$ values comes from the CEDF PC-PK1. If the PC-PK1
functional is excluded from consideration these values on average decrease by
a factor of 2 (compare panels (b) and (a) in Fig. 10). It is interesting that
in neutron-rich deformed $N\approx 190-236$ region the $\Delta S_{2n}(Z,N)$
values are on average comparable with those in known nuclei (see Fig. 10(b).
However, they still show increased magnitudes at above discussed locations of
nuclear chart.

##  IV $\alpha$-decay properties

In actinides and superheavy nuclei spontaneous fission and $\alpha$ emission
compete and the shortest half-life determines the dominant decay channel and
the total half-life. Only in the cases where the spontaneous fission half-life
is longer than the half-life of $\alpha$ emission can superheavy nuclei be
observed in experiment. In addition, only nuclei with half-lives longer than
$\tau=10$ $\mu$s are observed in experiments.

The $\alpha$ decay half-live depends on the $Q_{\alpha}$ values which are
calculated according to

$Q_{\alpha}=E(Z,N)-E(Z-2,N-2)-E(2,2)$ (18)

with $E(2,2)=-28.295674$ MeV [94] and $Z$ and $N$ representing the parent
nucleus.

Figure 12: The $Q_{\alpha}$ spreads $\Delta Q_{\alpha}(Z,N)$ as a function of
proton and neutron number. $\Delta Q_{\alpha}(Z,N)=|Q_{\alpha}^{\rm
max}(Z,N)-Q_{\alpha}^{\rm min}(Z,N)|$, where $Q_{\alpha}^{\rm max}(Z,N)$ and
$Q_{\alpha}^{\rm min}(Z,N)$ are the largest and smallest $Q_{\alpha}$ values
obtained with four employed CEDFs for the $(Z,N)$ nucleus.

To estimate theoretical uncertainties in the predictions of the $\alpha$-decay
half-lives, they were calculated using three phenomenological expressions,
namely,

  * •

the Viola-Seaborg semiempirical formula [95]

$log_{10}\tau_{\alpha}=\frac{aZ+b}{\sqrt{Q_{\alpha}}}+cZ+d$ (19)

employing two sets of parametrizations. The first one with the parameters
$a=1.66175$, $b=-8.5166$, $c=-0.20228$ and $d=-33.9069$ has been fitted in
Ref. [96]. This set and the results obtained with it are labelled further as
VSS-1989. Another set has been defined in Ref. [97] and its paramaters are:
$a=1.64062$, $b=-8.54399$, $c=-0.19430$ and $d=-33.9054$. The label VSS-2005
is used for it and its results.

  * •

phenomenological first modified Brown fit (mB1) [98]

$log_{10}\tau_{\alpha}=\frac{a(Z-2)^{b}}{\sqrt{Q_{\alpha}}}+c$ (20)

with the parameters $a=13.0705$, $b=0.5182$, $c=-47.8867$. This set and its
results are labeled further as MB-2016.

  * •

phenomenological Royer model [99]

$log_{10}\tau_{\alpha}=\frac{aZ}{\sqrt{Q_{\alpha}}}+b{A}^{\frac{1}{6}}{\sqrt{Z}}+c$
(21)

with the parameters $a=1.5864$, $b=-1.1629$ and $c=-25.31$ of Ref. [99]. Its
results are labeled further as Royer-2000.

Figure 13:  Calculated $log_{10}(\tau_{\alpha})$ values of the $\alpha$-decays
for even-even superheavy nuclei obtained with the VSS-2005 version of Viola-
Seaborg semi-empirical formula for four indicated CEDFs. Open squares are used
for the nuclei in which $\alpha$-decay is energetically forbidden. The white
line corresponding to $log_{10}(\tau_{\alpha})=1.0$ outlines the region of
nuclei in which the alpha-decay half-live is smaller than 10 s.

These phenomenological expressions employ different functional dependencies
(in particular, they show different dependencies on proton and mass numbers)
and are fitted to different sets of experimental data. This is expected to
lead to different predictions for $\tau_{\alpha}$ in high-$Z$ and neutron-rich
nuclei.

Figure 14: The $\tau_{\alpha}$ spreads $\Delta\tau_{\alpha}(Z,N)$ as a
function of proton and neutron number.
$\Delta\tau_{\alpha}(Z,N)=|\tau_{\alpha}^{\rm max}(Z,N)-\tau_{\alpha}^{\rm
min}(Z,N)|$, where $\tau_{\alpha}^{\rm max}(Z,N)$ and $\tau_{\alpha}^{\rm
min}(Z,N)$ are the largest and smallest $\tau_{\alpha}$ values obtained with
selected set of functionals and empirical formulas. Panel (a) shows these
spreads obtained with four employed CEDFs and VSS-2005 empirical formula. DD-
PC1 CEDF and four empirical formulas are used in panel (b). Panel (c) presents
combined spread of $\tau_{\alpha}$ obtained with four CEDFs and four empirical
formulas. Panel (d) is the subversion of panel (c) in which the PC-PK1
functional is excluded.

The $Q_{\alpha}$ values calculated with the DD-PC1, DD-ME2, NL3* and PC-PK1
functionals are presented in Fig. 11. One can see that for a fixed value of
$Z$ with increasing neutron number the $Q_{\alpha}$ values in general
decrease. They are positive in proton-rich nuclei as well as in the nuclei
located close to the $\beta$-stability line. The $Q_{\alpha}$ values
experience a substantial increase at shell closure with $N=184$121212Similar
increase is also seen in the vicinity of the $N=258$ shell closure in the
calculations with the NL3* and PC-PK1 CEDFs [see Fig. 11(c) and (d)]. (see
Fig. 11 in the present paper as well as Fig. 14 in Ref. [59]); note that the
effect of this shell closure is washed out in the $Z>110$ nuclei for the DD-
PC1 functional. With subsequent increase of neutron number the $Q_{\alpha}$
values become first smaller, then they become close (or equal) to zero and
with further increase of $N$ they get more and more negative. Note that
$\alpha$-decay is energetically not possible for $Q_{\alpha}\leq 0$ MeV. Thus,
very neutron-rich nuclei cannot decay by $\alpha$-emission.

Figure 15: The heights of primary fission barriers (in MeV) obtained in axial
RS-RHB and RA-RHB calculations as a function of proton and neutron numbers for
nuclei located between two-proton and two-neutron drip lines.

Note that general trends in the development of the $Q_{\alpha}$ values as a
function of proton and neutron number are similar in all functionals. The
major differences are related to the location of the two-neutron drip line,
the impact of the $N=184$ and $N=258$ shell closures and the location of the
boundaries between prolate and oblate nuclear shapes. These differences
between the functionals are summarized in Fig. 12 which shows the $Q_{\alpha}$
spreads $\Delta Q_{\alpha}(Z,N)$ as a function of proton and neutron number.
The largest spread in the predictions exists in the island centered around
$Z\sim 110,N\sim 198$ in which $\Delta Q_{\alpha}(Z,N)>3$ MeV. This spread
emerges from different predictions of the boundaries in the $(Z,N)$ plane
between prolate and oblate shapes (see Fig. 5) and coincides with the largest
spread in calculated ground state deformations (see Fig. 7). The next region
with the largest differences in the predictions is located between neutron
numbers N=236 and N=258 (see Fig. 12). However, these differences are not
critical because (i) this region is not expected to play a role in the
r-process, (b) expected $\alpha$-decay half-lives exceed 1020 s (see Fig. 13)
and (c) many of the nuclei in this region are not expected to decay by
$\alpha$-emission. High $\Delta Q_{\alpha}$ values ($\Delta
Q_{\alpha}(Z,N)\approx 1.5$ MeV) are observed near shell closure at $N=184$
and in very neutron-rich nuclei near two-neutron drip line. This is a
consequence of the difference in the predictions of the ground state
properties such as deformations in the nuclei near $N=184$ (see Ref. [59]) and
general deterioration of predictive power of nuclear models on approaching
neutron drip line (see Ref. [39]). In other regions of nuclear chart, $\Delta
Q_{\alpha}(Z,N)\leq 1.0$ MeV with smallest spreads seen in the $N<180$ nuclei.

Note that the inclusion of dynamical correlations (for example, by means of 5
dimensional collective Hamiltonian) can locally modify the binding energies
and $Q_{\alpha}$ values [63, 68, 86] but they have the largest impact on
transitional nuclei which represent only minor part of the nuclei under study.
For well deformed nuclei, the impact of dynamical correlations on $Q_{\alpha}$
values is rather modest [63]. Thus, their inclusion is not expected to change
drastically global picture for the behavior of $Q_{\alpha}$.

Calculated $\alpha$-decay half-lives $\tau_{\alpha}$ (in logarithmic scale)
obtained with the VSS-2005 empirical formula for four CEDFs are shown in Fig.
13. Other phenomenological formulas such as VSS-1989, MB-2016 and Royer-2000
give similar results; thus, they are not shown. For a given isotope chain the
calculated half-lives generally increase with increasing neutron number. This
trend is interrupted only at the $N=184$ and $N=258$ shell closures. The
consequence of this feature is the fact that traditional experimental
technique of detecting superheavy nuclei by $\alpha$-decay will not work in
neutron-rich nuclei because they can decay faster by spontaneous fission. Note
that $\alpha$-decay is energetically forbidden for a large group of very
neutron-rich nuclei located in the vicinity of neutron-drip line. In such
nuclei as well as in those which have very large $\tau_{\alpha}$ values, the
competition of spontaneous fission, neutron induced fission, $\beta$-decay,
and neutron emission will define the leading channel of decay in the r-process
calculations.

Figure 16: The distributions of the nuclei, in which the outer fission barrier
is higher than inner one in the RS-RHB calculations, in the $(Z,N)$ plane for
four employed CEDFs. Different types of symbols are used to indicate the
impact of octupole deformation on outer fission barriers of these nuclei. Blue
solid circles are used for the nuclei the heights of outer fission barriers of
which are not affected by the inclusion of octupole deformation. Solid red
circles/orange squares are used for the nuclei in which outer fission barrier
is affected by octupole deformation in the RA-RHB calculations and is
higher/lower than inner fission barrier.

It is important to remember that the typical timescale of the r-process is in
the seconds range [3, 2, 4]. Fig. 13 clearly illustrates that with few
exceptions the nuclei located to the right of white lines have $\alpha$-decay
half-lives exceeding 10 s. Thus, $\alpha$-decay half-lives of these nuclei are
longer than the typical timescale of the r-process and, as a consequence,
their alpha-decays will not affect the r-process simulations. These white
lines in Fig. 13 also outline the region of nuclear chart in which traditional
experimental measurements of superheavy nuclei based on the $\alpha$-decays
are possible: these are the regions located near and to the left of these
white lines.

Theoretical uncertainties in the predictions of $\alpha$-decay half-lives
given via the $\Delta(log_{10}(\tau_{\alpha}))$ spreads are summarized in Fig.
14. The comparison of the panels (a) and (b) clearly shows that these
uncertainties mostly emerge from the differences in the predictions of the
$Q_{\alpha}$ values by different functionals. These uncertainties exceed 50
orders of magnitude in the nuclei located in the vicinity of two-neutron drip
line and in some nuclei around $Z\approx 108,N\approx 198$ (see Fig. 14). The
uncertainties in $\tau_{\alpha}$ originating from different empirical formalas
[see Eqs. (19), (20) and (21)] are significantly smaller [see Fig. 14(b)]. For
absolute majority of the nuclei they are smaller than 5 orders of magnitude
and for proton-rich nuclei and the nuclei located close to the beta-stability
they are very small. They are larger than 10 orders of magnitude only in
neutron-rich nuclei located in close vicinity of two-neutron drip line.
Combined theoretical uncertainties in $\tau_{\alpha}$ emerging from the use of
four empirical formulas and four CEDFs are summarized in Fig. 14(c). One can
see that for almost half of nuclei they exceed 10 orders of magnitude; these
nuclei are located on neutron-rich side of nuclear chart. However, these
uncertainties are not very critical since the $\alpha$-decay lifetimes become
extremely large in such nuclei (see Fig. 13) so $\alpha$-decay can compete
neither with fission nor with $\beta$-decay. Note also that the removing of
the PC-PK1 functional from consideration does not change appreciably
theoretical uncertainties in the predictions of $\alpha$-decay half-lives
(compare Figs. 14 (c) and (d)).

##  V Fission properties

###  V.1 Primary fission barriers

The distributions of primary fission barriers131313The highest in energy
fission barrier (among inner and outer ones) is called primary and it plays an
important role in the r-process modeling (see Ref. [85]). (PFB) heights in the
$(Z,N)$ plane obtained with employed functionals are shown in Fig. 15. Fig. 16
presents the maps of the nuclei in the region under study in which outer
fission barriers are higher than inner ones in the RS-RHB calculations (see
Figs. 3 and 4 for more details). It also illustrates that the importance of
outer fission barriers in stabilization of nuclei in general decreases on
going from light actinides to superheavy nuclei because of increased
importance of Coulomb interaction (compare Figs. 4 and 3). Fig. 16 also
demonstrates the impact of octupole deformation (as obtained in RA-RHB
calculations) on the outer fission barriers and on their heights with respect
of inner ones. The lowering of outer fission barrier due to octupole
deformation indicates that asymmetric fission becomes dominant, while the
absence of the impact of octupole deformation on outer fission barrier height
tells that fission will be symmetric.

Fig. 16 shows that similar regions in the $(Z,N)$ plane, in which the outer
fission barriers are higher in energy than inner ones in the RS-RHB
calculations, appear in the calculations with all employed functionals.
However, these regions are substantially larger in the density dependent
functionals (such as DD-PC1 and DD-ME2) as compared with CEDFs NL3* and PC-
PK1. Octupole deformation does not affect outer fission barriers in the nuclei
located in the $Z\sim 110,N\sim 240$ region, in the $N\leq 120$ nuclei as well
as in the nuclei located not so far away from $N\approx 180$. On going away
from the latter two regions, octupole deformation first starts to reduce the
heights of outer fission barriers but they still remain higher in energy than
inner ones. Further transition away from these regions leads to the reduction
of the heights of outer fission barriers below the inner ones due to the
impact of octupole deformation.

Fig. 15 reveals a lot of similarities in the predictions of the global
structure of the maps of fission barrier heights obtained with employed
functionals. The highest PFBs are predicted in the islands of low-$Z$ nuclei
centered around spherical shell closures with $N=126$ and $N=184$ (and $N=258$
in the case of the PC-PK1 functional). Fission barriers reach 15 MeV in the
centers of these islands. Next island with high fission barriers exists around
$Z\approx 100,N\approx 150$. Left bottom part of this island coincides with
the region of actinides (see, for example, Fig. 7 in Ref. [62]) in which the
heights of fission barriers have been experimentally measured. Relativistic
mean field calculations with the NL3*, PC-PK1 and DD-PC1 functionals performed
by different groups rather well describe inner and outer fission barriers in
actinides [62, 63, 69, 64]. Note also that the spreads of the heights of inner
fission barriers obtained with these three functionals in the $Z\approx
100,N\approx 150$ island are relatively small (approximately 1 MeV or less)
for the majority of nuclei in this island (see Fig. 3b in Ref. [47]). On the
contrary, the DD-ME2 functional predicts somewhat higher fission barriers in
this island (see Fig. 15b) which leads to somewhat higher spreads $\Delta
E^{B}$ in the heights of primary fission barriers (see Fig. 17a). The island
of low fission barriers is seen near $Z\approx 108,N\approx 192$ in all
functionals. Then another island of high fission barriers centered around
$Z\approx 104,N\approx 216$ is formed. The highest fission barriers reaching
$10-11$ MeV in the center of this island are predicted by the DD-ME2
functional (see Fig. 15b). Somewhat lower fission barriers (with approximately
9 MeV height in the center of the island) are predicted by the DD-PC1
functional (see Fig. 15a). Fission barriers with height of approximately 6 MeV
appear in broad region of this island in the calculations with the NL3* and
PC-PK1 CEDFs (see Figs. 15c and d).

Figure 17: (a) The spreads $\Delta E^{B}$ of the heights of primary fission
barriers as a function of proton and neutron numbers. $\Delta
E^{B}(Z,N)=|E^{B}_{max}(Z,N)-E^{B}_{min}(Z,N)|$, where, for given $Z$ and $N$
values, $E^{B}_{max}(Z,N)$ and $E^{B}_{min}(Z,N)$ are the largest and smallest
heights of inner fission barriers obtained with employed set of four
functionals. (b-d) The spreads $\Delta E^{B}$ obtained for indicated pairs of
the functionals.

With increasing neutron number the predictions start to differ substantially.
The NL3* and PC-PK1 functionals predict extremely low fission barriers with
heights of around 2 MeV or less for the band of nuclei around $N\approx 240$
(see Figs. 15c and d). No such band is formed in the calculations with DD-ME2
and DD-PC1 functionals (see Figs. 15a and b). This could have a drastic impact
on the creation of superheavy elements in the r-process because the nuclear
flow during most of neutron irradiation step of the r-process follows the
neutron drip line and produces in tens of ms the heaviest drip line nuclei
(see the discussion in Sec. 4 of Ref. [12]). However, this nuclear flow will
most likely be terminated at $N\approx 240$ nuclei in the calculations with
NL3* and PC-PK1 since fission will be much faster than neutron capture. Thus,
no superheavy nuclei are expected to be formed in the r-process calculations
based on fission barriers obtained with these two functionals. This is similar
to the results of the r-process simulations based on non-relativistic models
(such as Skyrme DFT with HFB-14 functional, Thomas-Fermi (TF) model and Finite
Range Droplet Model (FRDM)) with low fission barriers in the vicinity of
neutron drip line [15, 12]. On the contrary, the formation of superheavy
elements in the r-process is more likely in the calculations based on the DD-
ME2 and DD-PC1 functionals since the ($Z,N$) region near neutron drip line is
characterized by relatively high fission barriers (see Figs. 15a and b). As a
consequence, neutron capture is expected to proceed faster than fission and
nuclear flow during neutron irradiation step of the r-process could extend to
higher proton numbers [15, 12]. To a degree this is similar to the r-process
simulations based on the models (such as Extended Thomas Fermi (ETFSI) model)
with high fission barriers near neutron drip line (see Ref. [15]).

Further increase of neutron number leads to a rise of the heights of primary
fission barriers (PFBs) and the formation of the band of nuclei near shell
closure at $N=258$ with the heights of PFBs exceeding 6 MeV (see Fig. 15b, c
and d). In some nuclei they even reach 12-15 MeV in the calculations with PC-
PK1 and NL3* CEDFs (see Fig. 15c and d). Note that this process is suppressed
in the calculations with CEDF DD-PC1 (see Fig. 15a) because of reduced impact
of the $N=258$ spherical shell gap on the ground state deformations (see Fig.
5a). Nuclear landscape extends substantially beyond $N=258$ in the
calculations with NL3* and PC-PK1 (see Fig. 15c and d). In this region we
again see the alteration of the regions of low (near $N\approx 268$) and high
(near $N\approx 280$ and above) PFBs. Note that in the $N>258$ region toroidal
shapes [100, 65] could become the lowest in energy solutions. This has been
verified in the calculations with DD-PC1 functional in Ref. [66] but has not
been checked for other functionals.

###  V.2 Theoretical uncertainties in primary fission barriers and their
sources

The spreads $\Delta E^{B}$ of the heights of primary fission barriers as a
function of proton and neutron numbers for four employed functionals are shown
in Fig. 17(a). One can see that on the average they are moderate (on the level
of 1.0-1.5 MeV) in the neutron poor region of actinides centered around
$(Z\sim 100,N\sim 140$). Then the spreads start to increase with a small
island of high $\Delta E^{B}\approx 4.0$ MeV values seen in superheavy nuclei
around $Z\approx 110,N\approx 164$. Further increase of neutron number leads
to the band of high $\Delta E^{B}\approx 4.0$ MeV values for the nuclei with
$N\approx 184$. The sources of these uncertainties are directly related to the
differences in the predictions of the ground state properties of these nuclei:
spherical ground states are predicted in these nuclei in the calculations with
PCPK1, DD-ME2 and NL3* functionals [see Figs. 5(b),(c) and (d)] while the DD-
PC1 functional prefers oblate ground states in superheavy $N\approx 184$
nuclei (see Fig. 5(a)). After crossing this band, theoretical uncertainties in
the heights of primary fission barriers substantially decrease and for the
majority of the nuclei located inside the triangle with the sides $N=188$ [for
$Z=96-120$], $Z=120$ [for $N=188-240$] and [$Z=96,N=188$ to $Z=120,N=240$]
they are in general better than 2.5 MeV and in many nuclei they are even
better than 1.5 MeV. However, the $\Delta E^{B}$ spreads start to increase
again on approaching two-neutron drip line. Here they form quite wide band of
the nuclei parallel to two-neutron drip line in which $\Delta E^{B}$ is close
to 4.0 MeV. Even higher spreads reaching 5.5 MeV are seen near shell closure
at $N=258$.

The analysis of the spreads $\Delta E^{B}$ allows to identify major sources of
theoretical uncertainties in the predictions of the heights of PFBs. These
could be reduces to two major contributors, namely, underlying single-particle
structure mostly affecting the ground state properties and nuclear matter
properties of employed CEDFs. To facilitate the discussion we will consider
the $\Delta E^{B}$ spreads for the pairs of selected functionals.

The lowest spreads exist for the pair of the NL3* and PC-PK1 functionals (see
Fig. 17c): $\Delta E^{B}\leq 0.5$ MeV for absolute majority of the nuclei and
only in specific regions of nuclear chart it is higher. Even in those regions
it is higher than 1 MeV only for limited set of the nuclei. These regions are:
(i) the actinides around $N=126$, (ii) the $Z=90$ and 92 actinides with
$N\approx 170$, (iii) superheavy nuclei in the vicinity of the $Z=120$ and
$N=184$ lines, (iv) very neutron-rich nuclei in the vicinity of two-neutron
drip line and (v) the band of the nuclei around $N\approx 246$. The
differences in the predictions of the heights of outer fission barriers are
responsible for the spreads in the region (i). At present, their source is not
clear. The spreads $\Delta E^{B}$ seen in the region (ii) are due both to
different proton and neutron dependencies of the impact of octupole
deformation on outer fission barriers in the NL3* and PC-PK1 functionals
(compare panels (c) and (d) of Fig. 16) and the fact that in some nuclei we
compare the heights of outer and inner fission barriers. In the region (iii),
the large $\Delta E^{B}$ values are due to the differences in the predictions
of the spherical shell closures at $Z=120$ and $N=184$ and the densities of
the single-particle states in their vicinities (see Fig. 1 in Ref. [59] and
the discussion in this reference). Slightly different isovector properties of
the NL3* and PC-PK1 functionals (see Table 2) may be responsible for the
divergence of their predictions in the region (iv). The large $\Delta E^{B}$
values in the region (v) are due to prolate-oblate-spherical shape coexistence
which takes place in slightly different regions of the $(Z,N)$ chart in these
two functionals [compare Figs. 5 (c) and (d)].

The comparison of the predictions of the NL3* and PC-PK1 functionals for the
fission barriers is quite illuminating since it shows in a global way a number
of important features. First, apart of the regions (i)-(v) the comparable
(typically within 0.5 MeV) predictions for the heights of PFBs are obtained on
a global scale by these two functionals despite the fact that they differ
substantially in the predictions of the ground states energies in neutron-rich
nuclei (see the discussion of Fig. 8). Thus, the description of the ground
state energies is to a degree decoupled from the description of fission
barriers; the latter depends on the relative energies of the saddle and the
ground state. As a consequence, good description of the ground state energies
does not guarantee good description of the fission barriers and vice versa.
Second, the differences in the predictions of the PFB heights seen in the
regions (iii) and (v) are related to the differences in the predictions of the
ground states properties, which in turn are defined by the differences in the
underlying single-particle structure. Third, if to exclude the regions (i)-(v)
from consideration it becomes clear that some differences in nuclear matter
properties such as the symmetry energy $J$ and its slope $L_{0}$ (see Table 2)
do not lead to important differences in the predictions for PFBs. Fourth,
comparable global predictions for the PFBs are obtained despite underlying
differences in the basic structure of the functionals and their fitting
protocols. The NL3* functional includes meson exchange of finite range, while
PC-PK1 CEDF does not have mesons and thus it has zero range interactions (see
Ref. [58]). The fitting protocol of the CEDF NL3* is based on 12 spherical
nuclei and includes empirical data on nuclear matter properties (see Ref.
[43]), while the one for PC-PK1 includes only data on binding energies (60
spherical nuclei) and charge radii (17 spherical nuclei) [50]. Note that the
NL3* and PC-PK1 functionals have 6 and 9 parameters, respectively. However,
the analysis of parametric correlations shows that in reality there are only 5
and 6 independent parameters in these two functionals [46, 58].

Next we consider the spreads $\Delta E^{B}$ obtained with the NL3*/DD-ME2 pair
of the functionals (see Fig. 17d). These two functionals have almost identical
fitting protocols (see Refs. [43, 42]). The only difference is the fact that
DD-ME2 fitting protocol uses 3 experimental data points on neutron skins as
compared with 4 in NL3* but the impact of this difference is expected to be
small. Thus, larger values of the $\Delta E^{B}$ spreads in the NL3*/DD-ME2
pair as compared with the ones in the NL3*/PC-PK1 pair are related to the
basic difference of these two functionals, namely, to the implementation of
density dependence. The DD-ME2 functional has explicit dependence of the
meson-nucleon coupling on the nucleonic density, while NL3* employs cubic and
quartic powers of the $\sigma$ meson for density dependence (see Sect. II of
Ref. [39] for details). In addition, the nuclear matter properties (in
particular, the symmetry energy $J$ and its slope $L_{0}$) of these two
functionals differ substantially (see Table 2) and this difference is expected
to contribute into the increase of the spreads $\Delta E^{B}$ obtained for the
NL3*/DD-ME2 pair as compared with those for the NL3*/PC-PK1 pair.

Finally, the $\Delta E^{B}$ spreads for the DD-PC1/DD-ME2 pair of the
functionals are presented in Fig. 17b. Nuclear matter properties of these two
functionals are close to each other and they are located within the limits of
the SET2b constraint set on the experimental/empirical ranges for the
quantities of interest derived in Ref. [51] (see Table 2). However, fitting
protocols of these two functionals are drastically different: CEDF DD-ME2 is
fitted to the properties of 12 spherical nuclei (see Ref. [43]) while DD-PC1
is defined by the properties of 64 deformed rare-earth nuclei and actinides
(see Ref. [44]). As a result of this difference in fitting protocols, the
largest $\Delta E^{B}$ spreads appear in the vicinity of spherical shell
closures at $N=184$ (with $\Delta E^{B}$ reaching 4.0 MeV) and $N=258$ (with
$\Delta E^{B}$ reaching 5.5 MeV) [see Fig. 17b]. Indeed, the impact of these
shell closures on the equilibrium deformation differs substantially in these
two functionals (compare Figs. 5(a) and (b) and see Fig. 7(b)) and this is a
reason for increased $\Delta E^{B}$ spreads.

###  V.3 The comparison with the results obtained in non-relativistic
calculations

It is interesting to compare the global trends of the heights of PFB in the
$(Z,N)$ plane obtained in the RHB calculations (see Fig. 15) with those
obtained in earlier non-relativistic calculations for which the maps similar
to those presented in Fig. 15 are available. Note that similar to our
calculations all these non-relativistic calculations have been performed only
for axial nuclear shapes.

Fission barriers obtained in Gogny DFT calculations with D1M* functional are
presented in Fig. 12 of Ref. [23]. These calculations cover the region from
two-proton drip line up to the nuclei with two-neutron separation energies of
$S_{2n}=4.0$ MeV. In these calculations, the fission barriers of the $N\leq
186$ nuclei typically exceed 6 MeV and in a number of these nuclei their
heights are close to 12 MeV. Then fission barriers in the $N\approx 190-210$
nuclei are lower than 4 MeV but they increase to approximately 8 MeV on
approaching $S_{2n}=4.0$ MeV line. The differences in the predictions of the
heights of PFB obtained in the CDFT and Gogny DFT calculations (compare Fig.
15 with Fig. 12 of Ref. [23]) are in part related to the differences in the
predictions of ground state properties (compare Fig. 5 in the present paper
with Fig. 5 in Ref. [23]).

Our results for fission barriers (Fig. 15) could also be compared with those
obtained in non-relativistic DFTs with the BCPM and HFB14 functionals and
FRLDM (see Fig. 7 in Ref. [18]). The calculations with HFB14 predict very low
fission barriers (with $E^{B}<4$ MeV) for the $Z\geq 110$ nuclei with
exceptionally low fission barriers ($E^{B}<2$ MeV) in many nuclei located in
the $N\approx 184-210$ region (see middle panel of Fig. 7 in Ref. [18] and
Fig. 13 in Ref. [85]). The RHB calculations predict in general higher fission
barriers (as compared with HFB14 ones), but similar island of low fission
barriers is seen near $Z\approx 108,N\approx 192$ in all functionals (see Fig.
15). However, this island is narrower as compared with the HFB14 one. Fission
barriers obtained with the BCPM functional and FRLDM are somewhat higher than
those obtained with HFB14 (compare top and bottom panel with middle panel of
Fig. 7 in Ref. [18]). However, they share the same general structure in the
$(Z,N)$ plane.

Fission barriers calculated in the DFT framework with Skyrme SLy6, SkI3, SV-
min and SV-bas functionals are presented in Fig. 5 of Ref. [21].
Unfortunately, the colormap used in this figure does not allow to extract the
details in the most interesting energy range of $6-10$ MeV141414Better
colormap for the fission barrier height distribution in the $(Z,N)$ plane
obtained with the Skyrme SV-min functional is used in Fig. 4 of Ref. [22]..
However, the region of low fission barriers (with $E^{B}<4$ MeV) similar to
that discussed above appear in all functionals for $N\approx 190-210,Z\approx
94-120$. Fission barriers obtained in the TF and ETFSI models for the
$(Z=84-120,N=140-236)$ and $(Z=84-115,N=140-216)$ regions of nuclear chart are
presented in Fig. 2 of Ref. [15]. Both of these models show the island of low
fission barriers centered around $Z\approx 110,N\approx 192$. In general, the
ETFSI results are close to above mentioned results obtained with Skyrme EDFs.

##  VI Conclusions

The systematic investigation of the ground state and fission properties of
even-even actinides and superheavy nuclei with $Z=90-120$ from the two-proton
up to two-neutron drip lines has been performed for the first time in the
framework of covariant density functional theory. Four state-of-the-art CEDFs
such as DD-PC1, DD-ME2, NL3* and PC-PK1 are used in this study. They represent
the major classes of the CDFT models which differ by basic assumptions and
fitting protocols. This allows a proper assessment of systematic theoretical
uncertainties for physical observables of interest. Obtained results provide a
necessary theoretical input for the r-process modeling in heavy nuclei and, in
particular, for the study of fission cycling. The main results can be
summarized as follows:

  * •

Quadrupole deformations of calculated ground states and related theoretical
uncertainties have been investigated. It turns out that four employed
functionals predict very similar deformations for the majority of the nuclei.
However, large theoretical uncertainties in quadrupole deformation exist for
some nuclei but they are well localized in the $(Z,N)$ plane. These
uncertainties are mostly due to the uncertainties in the predictions of the
underlying single-particle structure. They are dominated by the uncertainties
in the predictions of both spherical shell closures at $N=184$, $N=258$ and
$Z=120$ (for the $N<190$ nuclei) and deformed single-particle structures
leading to somewhat different boundaries in the $(Z,N)$ plane between the
regions with oblate and prolate shapes. The differences in nuclear matter
properties of employed functionals have only marginal impact on theoretical
uncertainties related to calculated quadrupole deformations.

  * •

Theoretical uncertainties $\Delta E(Z,N)$ in binding energies, emerging from
poorly defined isovector properties of CEDFs, increase drastically when
approaching the neutron drip line and in some nuclei they reach 50 MeV.
However, they reduce substantially (down to maximum value of $\Delta
E(Z,N)\approx 21$ MeV) when PC-PK1 functional is removed from consideration.
Two-neutron drip line of this functional is also located at substantially
higher neutron numbers as compared with the ones obtained with other
functionals. In addition, this functional is the major source of theoretical
uncertainties in two-neutron separation energies. Parametric correlations
leading to an over-parametrization of the isoscalar channel is a possible
reason for such a unique behaviour of the CEDF PC-PK1. Theoretical
uncertainties in two-neutron separation energies reveal clear importance of
the uncertainties in the $N=184$ and $N=258$ spherical shell closures and in
the location of the boundaries between the regions of prolate and oblate
shapes.

  * •

$\alpha$-decay properties, such as the $Q_{\alpha}$ values and the lifetimes
$\tau_{\alpha}$, and related theoretical uncertainties have been investigated
employing four empirical formulas and four CEDFs. While the predictive power
of the models is relatively high on proton-rich side of nuclear chart, it
starts to deteriorate on approaching $N=184$. It is especially low in the
nuclei around $Z\approx 108,N\approx 198$ [the region of the transition from
prolate to oblate ground states] and in very neutron-rich nuclei located in
the vicinity of two-proton drip line. However, the uncertainties in the latter
region are not very important since the $\alpha$-decay in these nuclei is not
expected to play any role in the r-process because of extremely large
lifetimes.

  * •

The distributions of the primary fission barriers in the $(Z,N)$ plane have
been investigated with four employed CEDFs. Globally, the highest fission
barriers are produced by the CEDF DD-ME2, and the lowest ones by the NL3* and
PC-PK1 functionals. The results obtained with DD-PC1 are located between these
two extremes but closer to the DD-ME2 ones. The presence of the band of nuclei
with $N\approx 240$ in the $(Z,N)$ plane with low fission barriers, obtained
in the calculations with the NL3* and PC-PK1 functionals, could have a drastic
impact on the creation of superheavy elements in the r-process. The nuclear
flow during most of neutron irradiation step of the r-process follows the
neutron drip line and produces in tens of ms the heaviest drip line nuclei.
However, this nuclear flow will most likely be terminated at $N\approx 240$
nuclei since fission will be much faster than neutron capture. On the
contrary, the formation of superheavy elements in the r-process is more likely
in the calculations based on the DD-ME2 and DD-PC1 functionals since the
$(Z,N)$ region near neutron drip line is characterized by relatively high
fission barriers and the band of nuclei with low fission barriers (similar to
the one at $N\approx 240$ in the NL3* and PC-PK1 functionals) is absent.

  * •

There are two major sources of theoretical uncertainties in the predictions of
the heights of PFBs, namely, underlying single-particle structure mostly
affecting the ground state properties and nuclear matter properties of
employed CEDFs. For example, the increase of theoretical uncertainties for the
ground states of the nuclei in the vicinity of the $N=184$ and $N=258$
spherical shell closures leads to an increase of theoretical uncertainties for
their fission barriers. The functionals with nuclear matter properties located
in the vicinity of empirical SET2b estimates [51] tend to produce higher
fission barriers as compared with the predictions of the functionals the
nuclear matter properties of which are located outside the limits of the SET2b
constraint set. The problem of finding the best functional for the description
of fission barriers is further complicated by the fact that the description of
the ground state energies is to a degree decoupled from the description of
fission barriers; the latter depends on the relative energies of the saddle
and the ground state. As a consequence, good description of the ground state
energies does not guarantee good description of the fission barriers and vice
versa.

This is first ever systematic attempt within the covariant density functional
theory to provide both the input for the r-process calculations which includes
the ground state and fission properties of actinides and superheavy nuclei and
the assessment of systematic theoretical uncertainties in the physical
quantities of interest. As such it follows the ideology of all previous non-
relativistic calculations of relevance for the r-process of heavy and
superheavy nuclei, which depend also on the fission processes, and assumes the
axial symmetry of nuclei. This is a reasonable approximation for the ground
state properties of the majority of nuclei; the only exception is transitional
nuclei which are soft in $\gamma$-deformation. However, the restriction to
axial symmetry leads to the fact that the calculated inner and outer fission
barriers represent the upper limits and can be potentially lowered when the
triaxiality is taken into account. The r-process simulations with the data
obtained in this study will allow to limit the region of the $(Z,N)$ plane
which has an impact on this process. The hope is that for this limited set of
nuclei, systematic refined calculations taking into account the dynamical
correlations and the triaxiality in the calculations of the part of the
$(\beta_{2},\gamma)$-plane covering ground state, inner fission barrier and
second minimum as well as triaixiality and octupole deformation in the
calculations of the part of the $(\beta_{2},\beta_{3},\gamma)$ plane covering
second minimum, outer fission barrier and the region beyond that will be
possible in the era of exascale computing.

Underlying single-particle structure and nuclear matter properties of CEDFs
emerge as the major sources of theoretical uncertainties. However, they affect
different physical observables in a different way. For example, theoretical
uncertainties in the ground state quadrupole deformations are defined mostly
by the uncertainties in the underlying single-particle structure. On the
contrary, both factors contribute into theoretical uncertainties for fission
barriers. The existence of appreciable theoretical uncertainties in the ground
state and fission properties calls for a better covariant energy density
functionals. The reduction of parametric correlations between the parameters
of CEDFs is one possible way in that direction [58, 55]. In addition,
experimental studies of superheavy elements in the vicinity of the $Z=120$ and
$N=184$ lines, planned at new facilities such as SHE factory [101], will
hopefully provide critical data which will allow to discriminate the
predictions of different models. Such information could be used for a better
constraint of the CEDFs and thus to the reduction of substantial theoretical
uncertainties in this region of nuclear chart which affect all physical
observables of interest and have a direct impact on the modeling of the
r-process.

##  VII Acknowledgement

This material is based upon work supported by the US Department of Energy,
National Nuclear Security Administration under Award No. DE-NA0002925, by the
US Department of Energy, Office of Science, Office of Nuclear Physics under
Award No. DE-SC0013037 and by Ghana Atomic Energy Commission, National Nuclear
Research Institute, Ghana.

##  Appendix A Neglect of dynamical correlations in the fission barrier
calculations

Some of non-relativistic calculations, mentioned in the last paragraph of Sec.
II, take into account dynamical correlations, but others not. Dynamical
correlations are not taken into account in our calculations due to following
reasons.

First, the analysis performed in Ref. [86] in the CDFT-based approach
indicates that in most of the nuclei dynamical correlations modify fission
barrier heights by less than 1 MeV but increase substantially computational
time. The only exceptions are the nuclei with soft potential energy surfaces
the ground state energy minimum of which is located at spherical shape. Note
that the absolute majority of the nuclei under consideration are deformed in
the ground states (see Sect. III). Thus, the errors introduced into the
fission barrier heights due to neglect of dynamical correlations are expected
to be smaller than the ones which are coming from the selection of CEDFs (see
Sect. V.2 below). In addition, theoretical uncertainties in fission barrier
heights defined by their spreads (see Eq. (1)) are not expected to be modified
much by the neglect of dynamical correlations. This is because in the majority
of the cases the topology of potential energy surface of a given nucleus
weakly depends on the underlying functional (see, for example, Figs. 7 and 8
in Ref. [47], Fig. 8 in Ref. [78] and supplemental material to Ref. [86]). As
a consequence, dynamical correlations are expected to be comparable for
different functionals and they will at least partially cancel each other in
Eq. (1).

Second, the inner fission barriers are lowered when the triaxiality is taken
into account (see Ref. [62] and references quoted therein) and the potential
energy surfaces of the ground states in many superheavy nuclei are soft in
$\gamma$-deformation (see Ref. [102] and Appendix B.) In the CDFT
calculations, the outer fission barriers can also be affected by triaxial
deformation151515The important role of triaxiality in the description of outer
fission barriers of actinides has also been discussed in the framework of
microscopic+macroscopic approach in Ref. [103]. via two mechanisms. In the
first one, the saddle of reflection-symmetric triaxial fission path becomes
lower in energy than the saddle of reflection-asymmetric (octupole-deformed)
axial fission path due to underlying shell structure [78]. In the second
mechanism, reflection-symmetric fission path and its saddle loose their axial
symmetry and attain some degree of triaxiality [69, 64]. The investigations of
the impact of dynamical correlations on fission barriers in triaxial
calculations are very rare and quite limited in coverage. In non-relativistic
frameworks, only limited set of actinides [104, 105, 106, 107] and superheavy
[108] nuclei have been studied so far. The impact of dynamical correlations on
fission barriers of restricted set of superheavy nuclei along the $Z=120$
isotopic and $N=174,184$ isotonic chains has been studied in the CDFT-based
framework in Ref. [86].

Third, dynamical calculations do not provide a unique answer because of
underlying assumptions and approximations [104, 106, 109, 36, 18, 110, 23].
For example, there exist substantial differences between the predictions based
on Adiabatic Time Dependent HFB (ATDHFB) and the Generator Coordinate method
(GCM) [based on Gaussian Overlap Approximation (GOA)] schemes [109, 36, 18,
23]. The differences between spontaneous fission half-lives $\tau_{SF}$
obtained in these two schemes could reach many orders of magnitude and
increase with the decrease of the fissility-related parameter $Z^{2}/A$ [which
is equivalent to the increase of neutron number within a given isotopic chain]
(see Fig. 2 in Ref. [23]). Large differences between experimental and
calculated $\tau_{SF}$ also exist; for example, in the U isotopes these
differences reach almost 20 orders of magnitudes when ATDHFB values for
$\tau_{SF}$ are used [23]. As illustrated in Refs. [104, 106] on the example
of the 250,264Fm and 240Pu nuclei, the inclusion of pairing fluctuations
within a least action approach improves the agreement between the predicted
$\tau_{SF}$ values and experiment. However, it remains to be seen whether that
is a general conclusion applicable to all nuclei. In addition, such
calculations are prohibitively expensive (in part, because of breaking of
axial symmetry) and thus are not scalable to global calculations.

In addition, the treatment of the ground state energy $E_{0}$ (which is also
tunneling energy for fission) relies on simplified approximations in the
majority of the publications (see discussion in Ref. [86]). In microscopic
calculations, tunneling energy is associated with the energy of collective
ground state defined either in GCM [111] or in five-dimensional collective
Hamiltonian (5DCH) [86]. The energy of collective ground state depends on
softness (both in quadrupole deformation $\beta_{2}$ and in triaxial
deformation $\gamma$) of collective energy surface in the vicinity of ground
state minimum. It differs from approximate values substantially [86]; this
could modify calculated $\tau_{SF}$ by several orders of magnitude [108, 86].
These extremely large theoretical uncertainties in $\tau_{SF}$, coming from
the selection of the method (ATDHFB versus GCM+GOA) and the treatment of the
ground state energy, are the reasons why we have not attempted to calculate
spontaneous fission half-lives in the present paper.

Fourth, there are some indications that the role of triaxiality can be reduced
in dynamical calculations for some nuclei. For example, it was shown in Ref.
[104] based on least-action calculations with Skyrme EDF SkM* that pairing
fluctuations act in the direction of restoration of axial symmetry along the
fission path in the 240Pu nucleus. This nucleus is characterized by relatively
modest decrease (approximately 2 MeV) of inner fission barrier height by
triaxiality in static calculations. Similar results have also been obtained
for 250,264Fm with similar formalism based on DD-PC1 CEDF in Ref. [106]. The
calculations of Refs. [107, 108] based on least-action principle also indicate
that the axial symmetry of fission pathway is restored in many nuclei; they
are based on the DFT approach with Gogny D1S functional [107] and on
macroscopic+microscopic method [108].

However, not in all nuclei the effect of triaxiality is eliminated by least-
action principle. For example, the least-action fission pathway in 264Fm is
still characterized by triaxiality (although it is somewhat reduced by
enhanced pairing as compared with static fission path) in the calculations
based on Skyrme SkM* functional [104]. Note that in this nucleus the
triaxiality has a large impact (slightly more than 4 MeV) on the height of
inner fission barrier in static calculations. The calculations of Ref. [107]
performed with Gogny D1S functionals also indicate that in some nuclei the
least-action fission pathway goes through triaxial saddles. More systematic
calculations161616These calculations are simplified as compared with quoted
DFT calculations since they use fixed single-particle spectrum for all nuclei
and neglect the deformations of higher order such as $\beta_{6}$ and
$\beta_{8}$. based on macroscopic+microscopic method show that the impact of
the triaxiality on the least action fission pathway (and thus on spontaneous
fission half-live) shows up in some nuclei with $Z=114$ and becomes much more
pronounced in the $Z\geq 120$ nuclei [108]. Note that the tendency towards
restoration of axial symmetry of the fission pathway in the least-action
calculations may somewhat be underestimated in Refs. [108, 107] because of the
neglect of pairing fluctuations.

The analysis of these publications suggests two possible situations in which
the least-action fission path will most likely be characterized by
triaxiality. In the first one, the decrease of the fission barrier by
triaxiality in static calculations is substantial being on the order of $3-4$
MeV [108, 104]. In the second one, the ground state is oblate (or possibly
soft in oblate-prolate direction [86]) so that the fission path across the
$\gamma$-plane is shorter than the one along the $\gamma=0^{\circ}$ line
[108]. As discussed in Sec. III and Appendix B only limited number of nuclei
satisfy such conditions. Thus, the restriction to axial symmetry should be
considered as a reasonable first approximation. However, one should keep in
mind that the values obtained for fission barriers represent upper limits
since their possible lowering due to triaxiality is neglected.

##  Appendix B Possible impact of triaxiality on inner fission barriers

Figure 18: Potential energy surfaces of the Sg ($Z=106$), Ds ($Z=110$), Fl
($Z=114$) and Og ($Z=118$) isotopes with neutron numbers $N=192,202,212$ and
222 obtained in the triaxial RHB calculations with the DD-PC1 functional.
Neutron number is increasing on going from left to right. The energy
difference between two neighboring equipotential lines is equal to 0.5 MeV.
The ground state minima and saddle points are shown by white circles and red
solid squares, respectively.

The restriction to axial symmetry is one of the approximations used in the
present study which is a consequence of the global character of the study (see
detailed discussion presented in the end of Sect. II and in Appendix A). In
order to better understand for which nuclei this approximation may be violated
(even in least-action calculations with pairing fluctuations included such as
those presented in Refs. [104, 106]) we consider the examples of potential
energy surfaces obtained in triaxial RHB calculations with the DD-PC1
functional. These PES calculated for the Sg ($Z=106$), Ds ($Z=110$), Fl
($Z=114$) and Og ($Z=118$) isotopes with neutron numbers $N=192,202,212$ and
222 are presented in Fig. 18. They represent the extension of the
calculations, executed in a more limited deformation space, the results of
which are discussed in Sect. XI of Ref. [66]. The summary of the heights
$E_{triax}^{B}$ of triaxial inner fission barriers and the decreases of the
fission barrier heights due to triaxiality $\Delta E^{gain}$ are presented in
Fig. 21 and Table II of this reference. Note that these superheavy nuclei are
selected in such a way that they cover the part of nuclear chart characterized
by both oblate and prolate ground states (see Fig. 5a).

The review of existing literature presented in Appendix A suggests two
possible scenarios in which the least-action fission path will most likely be
characterized by triaxiality. In the first one, the decrease of the fission
barrier by triaxiality in static calculations is substantial being on the
order of $3-4$ MeV [108, 104]. Such decreases are observed in 328Sg ($\Delta
E^{gain}=4.04$ MeV), 310Og ($\Delta E^{gain}=3.42$ MeV), and 320Og ($\Delta
E^{gain}=4.93$ MeV) (see Figs. 18(d), (m) and (n) and Table II in Ref. [66]).
In the second scenario, the ground state is oblate (or possibly soft in
oblate-prolate direction [86]) so that the fission path across the
$\gamma$-plane is shorter then the one along the $\gamma=0^{\circ}$ line
[108]. This condition is satisfied only in the 306Fl (see Fig. 18(i)) and
310,320Og (see Fig. 18(m) and (n)) nuclei. Based on general features discussed
in Refs. [104, 106], the analysis of PES of remaining nuclei (see Figs. 18 a,
b, c, e, f, g, h, j, k, l, o and p) suggests that least-action fission pathway
will be axial in these nuclei when pairing fluctuations are taken into
account.

Whether one or another scenario takes place depends on the underlying shell
structure (both at the ground state and saddle) defining the topology of
potential energy surfaces in the ($\beta_{2},\gamma$) plane. Fig. 18 shows
that the $N=192$ isotones are extremely soft in the $\gamma$-plane with clear
tendency for the formation of near-oblate triaxial ground state minimum in the
Fl and Og nuclei. However, static fission pathes from these minima are
characterized by low fission barriers so these nuclei are expected to be
unstable. Similar (but slightly less pronounced) situation is also seen for
the $N=202$ isotones. The increase of neutron number to $N=212$ and 222 leads
to a better localization of the ground state minimum at prolate shape and to
an increase of fission barrier heights.

## References

  * Martínez-Pinedo [2008] Martínez-Pinedo, Selected topics in nuclear astrophysics, Eur. Phys. J. Spec. Topics 156, 123 (2008).
  * Liccardo  _et al._ [2018] V. Liccardo, M. Malheiro, M. S. Hussein, B. V. Carlson, and T. Frederico, Nuclear processes in astrophysics: recent progress, Eur. Phys. J. A54, 221 (2018).
  * Kasen  _et al._ [2017] D. Kasen, B. Metzger, J. Barnes, E. Quataert, and E. Ramirez-Ruiz, Origin of the heavy elements in binary neutron-star mergers from a gravitational-wave event, Nature 551, 80 (2017).
  * [4] J. J. Cowan, C. Sneden, J. E. Lawler, A. Aprahamian, M. Wiescher, K. Langangke, G. Martínez-Pinedo, and F.-K. Thielemann, Making the heaviest elements in the universe: A review of the rapid neutron capture process, arkhiv: 1901.014110v1 [astro-ph.HE], subm. to Rev. Mod. Phys. .
  * Thielemann  _et al._ [2017] F. K. Thielemann, M. Eichler, I. V. Panov, and B. Wehmeyer, Neutron star mergers and nucleosynthesis of heavy elements, Ann. Rev. Nucl. Part. Sci. 67, 253 (2017).
  * Abbott  _et al._ [2017] B. P. Abbott _et al._ , Gw170817: Observation of gravitational waves from a binary neutron star inspiral, Phys. Rev. Lett. 119, 161101 (2017).
  * [7] B. P. Abbott _et al._ , Multi-messenger observations of a binary neutron star merger, Astrophys. J. Lett. 848, L12.
  * Just  _et al._ [2015] O. Just, A. Bauswein, R. A. Pulpillo, S. Goriely, and H.-T. Janka, Comprehensive nucleosynthesis analysis for ejecta of compact binary mergers, Mont. Not. R. Astr. Soc. 448, 541 (2015).
  * Sekiguchi  _et al._ [2015] Y. Sekiguchi, K. Kiuchi, K. Kyutoku, and M. Shibata, Dynamical mass ejection from binary neutron star mergers: Radiation-hydrodynamics study in general relativity, Phys. Rev. D 91, 064059 (2015).
  * Mendoza-Temis  _et al._ [2015] J. d. J. Mendoza-Temis, M.-R. Wu, K. Langanke, G. Martínez-Pinedo, A. Bauswein, and H.-T. Janka, Nuclear robustness of the $r$ process in neutron-star mergers, Phys. Rev. C 92, 055805 (2015).
  * Martínez-Pinedo  _et al._ [2007] G. Martínez-Pinedo, D. Mocelj, N. Zinner, A. Kelić, K. Langanke, I. Panov, B. Pfeiffer, T. Rauscher, K.-H. Schmidt, and F.-K. Thielemann, The role of fission in the r-process, Prog. Part. Nucl. Phys. 59, 199 (2007).
  * Goriely [2015] S. Goriely, The fundamental role of fission during r-process nucleosynthesis in neutron star mergers, Eur. Phys. J A51, 22 (2015).
  * Goriely  _et al._ [2013] S. Goriely, J.-L. Sida, J.-F. Lemaître, S. Panebianco, N. Dubray, S. Hilaire, A. Bauswein, and H.-T. Janka, New fission fragment distributions and $r$-process origin of the rare-earth elements, Phys. Rev. Lett. 111, 242502 (2013).
  * Eichler  _et al._ [2015] M. Eichler, A. Arcones, A. Kelic, O. Korobkin, K. Langanke, T. Marketin, G. Martinez-Pinedo, I. Panov, T. Rauscher, S. Rosswog, C. Winteler, N. T. Zinner, and F.-K. Thielemann, The role of fission in neutron star mergers and its impact on the r-process peaks, Astrophys. J. 808, 30 (2015).
  * Petermann  _et al._ [2012] I. Petermann, K. Langanke, G. Martínez-Pinedo, I. Panov, P.-G. Reinhard, and F.-K. Thielemann, Have superheavy elements been produced in nature?, Eur. Phys. J. A 48, 122 (2012).
  * Giuliani  _et al._ [2017] S. A. Giuliani, G. Martínez-Pinedo, L. M. Robledo, and M.-R. Wu, r-process calculations with a microsocpic description of the fission process, Acta Phys. Polonica B 48, 299 (2017).
  * Goriely and Arnould [1996] S. Goriely and M. Arnould, Waiting point approximation and canonical multi-event r-process revisited, Astr. and Astrophys. 312, 327 (1996).
  * 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).
  * Mamdouh  _et al._ [1998] A. Mamdouh, J. Pearson, M. Rayet, and F. Tondeur, Large-scale fission-barrier calculations with the etfsi method, Nucl. Phys. A 644, 389 (1998).
  * Goriely  _et al._ [2007] S. Goriely, M. Samyn, and J. M. Pearson, Further explorations of skyrme-hartree-fock-bogoliubov mass formulas. vii. simultaneous fits to masses and fission barriers, Phys. Rev. C 75, 064312 (2007).
  * Erler  _et al._ [2012] J. Erler, K. Langanke, H. P. Loens, G. Martinez-Pinedo, and P.-G. Reinhard, Fission properties for $r-$process nuclei, Phys. Rev. C 85, 025802 (2012).
  * Reinhard [2018] P.-G. Reinhard, Nuclear density-functional theory and fission of super-heavy elements, Eur. Phys. J 312, 13 (2018).
  * Rodríguez-Guzman  _et al._ [2020] R. Rodríguez-Guzman, Y. M. Humadi, and L. M. Robledo, Microscopic description of fission in superheavy nuclei with the parametrization d1m* of the gogny energy density functional, Eur. Phys. J. A56, 43 (2020).
  * Vretenar  _et al._ [2005] D. Vretenar, A. V. Afanasjev, G. A. Lalazissis, and P. Ring, Relativistic hartree-bogoliubov theory: Static and dynamic aspects of exotic nuclear structure, Phys. Rep. 409, 101 (2005).
  * LNP [2004] Extended density functionals in nuclear structure physics, Extended Density Functionals in Nuclear Structure Physics, Lecture Notes in Physics, edited by G. A. Lalazissis, P. Ring, and D. Vretenar (Springer-Verlag, Heidelberg, 2004) Vol. 641 (2004).
  * Bender  _et al._ [1999] M. Bender, K. Rutz, P.-G. Reinhard, J. A. Maruhn, and W. Greiner, Shell structure of superheavy nuclei in self-consistent mean-field models, Phys. Rev. C 60, 034304 (1999).
  * Litvinova and Afanasjev [2011] E. V. Litvinova and A. V. Afanasjev, Dynamics of nuclear single-particle structure in covariant theory of particle-vibration coupling: From light to superheavy nuclei, Phys. Rev. C 84, 014305 (2011).
  * Koepf and Ring [1989] W. Koepf and P. Ring, A relativistic description of rotating nuclei: the yrast line of 20ne, Nucl. Phys. A 493, 61 (1989).
  * Afanasjev and Abusara [2010a] A. V. Afanasjev and H. Abusara, Time-odd mean fields in covariant density functional theory: nonrotating systems, Phys. Rev. C 81, 014309 (2010a).
  * Hofmann and Ring [1988] U. Hofmann and P. Ring, A new method to calculate magnetic moments of relativistic mean field theories, Phys. Lett. B 214, 307 (1988).
  * Afanasjev and Ring [2000] A. V. Afanasjev and P. Ring, Time-odd mean fields in the rotating frame: Microscopic nature of nuclear magnetism, Phys. Rev. C 62, 031302(R) (2000).
  * Afanasjev and Abusara [2010b] A. V. Afanasjev and H. Abusara, Time-odd mean fields in covariant density functional theory: Rotating systems, Phys. Rev. C 82, 034329 (2010b).
  * Dobaczewski and Dudek [1995] J. Dobaczewski and J. Dudek, Time-odd components in the mean field of rotating superdeformed nuclei, Phys. Rev. C 52, 1827 (1995).
  * Schunck  _et al._ [2010] N. Schunck, J. Dobaczewski, J. McDonnell, J. Moré, W. Nazarewicz, J. Sarich, and M. V. Stoitsov, One-quasiparticle states in the nuclear energy density functional theory, Phys. Rev. C 81, 024316 (2010).
  * Hinohara  _et al._ [2012] N. Hinohara, Z. P. Li, T. Nakatsukasa, T. Nikšić, and D. Vretenar, Effect of time-odd mean fields on inertial parameters of the quadrupole collective hamiltonian, Phys. Rev. C 85, 024323 (2012).
  * 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 (2018).
  * Reinhard and Nazarewicz [2010] P. G. Reinhard and W. Nazarewicz, Information content of a new observable: The case of the nuclear neutron skin, Phys. Rev. C 81, 051303(R) (2010).
  * Dobaczewski  _et al._ [2014] J. Dobaczewski, W. Nazarewicz, and P.-G. Reinhard, Error estimates of theoretical models: a guide, J. Phys. G 41, 074001 (2014).
  * Agbemava  _et al._ [2014] S. E. Agbemava, A. V. Afanasjev, D. Ray, and P. Ring, Global performance of covariant energy density functionals: Ground state observables of even-even nuclei and the estimate of theoretical uncertainties, Phys. Rev. C 89, 054320 (2014).
  * 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 (2003).
  * Boguta and Bodmer [1977] J. Boguta and R. Bodmer, Nucl. Phys. A292, 413 (1977).
  * Lalazissis  _et al._ [2005] G. A. Lalazissis, T. Nikšić, D. Vretenar, and P. Ring, New relativistic mean-field interaction with density-dependent meson-nucleon couplings, Phys. Rev. C 71, 024312 (2005).
  * Lalazissis  _et al._ [2009] G. A. Lalazissis, S. Karatzikos, R. Fossion, D. P. Arteaga, A. V. Afanasjev, and P. Ring, The effective force nl3 revisited, Phys. Lett. B671, 36 (2009).
  * Nikšić  _et al._ [2008] T. Nikšić, D. Vretenar, and P. Ring, Relativistic nuclear energy density functionals: adjusting parameters to binding energies, Phys. Rev. C 78, 034318 (2008).
  * Typel and Wolter [1999] S. Typel and H. H. Wolter, Nucl. Phys. A656, 331 (1999).
  * Agbemava  _et al._ [2019a] S. E. Agbemava, A. V. Afanasjev, and A. Taninah, Propagation of statistical uncertainties in covariant density functional theory: Ground state observables and single-particle properties, Phys. Rev. C 99, 014318 (2019a).
  * Agbemava  _et al._ [2017] S. E. Agbemava, A. V. Afanasjev, D. Ray, and P. Ring, Assessing theoretical uncertainties in fission barriers of superheavy nuclei, Phys. Rev. C 95, 054324 (2017).
  * Agbemava  _et al._ [2016] S. E. Agbemava, A. V. Afanasjev, and P. Ring, Octupole deformation in the ground states of even-even nuclei: a global analysis within the covariant density functional theory, Phys. Rev. C 93, 044304 (2016).
  * Afanasjev and Agbemava [2016] A. V. Afanasjev and S. E. Agbemava, Covariant energy density functionals: Nuclear matter constraints and global ground state properties, Phys. Rev. C 93, 054310 (2016).
  * Zhao  _et al._ [2010] P. W. Zhao, Z. P. Li, J. M. Yao, and J. Meng, New parametrization for the nuclear covariant energy density functional with a point-coupling interaction, Phys. Rev. C 82, 054319 (2010).
  * Dutra  _et al._ [2014] M. Dutra, O. Lourenco, S. S. Avancini, B. V. Carlson, A. Delfino, D. P. Menezes, C. Providencia, S. Typel, and J. R. Stone, Relativistic mean-field hadronic models under nuclear matter constraints, Phys. Rev. C 90, 055203 (2014).
  * Reinhard  _et al._ [1986] P.-G. Reinhard, M. Rufa, J. Maruhn, W. Greiner, and J. Friedrich, Nuclear ground state properties in a relativistic meson field model, Z. Phys. A 323, 13 (1986).
  * Furnstahl  _et al._ [1997] R. Furnstahl, B. D. Serot, and H.-B. Tang, A chiral effective lagrangian for nuclei, Nucl. Phys. A 615, 441 (1997).
  * Roca-Maza  _et al._ [2011] X. Roca-Maza, X. Viñas, M. Centelles, P. Ring, and P. Schuck, Relativistic mean-field interaction with density-dependent meson-nucleon vertices based on microscopical calculations, Phys. Rev. C 84, 054309 (2011).
  * Nikšić  _et al._ [2017] T. Nikšić, M. Imbrišak, and D. Vretenar, “sloppy” nuclear energy density functionals. ii. finite nuclei, Phys. Rev. C 95, 054304 (2017).
  * Kumar  _et al._ [2017] B. Kumar, S. Singh, B. Agrawal, and S. Patra, New parameterization of the effective field theory motivated relativistic mean field model, Nucl. Phys. A 966, 197 (2017).
  * Yüksel  _et al._ [2019] E. Yüksel, T. Marketin, and N. Paar, Optimizing the relativistic energy density functional with nuclear ground state and collective excitation properties, Phys. Rev. C 99, 034318 (2019).
  * Taninah  _et al._ [2020] A. Taninah, S. E. Agbemava, A. V. Afanasjev, and P. Ring, Parametric correlations in energy density functionals, Phys. Lett. B 800, 135065 (2020).
  * Agbemava  _et al._ [2015] S. E. Agbemava, A. V. Afanasjev, T. Nakatsukasa, and P. Ring, Covariant density functional theory: Reexamining the structure of superheavy nuclei, Phys. Rev. C 92, 054310 (2015).
  * Lu  _et al._ [2015] K. Q. Lu, Z. X. Li, Z. P. Li, J. M. Yao, and J. Meng, Global study of beyond-mean-field correlation energies in covariant energy density functional theory using a collective hamiltonian method, Phys. Rev. C 91, 027304 (2015).
  * Agbemava and Afanasjev [2017] S. E. Agbemava and A. V. Afanasjev, Octupole deformation in the ground states of even-even $z\sim 96,n\sim 196$ actinides and superheavy nuclei, Phys. Rev. C 96, 024301 (2017).
  * Abusara  _et al._ [2010] H. Abusara, A. V. Afanasjev, and P. Ring, Fission barriers in actinides in covariant density functional theory: the role of triaxiality, Phys. Rev. C 82, 044303 (2010).
  * Prassa  _et al._ [2012] V. Prassa, T. Nikšić, G. A. Lalazissis, and D. Vretenar, Relativistic energy density functional description of shape transitions in superheavy nuclei, Phys. Rev. C 86, 024317 (2012).
  * Lu  _et al._ [2014] B.-N. Lu, J. Zhao, E.-G. Zhao, and S.-G. Zhou, Multidimensionally-constrained relativistic mean-field models and potential-energy surfaces of actinide nuclei, Phys. Rev. C 89, 014323 (2014).
  * Afanasjev  _et al._ [2018] A. V. Afanasjev, S. E. Agbemava, and A. Gyawali, Hyperheavy nuclei: Existence and stability, Phys. Lett. B 782, 533 (2018).
  * Agbemava  _et al._ [2019b] S. E. Agbemava, A. V. Afanasjev, A. Taninah, and A. Gyawali, Extension of the nuclear landscape to hyperheavy nuclei, Phys. Rev. C 99, 034316 (2019b).
  * Afanasjev  _et al._ [2020] A. V. Afanasjev, S. E. Agbemava, and A. Taninah, Exploring nuclear exotica at the limits, Acta Phys. Polonica B 13, 347 (2020).
  * Zhang  _et al._ [2014] Q. S. Zhang, Z. M. Niu, Z. P. Li, J. M. Yao, and J. Meng, Global dynamical correlation energies in covariant density functional theory: cranking approximation, Frontiers of Physics 9, 529 (2014).
  * Lu  _et al._ [2012] B.-N. Lu, E.-G. Zhao, and S.-G. Zhou, Potential energy surfaces of actinide nuclei from a multidimensional constrained covariant density functional theory: Barrier heights and saddle point shapes, Phys. Rev. C 85, 011301 (2012).
  * Afanasjev  _et al._ [2000] A. V. Afanasjev, P. Ring, and J. König, Cranked relativistic hartree-bogoliubov theory: formalism and application to the superdeformed bands in the $a\sim 190$ region., Nucl. Phys. A676, 196 (2000).
  * Kucharek and Ring [1991] H. Kucharek and P. Ring, Relativistic field theory of superfluidity in nuclei, Z. Phys. A 339, 23 (1991).
  * Ring [1996] P. Ring, Relativistic mean field theory in finite nuclei, Prog. Part. Nucl. Phys. 37, 193 (1996).
  * Karatzikos  _et al._ [2010] S. Karatzikos, A. V. Afanasjev, G. A. Lalazissis, and P. Ring, The fission barriers in actinides and superheavy nuclei in covariant density functional theory, Phys. Lett. B 689, 72 (2010).
  * Tian  _et al._ [2009] Y. Tian, Z. Y. Ma, and P. Ring, A finite range pairing force for density functional theory in superfluid nuclei, Phys. Lett. B 676, 44 (2009).
  * Berger  _et al._ [1991] J. F. Berger, M. Girod, and D. Gogny, Comp. Phys. Comm. 63, 365 (1991).
  * Afanasjev and Abdurazakov [2013] A. V. Afanasjev and O. Abdurazakov, Pairing and rotational properties of actinides and superheavy nuclei in covariant density functional theory, Phys. Rev. C 88, 014320 (2013).
  * Dobaczewski  _et al._ [2015] J. Dobaczewski, A. V. Afanasjev, M. Bender, L. M. Robledo, and Y. Shi, Properties of nuclei in the nobelium region studied within the covariant, skyrme, and gogny energy density functionals, Nucl. Phys. A 944, 388 (2015).
  * Abusara  _et al._ [2012] H. Abusara, A. V. Afanasjev, and P. Ring, Fission barriers in covariant density functional theory: extrapolation to superheavy nuclei., Phys. Rev. C 85, 024314 (2012).
  * Ring and Schuck [1980] P. Ring and P. Schuck, The Nuclear Many-Body Problem (Springer-Verlag, Berlin) (1980).
  * Bonche  _et al._ [2005] P. Bonche, H. Flocard, and P. H. Heenen, Solution of the skyrme hf+bcs equation on a 3d mesh, Comp. Phys. Comm. 171, 49 (2005).
  * Bürvenich  _et al._ [2004] T. Bürvenich, M. Bender, J. A. Maruhn, and P.-G. Reinhard, Systematics of fission barriers in superheavy elements, Phys. Rev. C 69, 014307 (2004).
  * Staszczak  _et al._ [2009] A. Staszczak, A. Baran, J. Dobaczewski, and W. Nazarewicz, Microscopic description of complex nuclear decay: Multimodal fission, Phys. Rev. C 80, 014309 (2009).
  * Möller  _et al._ [2009] P. Möller, A. J. Sierk, T. Ichikawa, A. Iwamoto, R. Bengtsson, H. Uhrenholt, and S. Åberg, Heavy-element fission barriers, Phys. Rev. C 79, 064304 (2009).
  * Myers and Świa¸tecki [1999] W. D. Myers and W. J. Świa¸tecki, Thomas-fermi fission barriers, Phys. Rev. C 60, 014606 (1999).
  * Arnould and Goriely [2020] M. Arnould and S. Goriely, Astronuclear physics: a tale of the atomic nuclei in the skies, arXiv:2001.11228v1 [astro-ph.SR (2020).
  * Shi  _et al._ [2019] Z. Shi, A. V. Afanasjev, Z. P. Li, and J. Meng, Superheavy nuclei in a microscopic collective hamiltonian approach: The impact of beyond-mean-field correlations on ground state and fission properties, Phys. Rev. C 99, 064316 (2019).
  * Sil  _et al._ [2004] T. Sil, S. K. Patra, B. K. Sharma, M. Centelles, and X. Viñas, Superheavy nuclei in a relativistic effective lagrangian model, Phys. Rev. C 69, 044315 (2004).
  * Bender  _et al._ [2001] M. Bender, W. Nazarewicz, and P.-G. Reinhard, Shell stabilization of super- and hyperheavy nuclei without magic gaps, Phys. Lett. B 515, 42 (2001).
  * Zhang  _et al._ [2005] W. Zhang, J. Meng, S. Zhang, L. Geng, and H. Toki, Magic numbers for superheavy nuclei in relativistic continuum hartree-bogoliubov theory, Nucl. Phys. A 753, 106 (2005).
  * Afanasjev  _et al._ [2015] A. V. Afanasjev, S. E. Agbemava, D. Ray, and P. Ring, Neutron drip line: Single-particle degrees of freedom and pairing properties as sources of theoretical uncertainties, Phys. Rev. C 91, 014324 (2015).
  * Eva [2015] Evaluated Nuclear Structure Data File (ENSDF) located at the website (http://www.nndc.bnl.gov/ensdf/) of Brookhaven National Laboratory. ENSDF is based on the publications presented in Nuclear Data Sheets (NDS) which is a standard for evaluated nuclear data. (2015).
  * Xia  _et al._ [2018] X. Xia, Y. Lim, P. Zhao, H. Liang, X. Qu, Y. Chen, H. Liu, L. Zhang, S. Zhang, Y. Kim, and J. Meng, The limits of the nuclear landscape explored by the relativistic continuum hartree–bogoliubov theory, At. Data Nucl. Data Tables 121-122, 1 (2018).
  * Afanasjev  _et al._ [2013] A. V. Afanasjev, S. E. Agbemava, D. Ray, and P. Ring, Nuclear landscape in covariant density functional theory, Phys. Lett. B 726, 680 (2013).
  * Wang  _et al._ [2012] M. Wang, G. Audi, A. H. Wapstra, F. G. Kondev, M. MacCormick, X. Xu, and B. Pfeiffer, Chinese Physics C36, 1603 (2012).
  * V. E. Viola and Seaborg [1966] J. V. E. Viola and G. T. Seaborg, J. Inorg. Nucl. Chem. 28, 741 (1966).
  * Sobiczewski  _et al._ [1989] A. Sobiczewski, Z. Patyk, and S. Cwiok, Deformed superheavy nuclei, Physics Letters B 224, 1 (1989).
  * Dong and Ren [2005] T. Dong and Z. Ren, New calculations of Î±-decay half-lives by the viola-seaborg formula., Eur. Phys. J. A26, 69 (2005).
  * Budaca  _et al._ [2016] A. Budaca, R. Budaca, and I. Silisteanu, Extended systematics of alpha decay half lives for exotic superheavy nuclei, Nucl. Phys. A 951, 60 (2016).
  * Royer [2000] G. Royer, Alpha emission and spontaneous fission through quasi-molecular shapes, J. Phys. G 26, 1149 (2000).
  * Wong [1973] C. Y. Wong, Toroidal and spherical bubble nuclei, Annals of Physics 77, 279 (1973).
  * Oganessian and Dmitriev [2016] Y. T. Oganessian and S. N. Dmitriev, Synthesis and study of properties of superheavy atoms. factory of superheavy elements, Rus. Chem. Rev. 85, 901 (2016).
  * Ćwiok  _et al._ [2005] S. Ćwiok, P.-H. Heenen, and W. Nazarewicz, Shape coexitence and triaxiality in the superheavy nuclei, Nature 433, 705 (2005).
  * Dobrowolski  _et al._ [2007] A. Dobrowolski, K. Pomorski, and J. Bartel, Fission barriers in a macroscopic-microscopic model, Phys. Rev. C 75, 024613 (2007).
  * 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(R) (2014).
  * Benrabia  _et al._ [2017] K. Benrabia, D. E. Medjadi, M. Imadalou, and P. Quentin, Triaxial quadrupole dynamics and the inner fission barrier of some heavy even-even nuclei, Phys. Rev. C 96, 034320 (2017).
  * 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).
  * Delaroche  _et al._ [2006] J.-P. Delaroche, M. Girod, H. Goutte, and J. Libert, Structure properties of even-even actinides at normal and superdeformed shapes analysed using the gogny force, Nucl. Phys. A 771, 103 (2006).
  * Gherghescu  _et al._ [1999] R. Gherghescu, J. Skalski, Z. Patyk, and A. Sobiczewski, Non-axial shapes in spontaneous fission of superheavy nuclei, Nucl. Phys. A 651, 237 (1999).
  * Schunck and Robledo [2016] N. Schunck and L. M. Robledo, Microscopic theory of nuclear fission: A review, Rep. Prog. Phys. 79, 116301 (2016).
  * Giuliani  _et al._ [2019] S. A. Giuliani, Z. Matheson, W. Nazarewicz, E. Olsen, P.-G. Reinhard, J. Sadhukhan, B. Schuetrumpf, N. Schunck, and P. Schwerdtfeger, Colloquium: Superheavy elements: Oganesson and beyond, Rev. Mod. Phys. 91, 011001 (2019).
  * Schindzielorz  _et al._ [1981] N. Schindzielorz, J. Erler, P. Klupfel, P. G. Reinhard, and G. Hager, Int. J. Mod. Phys. E 18, 773 (1981).