Hyperheavy nuclei: existence and stability

What are the limits of the existence of nuclei? What are the highest proton numbers $Z$ at which the nuclear landscape and periodic table of chemical elements cease to exist? These deceivably simple questions are difficult to answer especially in the region of hyperheavy ($Z\geq 126$) nuclei. We present the covariant density functional study of different aspects of the existence and stability of hyperheavy nuclei. For the first time, we demonstrate the existence of three regions of spherical hyperheavy nuclei centered around ($Z\sim 138, N\sim 230$), ($Z\sim 156, N\sim 310$) and ($Z\sim 174, N\sim 410$) which are expected to be reasonably stable against spontaneous fission. The triaxiality of the nuclei plays an extremely important role in the reduction of the stability of hyperheavy nuclei against fission. As a result, the boundaries of nuclear landscape in hyperheavy nuclei are defined by spontaneous fission and not by the particle emission as in lower $Z$ nuclei. Moreover, the current study suggests that only localized islands of stability can exist in hyperheavy nuclei.

The investigation of superheavy elements (SHE) remains one of the most important sub-fields of low-energy nuclear physics [1]. The element Og with proton number Z = 118 is the highest Z element observed so far [2]. Although future observation of the elements in the vicinity of Z ∼ 120 seems to be feasible, this is not a case for the elements with Z beyond 122. Considering also that the highest in Z spherical shell closure in SHE is predicted at Z = 126 in Skyrme density functional theory (DFT) [3], it is logical to name the nuclei with Z > 126 as hyperheavy [4,5]. The properties of such nuclei are governed by increased Coulomb repulsion and single-particle level density; these factors reduce the localization of shell effects in particle number [5].
Although hyperheavy nuclei have been studied both within DFTs [4,5,6,7,8,9] and phenomenological [10,11,12] approaches, the majority of these studies have been performed only for spherical shapes of the nuclei. This is a severe limitation which leads to misinterpretation of physical situation in many cases since there is no guarantee that spherical minimum in potential energy surface exist even in the nuclei with relative large spherical shell gaps (see discussion in Ref. [13]). In addition, the stability of hyperheavy nuclei against spontaneous fission could not be established in the calculations restricted to spherical shape. The effects of axial and triaxial deformations in hyperheavy nuclei are considered only in Refs. [6,14] and in Ref. [8,15], respectively. However, only few nuclei are studied in Refs. [6,14,15] and according to the present study the deformation range employed in Ref. [8] is not sufficient for Z ≥ 130 nuclei.
The investigation of hyperheavy nuclei is also intimately connected with the establishments of the limits of both the nuclear landscape and periodic table of elements. The limits of nuclear landscape at the proton and neutron drip lines and related theoretical uncertainties have been extensively investigated in a number of theoretical frameworks but only for the Z < 120 nuclei [16,17,18]. The atomic relativistic Hartree-Fock [19] and relativistic Multi-Configuration Dirac-Fock [20,21] calculations indicate that the periodic table of elements terminates at Z = 172 and Z = 173, respectively. However, at present it is not even clear whether such nuclei are stable against fission. In addition, Refs. [19,20,21] employ phenomenological expression for charge radii and its validity for the Z ∼ 172 nuclei is not clear.
To address these deficiencies in our understanding of hyperheavy nuclei the systematic investigation of even-even nuclei from Z = 122 up to Z = 180 is performed within the axial relativistic Hartree-Bogoliubov (RHB) framework employing the DD-PC1 covariant energy density functional [22]. This functional provides good description of the ground state and fission properties of known even-even nuclei [18,23]. To establish the stability of nuclei with respect to triaxial distortions a number of nuclei have been studied within the triaxial RHB [24] and relativistic mean field + BCS (RMF+BCS) [25] frameworks. The main goals of this study are (i) to understand whether the nuclei stable against fission could be present in the Z ≥ 126 region and (ii) to define the most important features of such nuclei.
The CDFT calculations are performed within the relativistic Hartree-Bogoliubov framework [26] employing the state-ofthe-art covariant energy density functionals the global performance of which with respect to the description of the ground state [13,18,27] and fission [23,24,25,28] properties is well   established. These functionals (DD-PC1 [22], DD-ME2 [29], PC-PK1 [30] and NL3* [31]) represent the major classes of covariant density functional models [18]. The absolute majority of the calculations employ the DD-PC1 functional which is considered to be the best relativistic functional today based on systematic and global studies of different physical observables [13,18,24,27,28,32]. Other functionals are used to assess the systematic theoretical uncertainties in the predictions of the heights of fission barriers around spherical minima. To avoid the uncertainties connected with the definition of the size of the pairing window [33], we use the separable form of the finiterange Gogny pairing interaction introduced in Ref. [34].
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 (and up to N B = 20 bosonic shells for the meson fields in meson exchange functionals) are taken into account. The comparison of the axial RHB calculations with N F = 20 and N F = 30 shows that in 208 Pb the truncation of basis at N F = 20 provides sufficient accuracy for all deformations of interest. 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. 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 nu-clei) 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 Figs. 1c and d for the 466 156 and 426 176 results); their description requires N F = 30. Thus, the truncation of basis is made dependent on the nucleus and typical profile of deformation energy curves or potential energy surfaces.
The deformation parameters β 2 and γ are extracted from respective quadrupole moments: via where R 0 = 1. 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. Note that physical observables are frequently shown as a function of the Q 20 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 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 [18]; 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. The nuclei up to Z = 138 are calculated using the basis with N F up to 26. On the contrary, with the exception of the 466 156 and 426 176 nuclei, the Z = 140 − 180 nuclei are calculated only with N F = 20. The major goals of the calculations for the Z = 140 − 180 nuclei are (i) to define the type of the LEMAS states, (ii) to find whether spherical or normal deformed states could be the LEMAS states of these nuclei and (iii) to calculate the fission barriers around spherical states.
The required size of the basis limits the applicability of triaxial calculations to typically |β 2 | < 2 range. The nuclei with the ground states located at the deformations below β 2 ∼ 1.0 are calculated in triaxial RHB framework [24], while a pair of nuclei with local minima at β ∼ 2.4, γ ∼ 60 • corresponding to toroidal shapes were calculated in triaxial RMF+BCS framework [25]. The later framework is more numerically stable at very large β 2 values. Because of high computational cost of the calculations with triaxiality included, only limited number of nuclei were studied in these frameworks. The role of octupole deformation in the nuclei shown in Fig. 6 has been studied in the axial reflection asymmetric RHB code of Ref. [32]. These calculations are performed with N F = 20. Fig. 1 illustrates the dependence of the deformation energy curves, obtained in axial RHB calculations, on the nucleus. The Z = 82 208 Pb nucleus is spherical in the ground state. The total energy of the nucleus is increasing rapidly with increasing oblate deformation. On the prolate side, it increases with the increase of quadrupole deformation up to β 2 ∼ 1.4 and then stays more or less constant. This leads to the existence of high (∼ 30 MeV) and very broad fission barrier which is responsible for the stable character of this nucleus.
The 354 134 nucleus shows completely different profile of the deformation energy curves (Fig. 1b). The LEMAS is located at β 2 ∼ −0.5 and the deformation energy curves on the oblate side are more flat in energy as compared with 208 Pb. The fission barrier for the β 2 ∼ −0.5 minimum is rather high (∼ 8. Further increase of proton number leads to drastic modifications of the deformation energy curves. In the 466 156 and 426 176 nuclei, the minimum appears at extreme β 2 ∼ −4.0 values. However, these minima are potentially unstable with respect to the transition to the prolate shape via γ-plane and subsequent fission since prolate shapes with corresponding quadrupole deformations are located at lower energies (compare dashed lines with solid ones in Figs. 1c and d). Note also that in the 466 156 nucleus there are excited local β 2 ∼ −1.2 and spherical minima which could be potentially stable against fission.
The evolution of the neutron density distributions with the change of the β 2 value are shown for the 466 156 nucleus in Fig.  2. The nucleus at spherical shape is characterized by the density depression in the central part of the nucleus; the maximum neutron density ρ = 0.0896 fm −3 is achieved at radial coordinate r = 6.55 fm while the density in the center is only ρ = 0.076 fm −3 . This depression is similar (but less pronounced) to the one predicted for the 292 120 superheavy nucleus in Refs. [3,35]. Our calculations show neither bubble nor semi-bubble shapes (in the language of Ref. [4]) for the lowest in energy solutions of spherical nuclei shown in Fig. 6 below. Note that proton density is roughly half of the neutron one and central density depression is somewhat more pronounced in proton subsystem as compared with neutron one. As illustrated in Fig. 2b, biconcave disk density distribution is formed at large oblate deformation of β 2 = −1.0. Further decrease of the β 2 values leads to the formation of toroidal shapes (Figs. 2c and d). It is observed that with the increase of absolute value of β 2 the radius of the toroid increases and the tube radius decreases.
The biconcave disk and toroidal shapes in atomic nuclei have been investigated in a number of the papers [14,15,37,38,36,39]. However, in absolute majority of the cases such shapes correspond to highly excited states either at spin zero [15,36] or at extreme values of angular momentum [37,38,40]. The latter substantially exceed the values of angular momentum presently achievable at the state-of-art experimental facilities [41]. The competition of such shapes at spin zero in superheavy eveneven Z = 120 isotopes with N = 166 − 190 and in the eveneven N = 184 isotones with Z = 106 − 124 has been investigated in constrained Skyrme-HFB calculations in Ref. [36]. It was concluded that investigated nuclei in toroidal shapes are unstable against returning to the shape of sphere-like geometry (Ref. [36]). Similar study for superheavy 316 122, 340 130, 352 134 and 364 138 nuclei has been performed in Skyrme Hartree-Fock calculations of Ref. [15]; only in 364 138 nuclei the toroidal solution is the lowest in energy. The Gogny HFB calculations of Ref. [14] showed that toroidal shapes represent the lowest in energy solutions at axial shape in the 416 164 and 476 184 nuclei. However, it is well known that triaxial deformation lowers the fission barriers in actinides and superheavy nuclei with Z ≤ 120 and N ≤ 184 [25,43,44,45,46,47]. These nuclei are either prolate or spherical in their ground states and thus the impact of triaxiality is limited: for example, the lowering of inner fission barriers in actinides due to triaxiality is typically on the level of 1-3 MeV. On the contrary, the impact of triaxiality on fission barriers gets much more pronounced in the nuclei with ground state oblate shapes and it generally increases with the rise of their oblate deformation. 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.
These features are illustrated in Fig. 4. The 360 130 nucleus is an example of the coexistence of spherical ground state and excited (at 0.8 MeV) oblate (with β 2 ∼ −0.5) minimum. The static fission paths from these minima are comparable in length and both of them have reduced (by ∼ 2 MeV) inner fission barriers as compared with axial RHB calculations (see supplementary Table 1). The effect of the reduction of inner fission barrier due to triaxiality becomes much more pronounced in the 432 134 nucleus. As compared with axial calculations, the presence of triaxiality leads to the shift of minimum from (β 2 ∼ 0.74, γ = 60 • ) to (β 2 ∼ 0.82, γ ∼ 37 • ) and the reduction of the fission barrier height from 8.  Table 1). The triaxiality leads to the γ-softness of potential energy surfaces so that these minima drift in the γ-plane by 10 − 15 • . However, it also leads to substantial reduction of fission barrier heights down to ∼ 2 MeV (see supplementary  Table 1). In axial RHB calculations, the 392 134 nucleus has superdeformed oblate ground state with β 2 = −0.79 and highly excited (at excitation energy of 2.69 MeV) oblate state with β 2 = −0.23. The fission barriers for these two minima are 10.24 and 7.55 MeV, respectively. The triaxiality substantially affects the position of first minimum so it drifts to β 2 = 0.88, γ = 39 • but has almost no effect on the second minimum. However, it has huge impact on the heights of their fission barriers which are reduced to 0.56 and 2.08 MeV, respectively (see supplementary  Table 1).
Supplementary Table 1 summarizes the results of more systematic triaxial RHB calculations. 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. As a result, the majority of these nuclei would be unstable with respect to fission. Similar trend of the evolution of fission barriers with proton number has also been seen in microscopic+macroscopic (mic+mac) calculations with Woods-Saxon potential and Skyrme DFT calculations with the SLy4 functional presented in Ref. [8]. Note that these calculations use smaller deformation plane (ranging from β 2 = −0.85 up to β 2 = 0.45) as compared with the one shown in Fig. 4 . The Skyrme DFT calculations provide higher fission barriers as compared with mic+mac and our RHB results. However, the SLy4 functional substantially overestimates fission barriers in actinides and SHE [8].
The situation however is substantially complicated by the fact that with increasing proton number toroidal shapes correspond to the lowest in energy solutions in axial RHB calculations (Fig. 3). Their large β 2 values and high Z and N values require increased basis which makes triaxial RHB and RMF+BCS calculations prohibitively time consuming. A priori we cannot exclude the stability of such shapes against fission or multifragmentation. This is illustrated by the calculations of the 354 134 (Fig. 5) and 348 138 (supplementary Fig. 3) nuclei, for which the N F = 20 basis provides acceptable numerical accuracy. In these nuclei, the oblate minimum with β 2 ∼ −2.5, β 4 ∼ −4.4 is unstable with respect to triaxial distortions (left panels of these figures). On the contrary, the excited β 2 ∼ −2.3, β 4 ∼ +1.5 minimum is stable with respect to triaxial distortions (see right panels of Fig. 5 and supplementary Fig.  3).
The triaxial RHB calculations for the |β 2 | ≤ 1.0 part of the deformation plane clearly indicate the general trend of the reduction of the stability of the minima located at these deformations with respect to fission with increasing proton number. The triaxial RMF+BCS calculations also indicate the potential stability of toroidal shapes located in the minima with β 2 < −2.0 with β 4 > 0. Unfortunately, the systematic triaxial calculations of the stability of such minima are beyond available computational power. Thus, their more detailed investigation is left for future.
Note also that toroidal nuclei are expected to be unstable against multifragmentation [42,51]. The most detailed investigation of the instabilities of toroidal nuclei with respect of socalled breathing and sausage deformations has been performed in Ref. [51]. The breathing deformation preserves the azimuthal symmetry of the torus and it is defined by the radius of torus and the radius of its tube. In our calculations, this type of deformation is related to the β 2 values (see discussion of Fig. 2 above).
The results of Ref. [51] clearly indicate the stability of toroidal nuclei with respect of breathing deformation both in liquid-drop model calculations and in Strutinsky type calculations. This is also the case in our calculations which show minima at large negative β 2 values in deformation energy curves presented as a function of β 2 (see Fig. 1b,c, and d). The sausage deformations make a torus thicker in one section(s) and thinner in another section(s); they are examplified by the density distributions shown in supplementary Fig. 2. The analysis of Ref. [51] clearly indicates the instability of toroidal nuclei with respect of sausage deformations in the liquid drop model. However, it was not excluded in Ref. [51] that the instability in the sausage degree of freedom may be counterbalanced by shell effects at some combinations of proton and neutron numbers and deformations. The situation here is similar to superheavy nuclei which are unstable in liquid drop model. The instability with respect of sausage deformations has not been studied so far in either Strutinsky type models [51] or in density functional theories. However, our results for the β 2 ∼ 2.5, β 4 ∼ −4.4, γ = 60 • solutions in the 354 134 and 348 138 nuclei show for the first time this type of instability also in the framework which takes shell effects into account. The analysis of the deformation energy curves obtained in axial RHB calculations reveals that hyperheavy nuclei could be stabilized at spherical shapes in some regions (see the insert to Fig. 1c). If the toroidal shapes in these nuclei are unstable against triaxial distortions or multifragmentation, these states represent the ground states. From our point of view, this is the most likely scenario. Otherwise, they are excited states frequently located at high excitation energies (Fig. 1c). It was verified that these spherical states are stable with respect to triaxial and octupole distortions. The largest island of stability of spherical hyperheavy nuclei is centered around Z ∼ 156, N ∼ 310 (Fig. 6a). In the calculations with the DD-PC1 functional the fission barriers reach 11 MeV for the nuclei located in the center of the island of stability. This is substantially larger than the fission barriers predicted in the CDFT for experimentally observed superheavy nuclei with Z ∼ 114, N ∼ 174 for which calculated inner fission barriers are around 4-5 MeV [24]. Smaller islands of stability of spherical hyperheavy nuclei are predicted at Z ∼ 138, N ∼ 230 and Z ∼ 174, N ∼ 410 (Fig. 6a). Since nuclei in these three regions have N/Z ≥ 1.64 they cannot be formed in laboratory conditions. The only possible environment in which they can be produced is the ejecta of the mergers of neutron stars [48].
Additional calculations have been performed with the DD-ME2 [29], PC-PK1 [30] and NL3* [31] functionals in order to evaluate systematic theoretical uncertainties [49] in the predictions of fission barriers for spherical hyperheavy nuclei. The DD-ME2 functional provides predictions comparable with the DD-PC1 one (Fig. 6a,b). In contrast, the PC-PK1 and NL3* functionals predict lower fission barriers and smaller regions of stability (Fig. 6c,d). Note that the nuclear matter properties and the density dependence are substantially better defined for density-dependent (DD*) functionals as compared with nonlinear NL3* and PC-PK1 ones [27]. As a consequence, they are expected to perform better for large extrapolations from known regions. The large fission barriers obtained in the densitydependent functionals will lead to substantial stability of spherical hyperheavy nuclei against spontaneous fission. This stability is substantially lower for the NL3* and PC-PK1 functionals.
Note that these spherical states are also relatively stable against α-decay (see supplementary Fig. 4). Theoretical uncertainties in the predictions of the α-decay half-lives due to the use of different empirical formulas for their calculations and the CEDFs are evaluated for the Z ∼ 156, N ∼ 310 region of spherical hyperheavy nuclei in supplementary Figs. 5 and 6, respectively. One can see that when combined these uncertainties could reach 10 orders of magnitude in the center of region. However, even with these uncertainties accounted the α-decay half-lives of many nuclei are substantial exceeding seconds, hours and days ranges. Considering empirical nature of the formulas employed more microscopic studies of the α-decay half-lives would be highly desirable. It is also important in future to investigate other competing decay modes such as cluster [52,53] and β [54,55] decays to fully establish the potential stability of spherical hyperheavy nuclei.
Existing atomic calculations suggest that the periodic table of elements ends at Z ∼ 172 [19,20,21]; this takes place when the 1s electron binding energy dives below −2mc 2 . However, these calculations employ the empirical formulas for the rootmean-square (RMS) nuclear charge radii. For example, the calculations of Ref. [19] employ the formula from Ref. [50] which underestimates the RMS nuclear charge radii as compared with the ones obtained in the RHB calculations. This is exemplified by the values of RMS nuclear charge radii in the 368 138, 466 156 and 584 174 nuclei which are 6.52 fm (6.91 fm), 7.10 fm (7.576 fm), 7.62 fm (8.312 fm) in the calculations with empirical formula of Ref. [50] (the RHB calculations with DD-PC1). Note that these nuclei represent the centers of the islands of stability of spherical hyperheavy nuclei (see Fig. 6a). Unfortunately, the impact of nuclear size changes on atomic properties and thus on the end of periodic table of elements has not been investigated in Refs. [19,20,21]. However, these differences in the RMS nuclear charge radii are substantial and new atomic calculations are needed to see how they can affect the end of periodic table of elements.
In  174, N ∼ 410). However, theoretical systematic uncertainties in the predictions of their fission barriers are substantial. These results clearly indicate that the boundaries of nuclear landscape in hyperheavy nuclei are defined by spontaneous fission and not by particle emission as in lower Z nuclei. Moreover, the current study suggests that only localized islands of stability can exist in hyperheavy nuclei.

