Extension of nuclear landscape to hyperheavy nuclei
S. E. Agbemava A. V. Afanasjev A. Taninah A. Gyawali





Department of Physics and Astronomy, Mississippi
State University, MS 39762


Department of Physics and Astronomy, Mississippi
State University, MS 39762
Yukawa Institute of Theoretical Physics, Kyoto University, 
                Japan


Department of Physics and Astronomy, Mississippi
State University, MS 39762


Department of Physics and Astronomy, Mississippi
State University, MS 39762





   The properties of hyperheavy nuclei and the extension of nuclear landscape
to hyperheavy nuclei are extensively studied within covariant density functional
theory. Axial reflection symmetric and reflection asymmetric relativistic 
Hartree-Bogoliubov (RHB) calculations are carried out. The role of triaxiality is 
studied within triaxial RHB and triaxial relativistic mean field + BCS frameworks. 
With increasing proton number beyond Z∼ 130 the transition from ellipsoidal-like 
nuclear shapes to toroidal ones takes place. The description of latter shapes 
requires the basis which is typically significantly larger than the one employed for 
the description of ellipsoidal-like shapes. Many hyperheavy nuclei with toroidal 
shapes are expected to be unstable towards multifragmentation. However, three 
islands of stability of spherical hyperheavy nuclei have been predicted for the first
time in Ref. <cit.>. Proton and neutron densities, charge radii, neutron skins 
and underlying shell structure of the nuclei  located in the centers of these islands 
have  been investigated  in detail. Large neutron shell gaps at N=228, 308 
and 406 define approximate centers of these islands in neutron number. On the
contrary,  large proton gap appear only at Z=154 in the (Z∼ 156, N∼ 310) 
island. As a result, this is the largest island of stability of spherical hyperheavy
nuclei found in the calculations. The calculations indicate the stability of the nuclei in 
these islands with respect of octupole and triaxial distortions. The shape evolution of 
toroidal shapes along the fission path and the stability of  such shapes with respect of 
fission have been studied. Fission barriers in neutron-rich superheavy nuclei  are studied 
in triaxial RHB framework; the impact of triaxiality on the heights of fission barriers is 
substantial in some parts of this region.  Based on the results obtained in the 
present work, the extension of nuclear landscape to hyperheavy  nuclei is provided.

21.10.Dr, 21.10.Pc,  21.10.Ft, 21.60.Jz, 21.60.Ka

Extension of nuclear landscape to hyperheavy nuclei
    A. Gyawali
    December 30, 2023
===================================================







§ INTRODUCTION



 
   One of the main focuses of modern low-energy  physics is the
limits of the existence of finite nuclei. New generation of facilities
such as FRIB, FAIR, RIKEN, and GANIL  will explore such limits in 
neutron-rich nuclei. SHE-factory and  similar facilities will attempt to 
extend the limits of our knowledge on superheavy nuclei. However, 
already now it is clear that there are significant restrictions on what could be 
achieved by these new facilities: many neutron-rich medium mass, 
heavy and superheavy nuclei will be beyond their experimental reach 
<cit.>. In such a situation, theoretical predictions became the 
only tool to investigate such limits. Indeed, a significant progress has 
been achieved in understanding of the limits of nuclear landscape for 
the Z<120 nuclei (see Refs. <cit.>) and 
more or less consistent picture has been obtained using the combination 
of  different theoretical tools. In addition, systematic theoretical 
uncertainties <cit.> and statistical errors 
<cit.> in the predictions of the properties 
of neutron-rich nuclei and the positions of two-proton and two neutron-drip 
lines have been evaluated.

  However, the nuclear landscape  is not restricted to the Z<120
nuclei. Although there were some attempts to investigate higher 
Z nuclei <cit.>, 
these systematic studies were restricted to spherical symmetry.  Our recent 
study (Ref.  <cit.>) based on systematic axial Relativistic Hartree-Bogoliubov
(RHB) calculations and triaxial RHB as well as triaxial relativistic mean field + 
BCS (RMF+BCS) calculations for a reasonable large
set of selected nuclei has invalidated many conclusions of these 
studies[The effects of axial and triaxial deformations have also been 
studied for a few hyperheavy nuclei in Refs. <cit.>
and Ref. <cit.>, respectively. Somewhat larger set of the nuclei was 
studied with triaxiality included in Ref. <cit.> but according
to Ref. <cit.> the deformation range employed in this work is 
not sufficient for Z ≥ 130 nuclei.].  In addition, it provided a new
vision on the properties of hyperheavy nuclei and on the extension 
of nuclear landscape to the Z>120 region. These results are briefly 
summarized below. The increase of proton number beyond Z=120 leads 
to the dominance of highly deformed and superdeformed oblate ground 
states.  However, these states with ellipsoidal-like shapes become unstable 
with respect of fission in  the Z∼ 130 region (see also Ref. <cit.> 
for the results obtained for fission barriers in non-relativistic theories). This 
triggers the transition to the states with toroidal shapes; the lowest in energy 
solutions in the Z=140-180 nuclei have such shapes in axial RHB calculations. 
It was illustrated that some of such states could be stable against fission. In 
addition,  some regions of stability of spherical hyperheavy nuclei have been 
predicted for the first time in Ref. <cit.>.  Although these states are 
highly excited with respect of the lowest in energy states with toroidal shapes
(as obtained in axial RHB  calculations),  they will become the ground states 
if toroidal states are not stable with respect of multifragmentation (which 
according to present understanding (see Ref.  <cit.>) is quite 
likely scenario).

   Note that only in hyperheavy nuclei the states with toroidal shapes
could become the lowest in energy. The toroidal shapes in atomic nuclei 
have been investigated in a number  of the papers 
<cit.>.
However, in absolute majority of the cases such shapes correspond to 
highly excited states either at spin zero <cit.> or at extreme
values of angular momentum <cit.>.  In the former 
case, such states are unstable against returning to the shape of sphere-like 
geometry (Ref. <cit.>).  In the latter case, calculated angular momenta 
at which toroidal shapes appear substantially exceed the values of angular 
momentum presently  achievable at the state-of-art experimental facilities 
<cit.>.









  




  
    The present manuscript aims at the extension of the investigations of 
the properties of hyperheavy nuclei and of nuclear landscape started 
in Ref. <cit.>.  The topics covered by this investigation are shortly
mentioned in the next paragraph and discussed in details in the sections
below.

   The manuscript is organized as follows. The details of theoretical
calculations are discussed in Sec. <ref>. Section 
<ref> is devoted to the analysis of the effects of the 
truncation of the basis on the results of calculations. Density profiles,
charge radii and neutron skins of spherical nuclei located in the
centers of the islands of stability and their dependence on the functional
are considered in Sec. <ref>.  Sec. <ref> 
discusses the shell closures in the islands of stability of spherical 
hyperheavy nuclei.  The stability of spherical nuclei in these islands
with respect of octupole and triaxial distortions is investigated in
Sec. <ref>.  Sec. <ref> analyses 
the impact of octupole deformation on stability of prolate superdeformed 
minima. Systematic analysis of the results of the calculations for the 
Z=138 isotopic chain is performed in Sec. <ref>.
The stability of toroidal shapes in selected nuclei and the evolution
of such shapes along their fission path are considered in Sec. 
<ref>. Sec. <ref> is devoted to the 
analysis of the impact of triaxial deformation on the fission barriers of 
neutron-rich superheavy nuclei. The extension of nuclear landscape to 
hyperheavy nuclei is discussed in Sec. <ref>. 
Finally, Sec. <ref> summarizes the results of our work.





§ THE DETAILS OF THE THEORETICAL CALCULATIONS




 The investigations of the properties of super- and hyperheavy nuclei have been 
performed in different theoretical frameworks. Systematic investigation of hyperheavy
nuclei across the nuclear landscape between two-proton and two-neutron drip lines 
is performed within the axial reflection symmetric relativistic Hartree-Bogoliubov  
(RHB) framework (see Ref. <cit.>).  The stability of prolate minima with 
β_2 ∼ 0.5 of superheavy and low-Z hyperheavy nuclei as well as of 
spherical minima of hyperheavy nuclei with respect of octupole deformation has 
been studied with reflection asymmetric RHB framework using OCT-RHB code of 
Ref. <cit.>.  Triaxial RHB (TRHB code) <cit.> and triaxial relativistic 
mean field + BCS  (TRMF+BCS code)  <cit.> frameworks have been employed 
for the study of fission barriers in superheavy nuclei and stability of hyperheavy 
nuclei with respect to triaxial distortions.  Note that the TRHB and TRMF+BCS codes do 
not include octupole deformation.   Considering very time-consuming nature of the 
calculations in the OCT-RHB, TRHB and TRMF+BCS codes, only restricted set of nuclei 
has been investigated in their frameworks.





 The absolute majority of the calculations has been performed with the DD-PC1 
covariant energy density functional <cit.>. This functional is considered 
to be the best relativistic functional today based on  systematic and global studies 
of different physical  observables related to the ground state properties and fission 
barriers <cit.>.  Other 
functionals such as DD-ME2 <cit.>, PC-PK1  <cit.> and NL3* 
<cit.>), representing other major classes of covariant density functional models 
<cit.>, are employed only for the study of some properties of spherical 
nuclei located in the centers of the islands of stability of hyperheavy nuclei (see Fig. 6 
in Ref. <cit.>). This is done for the assessment of systematic theoretical 
uncertainties in the predictions of their properties.  

   
  The constrained calculations in employed codes perform the variation of the
function


    E_RHB/RMF+BCS + ∑_λ,μ C_λ,μ (⟨Q̂ _λ,μ⟩
                                - q_λ,μ)^2.


Here (λ,μ) = (2,0), (λ,μ) = (2,0) and (3,0) as well as 
(λ,μ) = (2,0) and (2,2)
in the  RHB, OCT-RHB and TRHB (TRMF+BCS) calculations, respectively.  E_RHB 
and E_RMF+BCS are the total energies in the RHB and RMF+BCS calculations.
<Q_λ,μ>  stand for the expectation values of the respective multipole
moments which are defined as


    Q̂_20    =     2z^2 -x^2 -y^2, 
    Q̂_22    =     x^2 -y^2, 
    Q̂_30    =     z(2z^2 -3x^2 -3y^2).


C_λ,μ in Eq. (<ref>) are corresponding stiffness constants 
<cit.> and q_λ,μ are constrained values of respective moments. To 
provide the convergence to the exact value of the desired multipole moment, 
we use the method suggested in Ref. <cit.>.  Here the quantity q_λ,μ
is replaced by the parameter q_λ,μ^eff, which is automatically 
modified during the iteration in such a way that we obtain
⟨Q̂_λ,μ⟩ = q_λ,μ
for the converged solution. This method works well in our constrained
calculations. In the OCT-RHB code we also fix the (average) center-of-mass of 
the nucleus at the origin with the constraint


    ⟨Q̂_10⟩ =0


on the center-of-mass operator Q̂_10 to avoid a spurious motion
of the center of mass.
   
The deformation parameters β_2,  β_3  

and γ are extracted from respective multipole moments:


    Q_20   =   ∫ d^3r ρ(r⃗) (2z^2-x^2-y^2),
    
    Q_22   =   ∫ d^3r ρ(r⃗) (x^2-y^2), 
    
    Q_30   =   ∫ d^3r ρ(r⃗) z(2z^2 -3x^2 -3y^2),




via


    β_2    =   √(5/16π)4π/3 Z R_0^2√( Q_20^2 + 2Q_22^2)
    γ   =   arctan√(2)Q_22/Q_20
    β_3    =   √(7/16π)4π/3ZR_0^3 Q_30


where R_0=1.2 A^1/3. Note that Q_22=0 and γ=0 in
axially symmetric RHB calculations.




 

  The β_2 and  γ values have a standard meaning of the deformations 
of the ellipsoid-like density distributions only for |β_2| ≲ 1.0 values. At higher 
β_2 values they should be treated as dimensionless and particle normalized measures 
of the Q_20 and Q_22 moments. This is because of the presence of toroidal 
shapes at large negative β_2 values and of necking degree 
of freedom at large positive β_2 values (see Fig. <ref> below).


 Note that physical observables are frequently shown as a function of the 
Q_20,  Q_30  and Q_22 moments. However, from our point 
of view such way of presentation has a disadvantage that the physical 
observables of different nuclei related to the shape of the density 
distributions (such as deformations) are difficult to compare because 
the Q_20, Q_30 and Q_22 moments depend on particle 
number(s).

  For each nucleus under study, the deformation energy curves in  
the -5.0 < β_2 < 3.0 range are calculated in the axial reflection 
symmetric RHB framework <cit.>; such large range 
is needed for a reliable definition of the β_2 value
of the lowest in energy minimum for axial symmetry (LEMAS).
This LEMAS becomes the ground state if the higher order deformations
(triaxial, octupole) do not lead to the instability of these minima.
 
To avoid the uncertainties connected with the definition of the size of the 
pairing window <cit.>, we use the separable form of the finite-range 
Gogny pairing interaction introduced in Ref.  <cit.>. Its matrix elements 
in r-space have the form


    V( r_1, r_2, r_1', r_2')    =   
    
    = -  G δ( R-   R'   )P(r) P(r') 1/2(1-P^σ)


with R=( r_1+ r_2)/2 and r= r_1- r_2
being the center-of-mass and relative coordinates. The form factor P(r) is of 
Gaussian shape


    P(r)=1/(4 π a^2)^3/2e^-r^2/4a^2


The parameters of this interaction have been derived by a mapping of the ^1S_0 
pairing gap of infinite nuclear matter to that of the Gogny force D1S. The resulting
parameters are: G=728 MeV fm^3 and a=0.644 fm <cit.>. This pairing 
provides a reasonable description of pairing properties in heaviest nuclei (actinides
and light superheavy nuclei) in which pairing properties can be extracted from 
experimental data <cit.>.







§ THE EFFECTS OF THE BASIS TRUNCATION






   Considering that the results published in Ref. <cit.> and presented in 
   this manuscript represent the first studies of hyperheavy nuclei and toroidal shapes
   in such nuclei within the covariant density functional theory, a special attention has 
   been paid to the investigation of the impact of the truncation of the basis on the 
   results of the calculations at different β_2 values. The truncation of the basis 
   is performed in such a way that all states 
   belonging to the major shells up to N_F fermionic shells for the Dirac spinors 
   are taken into account. Note also that the results of the calculations depend on
   the deformation β_0 of the oscillator basis. The detailed investigation of
   the impact of the deformation of the basis on the results of the calculations for the 
   ^466156 nucleus is presented for the N_F=20 and N_F=30 fermionic 
   shells in Figs. <ref>c and d, respectively. One can see that the 
   deformation of basis β_0 = 0.5 β_2 typically leads to the lowest in energy 
   solutions. Here β_2 stands for the deformation of the nucleus.
   This truncation scheme is also characterized by the fastest convergence. Similar situation
   has been observed also in other nuclei. Thus, this deformation of the basis has
   been used in all calculations of Ref. <cit.> and the present manuscript. 
   

  Figs. <ref>a and b show the dependence of the results of 
    calculations on the number of fermionic shells N_F for the ^208Pb
    and ^466156 nuclei. Note that the deformation of basis  
    β_0 = 0.5 β_2 is used for both of these nuclei. In ^208Pb,
    the N_F=20 basis provides very accurate description of binding
    energies in the physically interesting range of quadrupole deformations.
    Only at β_2 < -3.5 there is some difference between the results 
    obtained with N_F=20 and N_F=30. However, this is not physically
    significant range of the β_2 values since binding energies at these 
    values exceed binding energy of the ground state by at least 250 MeV.
    
         However, in hyperheavy nuclei the required size of the basis depends 
      both on the nucleus and deformation range of interest. The N_F=20 
      basis is sufficient for the description of deformation energy curves in the 
      region of -1.8 < β_2 < 1.8  (see Fig. <ref>b).
      The deformation ranges -3.0 <β_2 <-1.8 and 1.8 <β_2 <3.0 
      typically require N_F=24 (low-Z and  low-N hyperheavy nuclei) or 
      N_F=26 (high-Z and high-N  hyperheavy nuclei). Even more 
      deformed ground states with β_2 ∼ -4.0 are seen in 
      high-Z/high-N hyperheavy nuclei (see Fig. <ref>b
      for the ^466156 results and Fig. 1 in Ref. <cit.> for
      the ^426176 results); their description requires N_F=30. Thus, in 
      our studies the truncation of basis is made dependent on the nucleus 
      and typical profile of deformation energy curves or potential energy 
      surfaces.


       
       Fig. <ref>  illustrates the evolution of nuclear shapes 
 along the lowest in energy solution of hyperheavy ^466156 nucleus
 obtained in axial RHB calculations  with  N_F=30. Starting from 
 spherical  shape at β_2=0.0, the increase of prolate deformation 
 leads to the emergence of hyperdeformed shapes at β_2=1.0, which
  evolve into the shapes consisting of two fragments connected by neck at higher
  β_2 values.  The separation of the fragments and the size of the neck 
  increases/decreases with increasing β_2 values. Figs. <ref>b, 
  c and d clearly indicate  increased dependence of the results on the parameters 
  N_F and β_0 of the basis for these shapes.
      
      The evolution  of the shapes for negative β_2 values is shown in
   Figs. <ref>a-e. Highly deformed oblate shape exists at 
   β_2 = - 0.5 deformation which transforms into biconcave disk shape
   at β_2 = -1.0.  Further decrease of the β_2
   values leads to toroidal shapes. Note that with the increase of absolute value 
   of β_2 the radius of the toroid increases and the tube radius decreases.
   Total energies and equilibrium deformations of toroidal shapes with 
   β_2 ∼ -1.45,  corresponding 
   to local minimum  seen in deformation energy curves  of Figs. <ref>b-d,
   are rather  well described with N_F=20 and show almost no dependence on the
   deformation of basis β_0. On the contrary, toroidal shapes with larger 
   (in absolute sense) β_2 values show substantial increase of the dependence 
   of total energies and equilibrium deformations on N_F and β_0.

       
   Thus, the present analysis clearly indicates that the N_F=20 basis is 
sufficient for a description of classical ellipsoidal shapes and some toroidal 
shapes with relatively low absolute values of β_2 even in hyperheavy 
Z<180 nuclei. On the contrary, significantly  larger basis is required for a 
description of more exotic shapes such as toroidal ones with large absolute 
β_2 values and two-fragment ones connected by neck.
      
  The equilibrium β_2 values of the Z=122-138 nuclei presented 