Acknowledgements
This material is based upon work supported by the Department of Energy National Nuclear Security Administration under Award No. de-na0002925 and by the US Department of Energy, Office of Science, Office of Nuclear Physics under Award No. de-sc0013037. MeV (with respect of the minimum) is shown by black cross. The energy difference between two neighboring equipotential lines is equal to 5 MeV and 2 MeV in left and right panels, respectively. The same energy minimum is used for colormap in both panels. Here β min , β saddle and E B ax are the equilibrium quadrupole deformation of the global (local) minimum, the quadrupole deformation and the energy of the saddle along respective fission path. The excited minima are indicated by asterisks (*). Their excitation energies are shown in brackets in column 3. The results of the triaxial RHB calculations are provided in the columns 6 − 8. Note that the allowance of triaxial deformation could shift the position of the local minimum in the deformation plane and in absolute majority of the cases shifts the positions of the saddle points. Thus, (β, γ) Figure 11: The evaluation of theoretical uncertainties in the predictions of α-decay half-lives emerging from the use of different empirical formulas. The calculations are performed with the DD-PC1 functional. The results obtained using Viola-Seaborg formula with the parameters defined in Ref. [2] and Ref. [3] are labelled as VS-1 and VS-2, respectively. Note that employed CEDFs describe well experimental α-decay half-lives of even-even actinides and superheavy nuclei when Viola-Seaborg formula is used [6]. The labels 'Royer' and 'mod-Brown' are used for the results calculated with Royer formula of Ref. [5] and modified Brown formula of Ref. [4], respectively. The results for a given isotope chain obtained with different empirical formulas are shown by the same type of line. Panel and Z = 160 (dotted lines) isotope chains of the Z ∼ 156, N ∼ 310 region of spherical hyperheavy nuclei. They are obtained using Viola-Seaborg formula with the parameters defined in Ref. [2]. The range of neutron numbers corresponds to the one shown in Fig. 6(a-d) of the manuscript.