in Fig. 2 of Ref. <cit.> have been calculated with fermionic 
bases including up to N_F=26 fermionic shells. The calculations for 
higher Z values require further increase of the size of fermionic 
basis (up to N_F=30 in high-Z/high-N nuclei). Such calculations are 
extremely time-consuming even in axial RHB framework and thus have not 
been undertaken. On the other 
hand, the type of the LEMAS can be established in the calculations with 
N_F=20. This is because even with N_F=20 the toroidal shapes 
with β_2 < - 1.4 represent 
the lowest in energy solutions at axial symmetry in the Z=140-180
part of nuclear landscape (see Fig. <ref>). This
figure clearly shows that classical ellipsoidal shapes are not 
energetically favored in hyperheavy nuclei. However, because of 
the limited size of the basis these β_2 values have to be 
considered as lower limits (in absolute sense). As illustrated in
Fig. <ref>, further increase of the size of fermionic basis 
will lead to the increase (in absolute sense) of the β_2 values 
of LEMAS and to more energetically favored status of toroidal shapes 
as compared with ellipsoidal ones.










§ DENSITY PROFILES, CHARGE RADII AND NEUTRON SKINS OF 
SPHERICAL HYPERHEAVY NUCLEI



  
  Three regions of spherical hyperheavy nuclei centered around 
(Z∼ 138, N∼ 230),  (Z∼ 156, N∼ 310) and 
(Z∼ 174, N∼ 410), which are expected to be reasonably 
stable against spontaneous fission and α-decay, have been 
predicted in Ref. <cit.>. The largest region with the 
highest fission barriers is centered at Z∼ 156, N∼ 310; 
other two regions are smaller with smaller fission barriers (see 
Fig. 6 in Ref. <cit.>). The CEDFs DD-PC1 and DD-ME2 predict
larger regions of stability and substantially higher fission 
barriers (reaching 10 MeV in some nuclei) as compared with the 
NL3* and PC-PK1 functionals.  The large fission barriers obtained 
in the density-dependent functionals will lead to substantial 
stability of spherical hyperheavy nuclei against spontaneous 
fission. This stability is significantly lower for the NL3* and 
PC-PK1 functionals.





   The nuclear matter properties and the density dependence 
are substantially better defined for density-dependent (DD*) 
functionals as compared with non-linear NL3* and PC-PK1 ones 
<cit.>. As a consequence, in general, they are expected 
to perform better for large extrapolations from known regions. 
In this context, it is also important to look on other features 
which may be critical in the discrimination of the predictions
of different functionals. Thus, this section is dedicated
to the analysis of charge radii, neutron skins and density
distributions of the nuclei located in the centers of this
potential islands of stability of hyperheavy nuclei. These are 
^368138, ^466156 and ^584174 nuclei. As a benchmark, 
we are using the ^208Pb and ^292120 nuclei. The properties
of latter nucleus were studied in details in Ref. <cit.>.





















  Charge radii r_ch and neutron skins r_skin of these 
nuclei are presented in Table <ref> and related
density distributions are shown in Fig. <ref>.
The predictions of different functionals for charge radii
of ^208Pb differ by less than 0.01 fm (see Table 
<ref>); this is also seen in proton density 
distributions (see Fig. <ref>a). On the
contrary, the spread in the predictions of neutron skin is
significant reaching 0.095 fm. Density dependent (DD) functionals
predict the lowest values for the neutron skin, while the NL3* 
functional predicts the highest value and the PC-PK1 results lie 
somewhere in between of the DD and NL3* ones. These features are 
reflected also in the neutron density distributions; they extend 
to higher radii in surface area and have lower density in the 
central region in the NL3* and PC-PK1 functionals as compared 
with the DD ones (see Fig. <ref>a). These differences 
between the functionals are realized when the neutron matter is 
moved from the surface region to the central and middle parts of 
the nucleus.  Comparable features are also seen in the ^292120 
nucleus (see Table <ref> and  Fig. <ref>b).

  However, these differences between the predictions of the
functionals become enhanced on going to the central nuclei
of the regions of potential stability of spherical hyperheavy 
nuclei. The spreads in the predictions of charge radii increase from
0.014 fm for the ^292120 nucleus to 0.052, 0.094 and 0.16 fm for 
the ^368138, ^466156 and ^584174 nuclei, respectively. 
The largest charge radius is always produced by the PC-PK1 
functional, while the smallest either by DD-PC1 or by DD-ME2. 
Note that latter two functionals give comparable results. 
 
   The spreads in the predictions on going from ^292120 nucleus
to higher Z/higher N nuclei increase also for neutron skins (see 
Table <ref>). For example, the difference in neutron skin 
of the ^584174 nucleus calculated with DD-PC1 and NL3* reaches 
0.176 fm. Similar to ^208Pb and ^292120 nuclei, density dependent 
(DD) functionals predict the lowest values for the neutron skin, while 
the NL3* functional predicts the highest value and the PC-PK1 results 
lie somewhere in between of the DD and NL3* ones.


  These results clearly indicate that the accuracy of the 
reproduction of charge radii and neutron skins by the CEDFs 
could be an important criteria in favoring or disfavoring
the predictions of one or another functional for the islands
of stability of spherical hyperheavy nuclei. Among considered
functionals, the DD-ME2 and DD-PC1 functionals provide the best
global description of charge radii (see Sect. X in Ref. 
<cit.>). However, the situation with neutron skins 
is more complex. Even for ^208Pb there is a significant controversy 
in the adopted experimental values of neutron skins (see discussion
in Sect. X  of Ref. <cit.> and in Ref. <cit.>). 
For example, the experiments based on hadronic probes provide neutron 
skin in ^208Pb around 0.2 fm or slightly smaller. However, these 
experimental data are extracted in model-dependent ways. Alternatively, 
a measurement using an electroweak probe has been carried 
out in parity violating electron scattering on nuclei (PREX) and it 
brings r_skin=0.33± 0.17 <cit.>. A central value of 0.33 
fm is particularly intriguing because it is around 0.13 fm higher 
than central values obtained in other experiments. Note that non-linear 
CEDFs typically give r_skin∼ 0.3 fm (see Table <ref>).
The electroweak probe has the advantage over experiments using hadronic 
probes that it allows a nearly model-independent extraction of the 
neutron radius that is independent of most strong interaction 
uncertainties <cit.>. Thus, the results obtained in future 
PREX-2 experiment <cit.> would be quite useful in
helping to discriminate the predictions.













































  General features of proton and neutron density 
distributions of hyperheavy nuclei seen in Fig. <ref> 
do not depend on employed functional. Apart of neutron density
in ^368138 nucleus (which is almost the same in the center of nucleus
and at its surface), both types of  densities are characterized by the 
density depression in the central part of the nucleus. 
Here we use the ratio ω = ρ^cent/ρ^surf of the 
density at the center ρ^cent to the maximum
density at the surface ρ^surf averaged over the set of
employed functionals to characterize this depression.
The central density depressions in neutron subsystems of the 
^466156 and ^584174 nuclei are rather modest with 
ω_ν= 0.814 and 0.86, respectively. Thus, neutron
densities are close to flat density distributions and could not be 
characterized as specific for semi-bubble 
nuclei (in the language of Ref. <cit.>).  However,  central 
depressions are significantly more pronounced in proton subsystems 
of hyperheavy nuclei with ω_π=0.753, 0.651 and 0.534 for 
the ^368138, ^466156 and ^584174 nuclei, respectively. 
Thus, they are close to the ones expected for semi-bubble nuclei 
(see Ref. <cit.>). Note that in a given nucleus the proton  
density is roughly half of the neutron one.

















   It is interesting that central depression is more pronounced in the 
^292120 superheavy nucleus (with ω_π∼ω_ν∼ 0.65)
as compared with higher-Z hyperheavy nuclei (see Fig. <ref>); 
the only exception is proton subsystem of the  ^584174 nucleus. The detailed 
analysis of Ref.  <cit.> strongly suggests that such central depression 
in the density distributions of the ^292120 nucleus is mostly due to underlying 
shell structure: it emerges due to the occupation of specific high-/low-j 
orbitals which place matter mostly in the surface/central region of the nucleus.  Indeed, the 
occupation of the neutron 3d_5/2, 3d_3/2 and 4s_1/2 orbitals on going
from N=172 to N=184 (the ^304120 nucleus) destroys this central
depression (see Fig. 2 in Ref. <cit.>). However, on going to hyperheavy
nuclei the role of strong Coulomb force (which pushes the matter to surface
region) in creation of central depression in density distributions is expected to 
become dominant <cit.>). The fact that the densities of the nuclei are 
similar within the regions of potentially stable spherical hyperheavy nuclei 
suggests reduced role of shell effects.











§ SHELL CLOSURES IN THE ISLANDS OF STABILITY
         OF SPHERICAL HYPERHEAVY NUCLEI




  The stability of spherical hyperheavy nuclei is  defined by underlying shell 
structure. However, in general the impact of shell gaps decreases with the 
increase of particle numbers (see discussion in Sect. III of Ref. <cit.>.
Fig. <ref> shows proton and neutron single-particle spectra 
of the ^466156  nucleus at spherical shape. Proton Z=154 and neutron 
N=308 shell gaps with the sizes of approximately 2 MeV are clearly visible 
for all employed functionals in this figure. Based on these spectra it is tempting to 
call the ^462154 nucleus as a doubly magic hyperheavy nucleus. However, 
the analysis of two-particle separation energies reveals more complicated situation.
The N=308 shell closure reveals itself via substantial drop of two-neutron 
separation energies at this particle number which exists for all proton numbers 
(see Fig.  <ref>a). This drop is also visible in calculated
α-decay half-lives (see Figs. 4 and 5 in supplemental material to
Ref. <cit.>). However, the impact of the Z=154 proton
shell closure on two-proton separation energies is substantially
smaller (see Fig. <ref>b) and it almost does
not exist for the N=308 nuclei.

  Although we have not performed detailed analysis of the separation
energies for the DD-ME2, PCPK1 and NL3* functionals, selected results
for α-decay half-lives presented in Fig. 6 of supplemental material
to Ref. <cit.> allow to perform the comparison with the DD-PC1 
ones. Similar to DD-PC1 the N=308 shell gap is expected to be seen in
two-neutron separation energies also for the DD-ME2 functional but its 
impact is somewhat smaller as compared with DD-PC1. On the contrary, 
the NL3* and PC-PK1 results presented for the Z=156 and Z=160 isotopic 
chains in Fig. 6 of supplemental material to Ref. <cit.> do not reveal 
the impact of the N=308 shell gap.  These observations suggest that 
in the center of the Z∼ 156, N∼ 310 island of stability of hyperheavy
nuclei total shell correction energies at spherical shape are more negative 
for the DD-PC1 and  DD-ME2 CEDFs as compared with the NL3* and PC-PK1
ones. This explains why fission barriers (and thus the size of the
islands of stability) [see Fig. 6 in Ref. <cit.>] are larger for the DD-PC1 
and  DD-ME2 CEDFs as compared with the NL3* and PC-PK1 ones.

 Figs.  <ref>  and fig.  <ref>  show single-particle 
 spectra for the ^366138 and ^580174 nuclei located in the vicinity of central 
 nuclei of two regions of spherical
 hyperheavy nuclei. Although the N=228 and N=406 neutron gaps with the size
 of approximately 2 MeV are seen in these figures, there are no substantial proton
 gaps at respective particle numbers. Similar to Fig. <ref>, these 
 neutron gaps are seen in two-neutron separation energies, but two-proton separation 
 energies are quite smooth as a function of proton number and do not reveal proton
 gaps. 
 
  The features observed for proton subsystems of the nuclei under discussion together 
with clear localization of the islands of stability of spherical hyperheavy nuclei in the 
(Z,N) plane strongly suggests that the shell  effects at deformed shapes leading to 
negative shell correction energies at some deformation and thus to fission barriers  play
also an important role in the stabilization of spherical hyperheavy nuclei in discussed 
regions.

   Some of discussed shell gaps appear as the gaps between the
members of the spin-orbit doublets. These are the N=308 shell gap 
between the 2i_13/2 and 2i_11/2 orbitals in the ^466156 
nucleus (Fig. <ref>b), the N=228 shell gap between the 
2h_11/2 and 2h_9/2 orbitals in the ^366138 nucleus 
(Fig. <ref>b) and the N=406 shell gap between 
the 2j_15/2 and 2j_13/2 orbitals in the ^580174 nucleus 
(Fig. <ref>b). The energy splitting between the spin-orbit 
partner orbitals depends on the profile of the density distribution in the surface 
region (see discussion in Ref. <cit.>). Indeed, the above mentioned 
gaps are similar in the DD-PC1/DD-ME2 and NL3*/PCPK1 pairs of the functionals 
reflecting the similarities and differences in their density distributions (see 
discussion in Sec. <ref>). Note that the gaps obtained in the  
DD-PC1/DD-ME2 functionals are larger than those seen in the NL3*/PC-PK1 
ones.









§ THE STABILITY OF SPHERICAL HYPERHEAVY NUCLEI WITH RESPECT OF
         OCTUPOLE AND TRIAXIAL DISTORTIONS




  Fig. 6 in Ref. <cit.> summarizes the heights of the fission 
barriers for the nuclei forming the islands of stability of spherical 
hyperheavy nuclei.  They represent the lowest in energy barriers 
amongst those obtained on oblate and prolate sides of spherical 
minimum in axial RHB calculations. However, one should investigate 
the stability of spherical minimum with respect of octupole and triaxial 
distortions to have a full understanding of the situation.  Such studies 
are very time-consuming and thus we present their results obtained with 
DD-PC1 CEDF only for the nuclei located in the centers of the islands 
of stability of spherical hyperheavy nuclei.

 Fig. <ref> shows the potential energy surfaces of 
 these nuclei in the (β_2, β_3) plane. One can see that in the 
 ^368138 and ^584174 nuclei, the barriers on the oblate and prolate 
 sides have the lowest values at β_3=0. The same is true also for the 
 barrier on the prolate side of spherical minimum in the ^466156 nucleus.
 However, potential energy surface is soft in octupole deformation in 
 the region of the barrier located on the oblate side of spherical 
 minimum in this nucleus. Thus, the saddle of the fission barrier
 is shifted from β_2=-0.196, β_3=0.0 (as obtained in axial
 RHB calculations without octupole deformation) to  
 β_2=-0.198, β_3=0.091  when octupole deformation is included
 in the calculations. This also leads to the decrease of the height of the
 barrier on the oblate side from 10.81 MeV down to  9.83 MeV. 
  However, this decrease has very little impact on the total stability 
 of this nucleus with respect of octupole deformation since resulting
 value of 9.83 MeV is only 120 keV lower than the height 
 E_B=9.95 MeV  of the barrier on the prolate side.  Note that 
 the barrier on the prolate side is the lowest one in the axial RHB
 calculations without octupole deformation. Thus, one can conclude 
 that above discussed nuclei are relatively stable with respect of 
 octupole distortions.  

  It is also necessary to mention that potential energy surfaces shown in Fig.  
<ref> do not suggest that cluster radioactivity from spherical 
^368138 and ^466156 hyperheavy nuclei plays an important role. Super-asymmetric 
fission channel leading to a cluster radioactivity shows itself as a narrow fission path 
in the (β_2,β_3) plane separate from main fission path (see example in Fig. 1 
of Ref. <cit.>). No such path is visible in Figs. <ref>a 
and b. However, because of convergence problems one cannot define whether
cluster radioactivity is important in the ^584174 nucleus.
 

 Fig. <ref> shows the potential energy surfaces obtained
in the TRHB calculations. In the ^584174 nucleus, the axial saddles are located at 
β_2 = 0.17, γ = 60^∘ (the barrier on the oblate side of spherical
minimum) and at  β_2 = 0.197, γ = 0^∘ (the barrier on the prolate 
side of spherical minimum). Their heights are 6.389 and 7.709 MeV, respectively.  The 
potential energy surface is somewhat soft in triaxial deformation so that the saddle of 
fission barrier is shifted to β_2=0.16, γ =36^∘ when the triaxiality is included 
in the calculations. However, the impact of triaxiality on the height of fission barrier is 
rather modest leading to its decrease (as compared with the lowest axial barrier located 
on oblate side) by only 0.26 MeV. Similar features are also seen in the ^368138 and 
^466156 nuclei. The reduction in the height of fission barrier is 0.71 and 0.58 MeV 
in these nuclei.  However, as compared with the heights of 5.88 and 10.81 MeV of the 
lowest axial barrier  (located at the oblate side  of spherical minumum), these are relatively 
modest reductions which do not decrease the stability of nuclei in a substantial way. 
 
  These features could be understood in the following way. The topology 
of potential  energy surfaces of the nuclei under study are similar to those 
of volcanos.  The central area around spherical minimum is similar to caldera, the rim
of which is represented by the fission barrier.  The area beyond the rim
(fission barrier) is fast down-sloping as a function of quadrupole deformation
β_2. The saddles of axial fission barriers (on oblate and prolate sides
of spherical minimum) are located at modest quadrupole deformation of 
β_2 ∼ 0.2. As a result, the distance between these two saddles
in the (β_2, γ) plane is relatively small, so that large
changes in binding  energy due to triaxiality for nearly constant β_2
values could not develop.  As a
consequence, the lowest fission barrier around spherical minimum obtained in 
axial RHB calculations is a  good approximation to the barrier obtained
in the TRHB calculations. The TRHB results discussed here clearly indicate
that spherical minima of the nuclei under study are relatively stable with respect of 
triaxial distortions.





§ PAIRING INTERACTION IN HYPERHEAVY NUCLEI




The magnitude of pairing interaction could be accessed via two calculated quantities: average pairing gap and pairing energy  <cit.>.  There are several definitions of average pairing gaps in literature (see discussion in Section IV  of 
Ref. <cit.>). The pairing gap


    Δ_ uv=∑_k u_kv_kΔ_k/∑_k u_kv_k ,


which is related to the average of the state dependent gaps over the pairing 
tensor, is used in the present study. The analysis of Ref. <cit.>
showed that the Δ_ uv gap is a better measure of pairing
correlations as compared with other definitions for average pairing gap.

 The pairing energy in the RHB calculations  is defined via


    E_pairing = -1/2 Tr (Δκ)
   

where Δ and κ are pairing field and pairing tensor, respectively 
<cit.>.  Note that E_pairing mixes particle-particle and particle-hole
channels of the model. As a consequence, its absolute value is typically by an 
order of magnitude larger than the gain in binding due to pairing correlations
which is equal to the difference of binding energies obtained in the calculations 
with and without pairing correlations (see Refs. <cit.>).

  Fig. <ref> shows the evolution of neutron and proton pairing 
energies E_pairing and pairing gaps Δ_uv as a function of
β_2 for the lowest in energy solution in the ^466156 nucleus.
The oscillating behavior of these quantities as a function of β_2
is due to the change of the density of the single-particle states in the 
vicinity of the Fermi level with deformation (see detailed discussion in 
Ref. <cit.>). In regions of high/low level density it is easier/more 
difficult for the quasi-particles to spread  around the Fermi surface and 
therefore the size of the pairing correlations depends strongly on the 
level density. As a consequence,  the low/high values of the |E_pairing| and
Δ_uv corresponds to low/high density of the single-particle states
in the vicinity of the Fermi level 
(see Ref. <cit.>). Indeed, in the nucleus under consideration 
there is a substantial reduction of neutron |E_pairing| and Δ_uv 
values near spherical shape which is attributable to the presence of large 
N=308 shell gap (see Fig. <ref>). Note also that there is a collapse 
of neutron pairing correlations at toroidal shapes in the β_2 range
from -4.2 up to -3.95. This range corresponds to the minimum of the 
deformation energy curve at toroidal shapes  and its vicinity (see Fig. 
<ref>d). The local minima in potential energy surfaces correspond 
to low density of the single-particle states in the vicinity of the Fermi level <cit.>; at above quoted deformations this density is low enough to trigger the collapse of pairing 
correlations in the RHB calculations. However, this collapse of pairing in neutron subsystem is not critical since the treatment of pairing by more sophisticated 
methods, which preserve exact particle number, would only lead to moderate 
additional binding due to pairing at quoted above deformations. 
 
  Particle number dependencies of proton and neutron pairing energies at 
spherical shape of the nuclei forming the (Z∼ 156, N∼ 310) island of 
stability of spherical  hyperheavy nuclei are shown in Fig. <ref>.
Apart of the N=308 isotopes, the absolute values of neutron pairing energies 
|E_pairing| are larger than 4 MeV. In the N=308 isotones, they are
smaller than 4 MeV and there is a collapse of neutron pairing due to large 
size of the N=308 shell gap in the nuclei with Z=158-174.
The treatment of pairing by the methods which include exact particle 
number projection (see, for example, Ref. <cit.>) will restore the pairing 
correlations in these nuclei.  However, in no way it will affect the conclusions
of the present work and of Ref. <cit.>. In reality, it is expected that
such methods most likely will increase the fission barriers around spherical
minimum because of additional binding due to pairing at spherical shape. 
Proton pairing energies are displayed in Fig. <ref>b; only 
in the Z=154, N=292-298 nuclei their absolute values are below 4 MeV.
However,  there is no proton pairing collapse in any of the nuclei shown in 
this figure. In other nuclei, proton pairing energies are quite large with 
|E_pairing| > 6.0 MeV and for the Z≥ 168 nuclei  the |E_pairing| 
values exceed 20 MeV.

  Figure <ref> shows the evolution of proton and neutron
pairing energies as a function of the β_2 and γ deformations
in the ^466156 nucleus. Similar to the calculations in axially symmetric 
case, there are substantial oscillations of the pairing energies as a function of 
deformation parameters which are due to underlying changes in the 
single-particle level density in the vicinity of the Fermi level. However, the 
topology of  these oscillations is more complex because of the presence of 
two deformation parameters. Note  that no pairing collapse is observed at 
any deformation point covered by these TRHB calculations.
    




§ THE IMPACT OF OCTUPOLE DEFORMATION ON THE STABILITY
OF PROLATE SUPERDEFORMED MINIMA




 It is well known fact that outer fission barriers exist in superheavy 
nuclei for many CEDFs when the calculations are restricted to axial 
reflection symmetric shapes <cit.>.  This is illustrated in 
Fig. <ref> which shows the competition in energy of 
two local minima, namely, prolate superdeformed minimum with β_2 ∼ 0.6 
and oblate one with β_2 ∼ -0.5. Although the prolate superdeformed minima 
are the lowest in energy, their stability depends on the properties of 
outer fission barriers. It turns out that in absolute majority of the Z>120 super- and 
hyperheavy nuclei such local minima either do not exist or the heights of outer fission 
barrier are less than 2 MeV. The latter is not sufficient for the 
stabilization of prolate superdeformed  minimum (see discussion in Ref. <cit.>). 
Note also that such local prolate minima do not exist in the results of 
axial reflection symmetric RHB calculations for hyperheavy nuclei with 
Z > 140 (see examples shown in Fig. 1 of Ref. <cit.>).


 For the cases in which the heights of outer fission barriers are higher 
than 2 MeV in axial reflection symmetric RHB calculations (as those 
shown in Fig. <ref>), we have performed the calculations 
in axial octupole deformed code RHB-OCT developed in Ref. <cit.>. 
The results of such calculations are illustrated in Fig. <ref>. 
One can see that the inclusion of octupole deformation leads to the 
instability of the prolate superdeformed β_2 ∼ 0.6 minima with respect of 
octupole deformation. In the ^316,320124 nuclei this instability will 
lead to asymmetric fission, while the competition of symmetric and 
asymmetric fissions is possible in the ^296,300122 nuclei. 
The systematic calculations for the cases in which the heights of 
outer fission barriers are higher than 2 MeV in axial reflection 
symmetric calculations clearly show that the inclusion of octupole 
deformation leads either to complete disappearance of outer fission
barrier (as seen in the cases of  ^296122 and ^316,320124
nuclei in Fig. <ref>) or to a significant reduction 
of the heights of outer fission barriers to the values which are 
substantially lower than 2 MeV (as seen for the ^300122 nucleus 
in Fig. <ref>).  Thus, the prolate superdeformed minima are 
expected to be unstable in all Z>120 nuclei. Note that outer 
fission barriers in super- and hyperheavy nuclei could also be 
affected by triaxiality (see Ref. <cit.>). However, the 
present analysis shows that prolate minima in the Z>120 nuclei 
are already unstable in axially symmetric calculations. This is 
a reason why nuclear landscape in the Z=122-130 nuclei is 
dominated by the oblate ground states (see Fig. 3 in Ref. 
<cit.>).




§ SYSTEMATIC ANALYSIS OF THE RESULTS OF CALCULATIONS FOR THE
Z=138 ISOTOPIC CHAIN.




  In order to illustrate the variation of the properties of 
the nuclei with neutron number, a detailed analysis of the 
results of the calculations for the Z=138 isotopic chain is 
presented in this section.

 Deformation energy curves of even-even Z=138 hyperheavy nuclei 
obtained in axial RHB calculations are shown in Fig. <ref>. 
The minimum of deformation energy curve is located at β_2 ∼ -4.0 
for proton-rich nuclei with N=186-190. This minimum appear at 
β_2 ∼ -2.5 for the N=194-230 nuclei. Further increase of 
neutron number leads to the growth (in absolute sense) of the 
β_2 values: the nuclei with N=246-262 have minima at 
β_2 ∼ -3.0. The minima of the deformation energy curves are 
located at β_2 ∼ -3.5 for neutron-rich nuclei with N=266-326.
As discussed in details in Ref. <cit.>, the nuclei have
toroidal shapes in such minima. The competing local minima with
β_2 ∼ -0.8 (which corresponds to oblate ellipsoidal shape, 
see Ref. <cit.>) are located at high excitation energies 
of 10-30 MeV. This excitation energy depends on the nucleus but 
its maximum is reached at N=286-290 and minimum at N∼ 206.
Note the complex pattern of deformation energy curves at low 
deformation with a number of coexistent local minima. 
 
   Similar to few nuclei discussed in Ref. <cit.>, these minima 
in deformation energy curves with toroidal shapes are potentially 
unstable with respect of the transition to prolate shape via 
γ-plane and subsequent fission since prolate 
shapes with corresponding quadrupole deformations are located 
at lower energies (compare dashed lines [which represent mirror 
reflection of the positive β_2 part of deformation energy curve 
onto  negative β_2 values] with solid ones in Figs. 1c and d of 
Ref. <cit.>).  However, as discussed in Ref. <cit.>
and in Sec. <ref> below, some of these minima 
could be stable.

  Proton and neutron chemical potentials for the solutions, 
displayed in Fig. <ref>, are shown in Fig. <ref>. 
They behave differently as a function of β_2. Neutron chemical 
potential on average is almost flat as a function of β_2. The 
magnitude of the fluctuations of the  λ_n(β_2) values 
with respect of average values decreases on approaching two-neutron 
drip line. On the contrary, there is a pronounced slope in the 
λ_p(β_2) values: they on average decrease with 
increasing absolute value of β_2. Note that this slope 
is especially pronounced in proton-rich nuclei. As a consequence, 
in the nuclei with N=186-210 there is the range of the 
β_2 values in which the proton chemical potential is 
positive (see top row of Fig. <ref>) and the shape 
of the nucleus in respective local minimum is oblate ellipsoidal. 
Even-even nuclei with λ_p>0 are typically expected 
to be unstable with respect of proton emission[
The discussion of the two-proton drip line in terms of 
proton chemical potential has its own meaning. Strictly
speaking the two-proton drip line is reached when two-proton 
separation energy S_2n becomes negative. Alternatively (but 
less strictly) the position of the two-proton drip line is 
defined via the proton chemical potential λ_p=dE/dZ
as a point (nucleus) of the transition from negative λ_p 
(“bound” nuclei) to positive λ_p (“unbound” nuclei)
values. This definition depends on the employed pairing model. 
In addition, it presents a linear approximation in a Taylor 
expansion and, therefore, it ignores nonlinear effects like shape 
changes on going from the (Z-2,N) to the (Z,N) nucleus and 
their contribution to S_2p. However, even in the case of 
two-neutron drip line (which extremely sensitively depends on 
the fine details of the calculations) this definition leads in 
approximately two-thirds of the cases to the same two-neutron 
drip line as obtained in the definition of the two-neutron 
drip line via the separation energies <cit.>. In the 
remaining one-third of the cases, it leads to a two-neutron
drip line which is two neutrons short of the two-neutron drip
line defined via the separation energies; the nucleus which
is unbound (as defined via the chemical potential) has in
most of the cases a low positive value of λ_n ∼  0.05 MeV.
These results were obtained in Ref. <cit.> from the 
calculations of Refs. <cit.> by analyzing 
the two-neutron drip line positions of 60 isotopic chains for 
four different CEDF’s.].

 On the contrary, for the N=186-210 nuclei proton chemical 
potential is negative below β_2 ∼ -1.5 and it becomes 
more negative with the increase of absolute value of β_2.
As a consequence, toroidal shapes in these nuclei are expected 
to be stable with respect of particle emission. This feature 
in the behavior of the proton chemical potential as a function
of β_2 leading to instability/stability of ellipsoidal/toroidal 
shapes in the same nucleus with respect of particle emission 
is a source of unusual shift in the position of two proton-drip 
line towards more proton rich nuclei (as compared with general 
trend seen in the (Z,N) plane for the Z < 120 nuclei) which 
is clearly visible in Fig. <ref> below. 
Note that such shift is absent for two-neutron drip line most 
likely because of above mentioned flatness of neutron chemical 
potential as a function of β_2. 





§ TOROIDAL SHAPES: STABILITY AND SHAPE EVOLUTION
ALONG THE FISSION PATH.




    The investigations of Ref. <cit.> showed that some toroidal shapes 
 could be stable with respect of triaxial distortions. Fig. <ref> shows 
 potential energy surfaces  of the ^354134 and ^348138 nuclei around 
 minima of such configurations  located at β_2  ∼ 2.3, β_4 ∼ +1.5, 
 γ = 60^∘. The saddle points of the first fission barriers of these 
 configurations are located  at 4.4 and 8.54 MeV, respectively. However, physical 
 reasons for such stability of toroidal shapes have not been discussed in Ref. 
 <cit.>.
 
   To understand these reasons the evolution of toroidal shapes along the 
fission path of the configuration in the ^354134  nucleus (shown in Fig.
<ref>) is displayed in Fig. <ref>. The toroid and its 
tube are fully symmetric at the minimum (Fig. <ref>a-c).  The deviations 
from axial symmetry lead to the distortions which are already seen at 
γ=50^∘; the toroid is stretched out in the direction of the
axis of symmetry  and squeezed in perpendicular direction (Fig. <ref>f)).
However, the tubes of the toroid still remain cylindrical (Fig. <ref>d,e)).
Figs. <ref>g,h,i show the density distributions at the deformations
corresponding to the saddle point. One can see further increase of the asymmetry
of torus and the decrease of the area of toroid hole.   Thus, one can 
conclude that the barrier against fission emerges because these deviations 
from symmetrical shape of toroid cost the energy.

    Further decrease of the γ- and β_2 deformations leads to
increasing  distortion of the geometry of toroid ((Fig. <ref>l)) the
outer shape of which looks now similar to parallelogram and the shape of 
toroid  hole comes closer to square.  In addition, Figs. <ref>j,k  reveal 
visible deviations from cylindrical shape of the tube of toroid.  However, 
these changes reduce the total energy of the configuration as compared
with the one at the saddle point. 

 Above discussed changes in shapes and total energies along the fission 
path are the consequences of a delicate balance of liquid drip and shell correction 
energy contributions.






§ THE IMPACT OF TRIAXIAL DEFORMATION ON THE FISSION BARRIERS OF
              NEUTRON-RICH SUPERHEAVY NUCLEI





    Although oblate minima of high-Z (Z>120) superheavy and low-Z 
hyperheavy nuclei are relatively stable with respect of axial reflection symmetric and 
asymmetric deformations (see Ref. <cit.> and Sec. <ref> 
in the present manuscript), that is not necessary the case with respect of triaxial 
deformation. The impact of triaxiality on the fission paths and the differences in 
the stability of super- and hyperheavy elements is illustrated in Fig. 
<ref> on the example of superheavy ^268Sg and 
^332Ds nuclei and hyperheavy ^360130 and ^432134 nuclei.
 
  Both in super- and hyperheavy nuclei the potential energy surfaces (PES) represent
the canyon in which some local valleys and hills are located. However, there are
two principal differences between super- and hyperheavy nuclei. The canyon  
is quite narrow in superheavy nuclei which prevents the formation of local minima
at large oblate deformation and limits the role of triaxial deformation. One can see
that normal deformed minima are prolate in superheavy ^268Sg and 
^332Ds nuclei and fission paths from these minima is located not far away
from the γ=0^∘ axis. In addition, the bottoms of the canyons  in PES are
on average flat.

   On the contrary, in hyperheavy nuclei the walls of the canyon with very rapid 
raise of energy with deformation are located at larger separation (so only right wall 
is seen in the bottom  panels of Fig. <ref>) as compared with
superheavy nuclei and the mountain centered around β_2 ∼ 0 is formed 
in this canyon.  The slope of the mountain in the direction of the β_2-deformation 
at  γ=0^∘ is very high. This indicates higher instability of hyperheavy 
nuclei against fission as compared with superheavy ones.  The larger separation of 
the canyon walls leads to an increased role of triaxiality in hyperheavy nuclei: local 
minima are formed either at oblate superdeformation (see example of the ^360130 
nucleus in Fig. <ref>) or at very large γ-deformation (see 
example of the ^432134 nucleus in Fig.  <ref>).   In addition,
the fission paths from these minima proceed at larger γ-deformations
as compared with superheavy nuclei.  Not only the fission through the γ-plane 
gets more energetically favored, but also  the fission path through γ-plane 
becomes much shorter than the  one through the γ=0^∘ axis (see
also the discussion in Ref. <cit.>).
    
   The general conclusion is that the barriers along the fission paths 
emerging from the oblate minima located within the -1.0 < β_2 ≤ 0.0 range 
decrease with increasing proton number (see Ref. <cit.>). As a result (see
discussion in  Sect. <ref> below), the fission barriers for  oblate 
ellipsoidal shapes become consistently lower than 2 MeV  above some (Z, N) line 
in nuclear landscape so the nuclei in ellipsoidal shapes cease to exist for these
particle numbers.

  However, in order to delineate this borderline additional information on the impact
of triaxiality on the fission barrier heights of the superheavy Z=106-118 nuclei 
located between two-proton and two-neutron drip lines is needed.  So far, such 
information is available only for  actinides and superheavy nuclei with Z ≤ 120 
and N≤ 184  <cit.> .
These nuclei are either 
prolate or spherical in their ground states and thus the impact of triaxiality is limited: 
the lowering of inner fission barriers in actinides due to triaxiality is typically on the 
level of 1-3 MeV. This reduction is substantially smaller as compared with the one
typically seen in hyperheavy nuclei <cit.>.

  Unfortunately, even nowadays fully systematic triaxial  RHB calculations are 
extremely computationally demanding.  Thus, in order to get at least rough 
outline of the impact of low fission barriers on potential topology of nuclear 
landscape, we have performed triaxial RHB calculations for selected 
nuclei with Z=106, 110, 114 and 118 and with N=152 + Δ N, where 
Δ N = 0, 10, 20, ... and only the nuclei between two-proton and 
two-neutron drip lines are considered here. Note that in some nuclei
there are two fission barriers. If the outer fission barrier is lower than 2 MeV
in axial RHB calculations, it is ignored and the TRHB calculations are focused on 
the inner fission barrier.  If the outer fission barrier is higher than 2 MeV, then 
we first perform TRHB calculations in order to see whether triaxiality lowers outer
fission barrier below 2 MeV. If that is no a case, we carry out RHB-OCT calculations
in order to see whether octupole deformation lowers outer fission barrier below 2 
MeV. The results of such calculations, which provide information on the highest
fission barrier in the nuclei under consideration, are summarized in Table 
<ref>. Together with the results of the TRHB calculations presented in 
Refs. <cit.>  for superheavy Z=122 and hyperheavy 
Z=126, 130 and 134 nuclei they allow in an approximate way to delineate the 
impact of fission on the boundaries 
of nuclear landscape[The boundaries of nuclear landscape in heavy nuclei
with ellipsoidal shapes in the ground states
are defined by spontaneous fission and not by the particle emission as in lower 
Z nuclei (see Ref. <cit.> and the discussion in Sec. <ref>).
This fact has been ignored in many studies of nuclear landscape in 
neutron-rich actinides and superheavy nuclei (such as Refs.  <cit.>) 
since the boundaries of nuclear landscape were defined via two-neutron 
and two-proton separation energies. The reasons for such a choice are obvious: 
such calculations require only binding energies of the ground states which are 
relatively easy to compute. On the contrary, the calculations of fission barriers in 
triaxial  DFT codes are by approximately three orders of magnitude more 
numerically time-consuming. ] in Sec. <ref> below.

   Note that these are  first ever triaxial calculations for fission barriers in
neutron-rich (N>200) 
superheavy nuclei and as such they can be very useful for understanding
the role of the fission in the r-process calculations. The fission of heavy systems
affects this process via fission recycling <cit.>. However, so far 
the fission barriers for such nuclei were calculated only in theoretical frameworks 
restricted to axial symmetry <cit.>. 

   The results for inner fission barrier heights E^B_triax and the decrease
 of the height of inner fission barrier due to triaxiality Δ E^gain, presented
 in Table <ref>,  are summarized in Fig. <ref>.  The 
 E^B_triax values show oscillatory behavior as a function of neutron number N
 with maxima seen at N∼ 162 and N∼ 222 (at N=172 and 232 in the 
 Z=118  nuclei) and minima at N=192 and N∼ 252 (at N=202 and 272 in 
 the Z=118  nuclei). More erratic behavior is  seen for the Δ E^gain values; 
 note that the triaxiality has no effect on the heights of inner fission barriers for  
 Δ E^gain=0 MeV.  In some nuclei the triaxiality lowers inner fission barrier
 by more than 4 MeV. These features seen in E^B_triax  and  Δ E^gain 
 curves are the consequences of underlying shell structure and its evolution with proton 
 and neutron numbers (see discussion in Ref. <cit.>). A general trend of the 
 lowering of the height of inner fission barrier with the increase of proton number is 
 clearly seen in Table <ref> and Fig. <ref>a. 
 
   Table <ref> clearly indicates three regions of instability based on the fission barrier heights;
in these regions the height of the highest fission barrier is below 2 MeV.  These are two islands of 
instability centered around (Z∼ 114, N∼ 192) and (Z∼ 118, N∼ 202) 
which are shown in Fig. <ref>b. In addition, very neutron rich nuclei near and above 
N=252 in the Z=114 isotopic chain as well as near and above N=262 in the 
Z=118 isotopic chain are unstable with respect of fission. Such instability against fission 
is also seen in very neutron rich hyperheavy nuclei (see supplemental material to Ref. <cit.>). 
These factors together lead to a substantial reduction of the region of potentially stable elliposoidal-like 
nuclei in the N≥ 258 region (compare panels (a) and (b) in Fig. <ref> below).
 




§ EXTENSION OF NUCLEAR LANDSCAPE TO HYPERHEAVY NUCLEI




  One of important goals of the present manuscript is the 
extension of nuclear landscape to the limits of extreme Z values.
There are numerous studies of the limits of nuclear landscape at 
the neutron and proton drip lines for the Z<120 nuclei carried 
out in different theoretical frameworks (see Refs. 
<cit.> and,
in particular,  the compilation presented in  Sec. VIII of Ref. 
<cit.>). The studies of Refs. <cit.>
also define systematic theoretical uncertainties in the position of 
two-proton and two-neutron drip lines. On the contrary, 
nothing was known about the nuclear landscape in hyperheavy nuclei 
and its limits before our previous publication in Ref.  
<cit.>. The goal of this section is to present a comprehensive 
summary on the structure and limits of nuclear landscape with 
special emphasis on the region of hyperheavy nuclei.

  The results of Ref. <cit.> and the present study clearly 
show that critical distinction between the parts of nuclear chart  
are related  to the dominance of two different types of shapes: 
ellipsoidal-like and toroidal ones. Note that concave disk shapes, 
appearing at large oblate deformation, belong to ellipsoidal-like shapes. 
Fig. 1a shows the region of nuclear chart which is 
dominated by ellipsoidal-like shapes. Note that for Z<120 LEMAS 
obtained in reflection-symmetric RHB calculations  typically 
correspond to the ground states since only few nuclei in their 
ground states  are affected by γ-deformation (see Ref. 
<cit.> for the results obtained in microscopic+macroscopic method) 
and octupole deformation shows up in the ground states of the nuclei in 
few localized regions <cit.>. 

  The situation changes in the Z>120 nuclei which are typically  soft with 
respect of triaxial deformation up to the point that many ground 
states possess triaxial deformation (see Table 1 in supplemental
material of Ref. <cit.> and Table <ref> in 
the present manuscript). This softness also leads to a substantial 
reduction of the heights of the fission barriers in many nuclei. 
If the barrier height is less than 2 MeV, the nucleus is typically 
considered unstable against fission (see discussion in Ref. 
<cit.>). If we take this 
fact into account, the region of nuclear chart with ellipsoidal-like 
shapes will be considerably reduced at high-Z values; this is 
illustrated in Fig. 1b. Note that in some nuclei eliminated on 
transition from panel (a) to panel 
(b) the local minima (which are otherwise excited ones) with toroidal 
shapes could become the lowest in energy solutions if ellipsoidal-like 
shapes are unstable with respect of fission. Finally, two-proton 
and two-neutron drip lines for toroidal shapes are added in Fig. 1c.

  While there is the coexistence of ellipsoidal-like and toroidal 
shapes in the Z=120-140 part of nuclear chart (see Fig. 1 and 
Fig. 3 in Ref. <cit.>), with increasing proton number 
beyond Z=140, the LEMAS always have toroidal shapes (see 
discussion in Sec. <ref> and in Ref. 
<cit.>). The nuclear chart extended up to Z=180 displays 
the two-proton and two-neutron drip lines for toroidal nuclei 
outlining the potential  limits of nuclear landscape (Fig.  
<ref>).

The transition from ellipsoidal to toroidal shapes is driven by Coulomb
repulsion and has a lot of similarities to Coulomb frustration phenomenon
seen in nuclear pasta phase of neutron stars.  Fig. <ref> shows
the calculated Coulomb energies E_Coul as a function of the β_2
values. One can see that in all nuclei the largest Coulomb energy is 
calculated at spherical shape which is the most compact shape for a given 
volume. The deviation from sphericity decreases the 
Coulomb energy and for a given absolute value of β_2 this effect  
is especially pronounced for negative β_2 values. Moreover, 
the magnitude of E_Coul and its slope with deformation
for negative β_2 values drastically increases with increasing proton 
number (see Fig. <ref>).  This is also quantified in Fig. <ref> 
by the Δ E_Coul quantity which is the 
reduction of the Coulomb energy due to the transition from spherical
shape to typical toroidal one. The Δ E_Coul increases from 
346 MeV in ^208Pb to 721, 874 and 1126 MeV in the ^354134,
^466156 and ^426176 nuclei, respectively. This clearly tells
that toroidal shapes are more energetically favored by Coulomb interaction 
than spherical (or ellipsoidal-like) ones. However, only in hyperheavy nuclei 
the Coulomb energy becomes strong enough to trigger the transition to toroidal 
shapes in the lowest in energy solutions of axial RHB calculations (see Fig. 1 
in Ref. <cit.>).


    It is clear that nuclear landscape shown in Fig. 
<ref> is not complete because 
it does not take into account the potential instabilities 
of toroidal shapes with respect of different types of 
distortions. Although it was shown in Ref. <cit.> 
that some toroidal nuclei are potentially stable with 
respect of triaxial distortions, this was illustrated 
only for two nuclei. The underlying mechanism of
their stability is discussed in Sec. <ref>. 
The problem is that with increasing proton number the extreme 
sizes of fermionic basis (up to N_F=30 for nuclei at Z=156 
and beyond, see discussion in Sec. <ref>) 
are required for the description of toroidal nuclei. Neither 
triaxial nor axial reflection asymmetric calculations are 
possible nowadays for such sizes of basises. 

 The investigations of Ref. <cit.> suggest that 
toroidal shapes are expected to be unstable with respect of 
so-called sausage deformations which make a torus thicker in one 
section(s) and thinner in another section(s). They are 
expected to lead to multifragmentation[There are 
also some experimental indications of the role of 
multifragmentation in toroidal nuclei, but they are restricted 
to a single  ^86Kr + ^93Nb reaction <cit.>.].  However, these 
investigations are performed in the liquid drop model which neglects
potential stabilizing role of the shell effects. In 
addition, they do not extend beyond the region
of superheavy nuclei. To clarify the situation
the DFT studies of potential stability of toroidal shapes
with respect of triaxial distortions similar to those performed 
for ^354134 and ^348138 nuclei in Ref. <cit.> 
(see also Sec. <ref> in the present 
manuscript) are needed for higher Z values. Unfortunately, 
as mentioned above they are not possible nowadays because of 
extreme sizes of basises.

   If the toroidal shapes are unstable (and, from our point of 
view, the likelihood of this scenario is high in high-Z nuclei),
then the spherical shapes become the ground states in the
islands of potential stability of spherical hyperheavy nuclei
(see Ref. <cit.>). These islands centered around
(Z∼ 138, N∼ 230),  (Z∼ 156, N∼ 310) and 
(Z∼ 174, N∼ 410) are shown in 
Fig. <ref>.

    The analysis of Figs. <ref> and 
<ref> clearly indicates 
that the classical structure of the nuclear landscape in 
which spherical shell closures at different particle numbers
play a defining role disappears in the Z>120 nuclei.
This is because the ground states are either oblate or
toroidal in axial RHB calculations.

  The extrapolation to unknown regions is definitely associated
with theoretical uncertainties <cit.> which are especially large for
the position of two-neutron drip line 
<cit.> and fission 
barriers <cit.>. In the CDFT 
framework, systematic theoretical
uncertainties due to the form of the CEDF are substantially
larger than statistical errors <cit.>. So far, such
systematic uncertainties 
have been estimated only for fission barriers in the regions of 
potentially stable spherical hyperheavy nuclei and for the sizes
of these regions (see Fig. 6 in Ref. <cit.>). Their more 
global  evaluation is extremely time consuming and at this stage,
when we try to understand the general features of hyperheavy
nuclei, is not necessary. These uncertainties will definitely
affect the stability of ellipsoidal shapes with respect of
fission and, as a consequence, the boundary of the transition from 
ellipsoidal to toroidal shapes and the two-proton and two-neutron 
drip lines for toroidal  nuclei.  However, they will not affect the 
general features.






















§ CONCLUSIONS




      In conclusion, the detailed investigation of the properties of hyperheavy 
nuclei has been performed in the framework of covariant density functional
theory. The following conclusions have been obtained:




  * The stability of spherical hyperheavy nuclei located in the centers
of the  (Z∼ 138, N∼ 230),  (Z∼ 156, N∼ 310) and 
(Z∼ 174, N∼ 410) islands of stability with respect of triaxial
and octupole distortions has been established in the RHB+OCT 
and TRHB calculations.


  * Proton and neutron densities, central depressions in these densities,
charge radii, and neutron skins of the nuclei  located in the centers of 
these islands of stability have been investigated  in detail. Obtained 
results clearly indicate that the accuracy of the reproduction of charge 
radii and neutron skins by the CEDFs could be an important criteria in 
favoring or disfavoring the predictions of one or another functional for 
the islands of stability of spherical hyperheavy nuclei. Among considered 
functionals, the DD-ME2 and DDPC1 functionals provide the best global 
description of charge radii and predict the highest fission barriers in these 
regions.  The results obtained in future PREX-2 experiment on neutron skin 
in ^208Pb <cit.> would be quite useful in helping to 
discriminate the predictions of different functionals for the islands of  
stability of spherical hyperheavy nuclei.
 

  * Underlying shell structure of the nuclei  located in the centers of these islands 
of  stability has been investigated  in detail. Large neutron shell gaps at 
N=228, 308 and 406 have a sizeable impact on two neutron-separation 
energies. On the other hand, large proton gap appear only at Z=154 in
the (Z∼ 156, N∼ 310) island of stability of spherical hyperheavy 
nuclei.  As a result, this is the largest island of stability of spherical superheavy
nuclei found in the calculations. No significant proton gaps are seen in other 
two islands of stability. Taking into account clear localization of the islands of 
stability of spherical hyperheavy nuclei in the (Z, N) plane these features 
strongly suggest that the shell effects at deformed shapes leading to negative 
shell correction energies at some deformation and thus to fission barriers play 
an  important role in the stabilization of spherical hyperheavy nuclei.



  * The shape evolution of toroidal shapes along the fission path and the stability 
of  such shapes with respect of fission have been studied. In considered cases,
the analysis shows the transition from symmetrical toroid (at the local minimum) 
to the asymmetric one (at the saddle point). This transition cost the energy which 
is a physical reason for the formation of fission barrier and, thus, for the stability 
of such shapes.


  * The topology of potential energy surfaces for ellipsoidal shapes of the super- and
hyperheavy nuclei has been compared. In both types of the nuclei the PES has
the form of  the 
canyon in which some local valleys and hills are located.  The canyon is quite 
narrow in superheavy nuclei which prevents the formation of local minima at 
large oblate deformation and limits the role of triaxial deformation. On the contrary,
this canyon is much wider in hyperheavy nuclei with a mountain, centered at 
β_2 ∼ 0, formed in it. This leads to the formation of local minima at oblate 
superdeformation,  increased role of triaxiality and higher instability of hyperheavy 
nuclei against fission as compared with superheavy ones.


  * The extension of nuclear landscape  to hyperheavy nuclei  with proton numbers up 
to Z=180 has been performed.   With increasing proton number beyond Z ∼ 130 
the transition from ellipsoidal-like nuclear shapes to toroidal shapes takes place in axial
RHB calculations. The ellipsoidal ground states are affected by above-mentioned 
increased instability against fission.  Many hyperheavy nuclei with toroidal shapes (as 
the lowest in energy solutions in axial RHB calculations) are expected to be unstable 
towards multi-fragmentation.  However, it is difficult to quantify their stability or instability
since the description of toroidal shapes requires the basis which is typically significantly 
larger than the one employed for the description of ellipsoidal-like shapes. This makes
the calculations with octupole or triaxial deformation included impossible for toroidal 
shapes with extreme β_2 values. Nevertheless,  three islands of stability of spherical 
hyperheavy nuclei are predicted. The nuclei in these  islands will become the ground states in the case of instability of relevant toroidal states.




   Detailed investigation of possible mechanisms of the creation of spherical 
and toroidal hyperheavy nuclei represents an interesting topic but goes beyond the 
scope of the present manuscript. The nuclei in the (Z∼ 138, N∼ 230),  
(Z∼ 156, N∼ 310) and (Z∼ 174, N∼ 410) islands of stability of
spherical hyperheavy nuclei have neutron to proton ratios of N/Z ∼ 1.67, 
N/Z ∼ 1.99 and N/Z ∼ 2.36, respectively. Thus, they cannot be formed 
in laboratory conditions and the only possible environment in which they can 
be produced is the ejecta of the mergers of neutron stars <cit.>.
In a similar fashion, the regions of neutron stars with nuclear pasta phases 
<cit.> may be a breading ground for the formation of 
toroidal nuclei in the ejecta of the merger of neutron stars.  
The two-proton drip line for toroidal nuclei is characterized by neutron to proton ratio 
of N/Z ∼ 1.25. Thus, the stability and/or multi-fragmentation of toroidal nuclei located 
in the vicinity of two-proton drip line could possibly be studied in  nucleus-nucleus collisions 
of stable nuclei or the nuclei located close to the beta-stability line. This is similar to what 
has alrelsady been done  in the ^86Kr + ^93Nb reaction at incident energies ranging 
from 35 to 95 MeV/nucleon in Ref. <cit.>.


   Any extrapolation beyond known regions of nuclear chart in which the functionals have 
been fitted is associated with theoretical uncertainties <cit.>. This 
is especially  true for the present study  with its significant extrapolations  in proton and 
neutron numbers. Despite the fact that our study is mostly based on the DD-PC1 CEDF which, according to
the results of global studies, is considered to be  the best relativistic functional, in no
way the extrapolations based on it should be considered as completely safe.  This is also
true for any relativistic or non-relativistic functional.  However, when calculated effects 
are substantially larger than the expected theoretical uncertainties one can speak about reliable theoretical 
predictions. For example, the predicted transition from ellipsoidal to toroidal shapes with 
increasing proton number is a solid prediction.  Note that according to Ref. <cit.>
it appears also in  Gogny DFT as exemplified by the calculated results for two hyperheavy
nuclei. On the other  hand, the borderline in the (Z,N) plane between these two types of 
the shapes is expected to depend on the details of the functional since it is defined by the 
fission properties which are subject of appreciable theoretical uncertainties <cit.>. 
However, the lowering of the fission barrier heights for ellipsoidal shapes with increasing 
proton number Z, which defines this boundary,  appears both in relativistic and non-relativistic 
models (see Sect. 
<ref> in the present paper and Refs. <cit.>).  In addition, the 
size of the regions of possible stability of spherical hyperheavy nuclei  and the stability of the
nuclei in these regions depend on the functional (see Ref. <cit.>) and on the details 
of underlying shell structure (see Sect. <ref> in the present paper).  The latter
is subject of appreciable theoretical uncertainties when extrapolations are performed to
unknown regions in the (Z,N) plane which are located far away from known part of 
nuclear chart in which the functionals have been fitted <cit.>.  In addition, 
there could be hidden biases in the CDFT which could affect model predictions.  In such a 
situation, detailed investigations of hyperheavy nuclei in the framework of non-relativistic 
density functional theories based on the Skyrme and Gogny forces would be extremely 
useful for an understanding of general structure of hyperheavy nuclei, the transition from 
ellipsoidal to toroidal shapes with increasing proton number and possible existence of the 
islands of stability of spherical hyperheavy nuclei.  They will also allow to compare the 
predictions obtained in relativistic and non-relativisitc models.



 




§ ACKNOWLEDGMENTS



 This material is based upon work supported by the U.S. Department 
of Energy, Office of Science, Office of Nuclear Physics under Award 
No. DE-SC0013037 and by the Department of Energy National Nuclear 
Security Administration under Award No. DE-NA0002925.


