Extension of nuclear landscape to hyperheavy nuclei

TL;DR

This study extends the nuclear landscape to hyperheavy nuclei using covariant density functional theory, revealing new islands of stability and the impact of shape and triaxiality on nuclear properties.

Contribution

First prediction of three islands of stability for spherical hyperheavy nuclei and detailed analysis of their shell structure and stability factors.

Findings

Identification of three stability islands at Z~154, N=228, 308, 406

Large neutron shell gaps define island centers

Triaxiality significantly affects fission barrier heights

Abstract

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. [1]. 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. 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.

Tables

Table 1. Table 1: Charge radii r c ​ h subscript 𝑟 𝑐 ℎ r_{ch} [in f ​ m 𝑓 𝑚 fm ] and neutron skins r s ​ k ​ i ​ n subscript 𝑟 𝑠 𝑘 𝑖 𝑛 r_{skin} [in f ​ m 𝑓 𝑚 fm ] of the density distributions shown in Fig.  4 .

Z	N		DD-PC1	DD-ME2	NL3*	PC-PK1
Pb	126	$r_{ch}$	5.513	5.518	5.509	5.519
		$r_{skin}$	0.202	0.193	0.288	0.257
120	172	$r_{ch}$	6.272	6.282	6.276	6.286
		$r_{skin}$	0.104	0.091	0.164	0.141
138	230	$r_{ch}$	6.759	6.765	6.799	6.811
		$r_{skin}$	0.198	0.188	0.283	0.249
156	310	$r_{ch}$	7.330	7.326	7.402	7.420
		$r_{skin}$	0.290	0.295	0.427	0.364
174	410	$r_{ch}$	7.927	7.930	8.071	8.087
		$r_{skin}$	0.440	0.466	0.616	0.520

Table 2. Table 2: The heights of the fission barriers along the fission paths from different minima obtained in axial and triaxial RHB calculations. The columns 3 − 5 3 5 3-5 show the results of the axial RHB calculations. Here β m ​ i ​ n subscript 𝛽 𝑚 𝑖 𝑛 \beta_{min} , β s ​ a ​ d ​ d ​ l ​ e subscript 𝛽 𝑠 𝑎 𝑑 𝑑 𝑙 𝑒 \beta_{saddle} and E a ​ x B superscript subscript 𝐸 𝑎 𝑥 𝐵 E_{ax}^{B} are the equilibrium quadrupole deformation of the global minimum, the quadrupole deformation and the energy of the saddle along respective fission path. The results of the triaxial RHB calculations are provided in the columns 6 − 8 6 8 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, ( β , γ ) m ​ i ​ n subscript 𝛽 𝛾 𝑚 𝑖 𝑛 (\beta,\gamma)_{min} , ( β , γ ) s ​ a ​ d ​ d ​ l ​ e subscript 𝛽 𝛾 𝑠 𝑎 𝑑 𝑑 𝑙 𝑒 (\beta,\gamma)_{saddle} and E t ​ r ​ i ​ a ​ x B superscript subscript 𝐸 𝑡 𝑟 𝑖 𝑎 𝑥 𝐵 E_{triax}^{B} show the deformations of the minima, the deformations of saddle points and their energies obtained in triaxial RHB calculations. The neutron numbers of the nuclei in which superdeformed minimum with β 2 ∼ 0.5 similar-to subscript 𝛽 2 0.5 \beta_{2}\sim 0.5 is lower than normal-deformed or spherical ones are marked by an asterisk. With exception of these nuclei, the values in parentheses show either the deformation of superdeformed minimum or the deformation of the saddle of outer fission barrier or the height of outer fission barrier. Note that these values are shown only when the lowest height of the outer fission barrier obtained in the TRHB and RHB+OCT calculations is higher than 2 MeV. The column 3 shows the Δ ​ E g ​ a ​ i ​ n = E a ​ x B − E t ​ r ​ i ​ a ​ x B Δ superscript 𝐸 𝑔 𝑎 𝑖 𝑛 superscript subscript 𝐸 𝑎 𝑥 𝐵 superscript subscript 𝐸 𝑡 𝑟 𝑖 𝑎 𝑥 𝐵 \Delta E^{gain}=E_{ax}^{B}-E_{triax}^{B} quantity which is the decrease of the height of respective fission barrier due to triaxiality.

		Axial RHB			Triaxial RHB
$Z$	$N$	${\beta }_{min}$	${\beta }_{saddle}$	$E_{ax}^B$	${(\beta ,\gamma )}_{min}$	${(\beta ,\gamma )}_{saddle}$	$E_{triax}^B$	$\Delta E^{gain}$
1	2	3	4	5	6	7	8	9
106 (Sg)	152	0.29	0.57	10.09	0.29, 0	0.62, 12.4	7.04	3.05
	162	0.26	0.65	10.70	0.26, 0	0.68, 8.4	7.48	3.32
	172	0.14	0.69	5.31	0.14, 0	0.71, 3.5	2.95	2.36
	182	-0.05 (0.49)	0.27 (0.73)	4.25 (3.42)	0.05, 60.0 (0.49, 0)	0.47, 23.6 (0.81, 8.0)	3.70 (2.47)	0.55 (0.95)
	192	0.39	0.59	2.42	0.40, 0	0.61, 6.3	2.20	0.22
	202	0.28	0.59	4.27	0.29, 9.5	0.66, 13.2	3.73	0.54
	212	0.25	0.54	7.23	0.25, 0	0.69, 10.3	4.99	2.24
	222	0.25	0.55	8.53	0.250, 0	0.70, 8.7	4.49	4.04
	232	0.23	0.65	6.73	0.23, 0	0.62, 10.4	4.63	2.10
	242	0.13 (0.45)	0.25 (0.65)	3.77 (5.08)	0.13, 0 (0.45, 0)	0.25, 0 (0.60, 6.9)	3.77 (2.37)	0.0 (2.71)
	252	-0.06 (0.45)	0.25 (0.69)	6.40 (5.68)	0.06, 59.1 (0.45, 0)	0.42, 26.7 (0.69,10.3)	3.86 (3.20)	2.54 (2.48)
110 (Ds)	162	0.24	0.66	8.98	0.242,0	0.65, 6.9	6.20	2.78
	172	0.15 (0.46)	0.30 (0.70)	3.97 (5.75)	0.15, 0 (0.45, 0)	0.30, 0.0 (0.74, 5.6)	3.92 (3.24)	0.05 (2.51)
	182	-0.14	0.26	3.92	0.139, 60	0.36, 35.6	2.54	1.38
	192	0.41	0.58	2.52	0.417,0	0.62, 5.8	2.52	0.0
	202	0.38	0.56	4.56	0.385,0	0.73, 14.4	2.90	1.66
	212	0.26	0.51	5.95	0.262,0	0.54, 15.7	4.28	1.67
	222	0.24	0.54	7.02	0.243,0	0.36, 28.5	5.90	1.12
	232	0.22	0.60	6.27	0.220,0	0.61, 7.0	4.51	1.76
	242	0.14 (0.47)	0.27 (0.66)	3.92 (5.06)	0.17, 0 (0.46, 0)	0.27, 0 (0.68, 6.0)	3.92 (2.51)	0.0 (2.55)
	252	0.44	0.70	4.79	0.444,0	0.72, 11.6	2.27	2.52
114 (Fl)	162	0.23 (0.52)	0.40 (0.66)	5.95 (6.00)	0.23, 0 (0.54, 0)	0.40, 0 (0.57, 5.1)	5.95 (3.85)	0.0 (2.15)
	172	0.15 (0.50)	0.32 (0.73)	5.22 (4.76)	0.15, 0 (0.51, 0)	0.32,0 (0.73, 6.0)	5.22 (3.33)	0.0 (1.43)
	182	-0.14	0.26	4.26	0.14, 60	0.41, 38.5	2.72	1.54
	192	-0.38	0.15	3.79	0.40, 42	0.42, 33.4	1.35	2.44
	202	0.38	0.54	2.79	0.38, 0	0.54, 3.7	2.71	0.08
	212	0.27	0.49	4.27	0.28, 0	0.40, 23.1	3.31	0.96
	222	0.24	0.41	5.56	0.24, 0	0.36, 26.3	4.95	0.61
	232	0.21	0.35	4.68	0.21, 0	0.35, 0	4.68	0.0
	242	0.14	0.27	3.73	0.142,0	0.27,0	3.73	0.0
	252*	0.458 (0.0)	0.72 (0.23)	3.31 (4.68)	0.459,0 (0.0, 0)	0.63, 8.0 (0.46, 23.9)	2.27 (1.40)	1.04 (3.28)
118 (Og)	172	0.0	0.34	5.11	0.000,0	0.49, 28.2	4.47	0.64
	182*	0.57 (-0.27)	0.84 (0.26)	4.02 (5.32)	0.58,0 (0.27, 60)	0.66, 9.5 (0.42, 36.4)	2.41 (3.57)	1.61 (1.75)
	192	-0.39	0.15	5.79	0.40, 74	0.33, 43.3	2.37	3.42
	202	-0.43	0.07	6.24	0.43, 60	0.46, 45.0	1.31	4.93
	212	0.29	0.44	3.03	0.30, 11	0.64, 21.3	2.45	0.58
	222	0.24	0.39	5.19	0.24, 0	0.37, 27.2	4.18	1.01
	232	0.22	0.35	4.89	0.22, 0	0.35, 0	4.56	0.33
	242	-0.20	0.31	3.70	0.21, 60.0	0.47, 32.4	2.88	0.82
	252	-0.19	0.19	5.06	0.20, 60	0.41, 38.5	2.98	2.08
	262	-0.23	0.15	4.95	0.24, 59.9	0.32, 41.6	0.80	4.15
	272	-0.49	-0.01	5.98	0.49, 59.6	0.48, 51.7	0.44	5.54

Equations

$$E_{RHB/RMF+BCS}+\sum _{\lambda ,\mu }{C}_{\lambda ,\mu }(⟨{\hat{Q}}_{\lambda ,\mu }⟩-q_{\lambda ,\mu }{)}^2.$$

$${\hat{Q}}_{20}$$

$${\hat{Q}}_{22}$$

$${\hat{Q}}_{30}$$

$$⟨{\hat{Q}}_{10}⟩=0$$

$$Q_{20}$$

$$Q_{22}$$

$$Q_{30}$$

$${\beta }_2$$

$$\gamma$$

$${\beta }_3$$

$$V({\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{r}}_1^{\prime },{\mathbf{r}}_2^{\prime })$$

$$=-G\delta (\mathbf{R}-$$

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

$${\Delta }_{\mathrm{uv}}=\frac{{\sum }_k{u}_k{v}_k{\Delta }_k}{{\sum }_k{u}_k{v}_k},$$

$$E_{pairing}=-\frac{1}{2}Tr(\Delta \kappa )$$

Full text

I 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 Afanasjev and Agbemava (2016). 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. Erler et al. (2012a); Afanasjev et al. (2013); Agbemava et al. (2014)) and more or less consistent picture has been obtained using the combination of different theoretical tools. In addition, systematic theoretical uncertainties Erler et al. (2012a); Agbemava et al. (2014); Afanasjev et al. (2015) and statistical errors Gao et al. (2013); Kortelainen et al. (2013); Agbemava et al. (2019) 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 Dietrich and Pomorski (1998); Dechargé et al. (1999); Bender et al. (2001); Denisov (2005); Gambhir et al. (2015); Ismail et al. (2016), these systematic studies were restricted to spherical symmetry. Our recent study (Ref.  Afanasjev et al. (2018)) 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 studies111The effects of axial and triaxial deformations have also been studied for a few hyperheavy nuclei in Refs. Berger et al. (2001); Warda (2007) and Ref. Staszczak and Wong (2009), respectively. Somewhat larger set of the nuclei was studied with triaxiality included in Ref. Brodziński and Skalski (2013) but according to Ref. Afanasjev et al. (2018) the deformation range employed in this work is not sufficient for $Z\geq 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\sim 130$ region (see also Ref. Brodziński and Skalski (2013) 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. Afanasjev et al. (2018). 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.  Wong (1973)) 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 Warda (2007); Staszczak and Wong (2009, 2014); Ichikawa et al. (2014); Kosior et al. (2017); Nazarewicz et al. (2002). However, in absolute majority of the cases such shapes correspond to highly excited states either at spin zero Staszczak and Wong (2009); Kosior et al. (2017) or at extreme values of angular momentum Staszczak and Wong (2014); Ichikawa et al. (2014); Staszczak et al. (2017). In the former case, such states are unstable against returning to the shape of sphere-like geometry (Ref. Kosior et al. (2017)). 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 Afanasjev et al. (1999).

The present manuscript aims at the extension of the investigations of the properties of hyperheavy nuclei and of nuclear landscape started in Ref. Afanasjev et al. (2018). 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. II. Section III 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. IV. Sec. V 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. VI. Sec. VIII 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. IX. The stability of toroidal shapes in selected nuclei and the evolution of such shapes along their fission path are considered in Sec. X. Sec. XI 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. XII. Finally, Sec. XIII summarizes the results of our work.

II 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. Agbemava et al. (2014)). The stability of prolate minima with $\beta_{2}\sim 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. Agbemava et al. (2016). Triaxial RHB (TRHB code) Agbemava et al. (2017) and triaxial relativistic mean field + BCS (TRMF+BCS code) Abusara et al. (2010) 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 Nikšić et al. (2008). 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 Agbemava et al. (2015, 2014, 2017); Afanasjev and Agbemava (2016); Prassa et al. (2012); Lu et al. (2012); Agbemava et al. (2016). Other functionals such as DD-ME2 Lalazissis et al. (2005), PC-PK1 Zhao et al. (2010) and NL3* Lalazissis et al. (2009)), representing other major classes of covariant density functional models Agbemava et al. (2014), 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. Afanasjev et al. (2018)). 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

[TABLE]

Here $(\lambda,\mu)=(2,0)$, $(\lambda,\mu)=(2,0)$ and $(3,0)$ as well as $(\lambda,\mu)=(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_{\lambda,\mu}>$ stand for the expectation values of the respective multipole moments which are defined as

[TABLE]

$C_{\lambda,\mu}$ in Eq. (1) are corresponding stiffness constants Ring and Schuck (1980) and $q_{\lambda,\mu}$ 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. Bonche et al. (2005). Here the quantity $q_{\lambda,\mu}$ is replaced by the parameter $q_{\lambda,\mu}^{eff}$, which is automatically modified during the iteration in such a way that we obtain $\langle\hat{Q}_{\lambda,\mu}\rangle=q_{\lambda,\mu}$ 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

[TABLE]

on the center-of-mass operator $\hat{Q}_{10}$ to avoid a spurious motion of the center of mass.

The deformation parameters $\beta_{2}$, $\beta_{3}$ and $\gamma$ are extracted from respective multipole moments:

[TABLE]

via

[TABLE]

where $R_{0}=1.2A^{1/3}$. Note that $Q_{22}=0$ and $\gamma=0$ in axially symmetric RHB calculations.

The $\beta_{2}$ and $\gamma$ values have a standard meaning of the deformations of the ellipsoid-like density distributions only for $|\beta_{2}|\lesssim 1.0$ values. At higher $\beta_{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 $\beta_{2}$ values and of necking degree of freedom at large positive $\beta_{2}$ values (see Fig. 2 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<\beta_{2}<3.0$ range are calculated in the axial reflection symmetric RHB framework Agbemava et al. (2014); such large range is needed for a reliable definition of the $\beta_{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 Karatzikos et al. (2010), we use the separable form of the finite-range Gogny pairing interaction introduced in Ref.  Tian et al. (2009). Its matrix elements in $r$-space have the form

[TABLE]

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

[TABLE]

The parameters of this interaction have been derived by a mapping of the 1S0 pairing gap of infinite nuclear matter to that of the Gogny force D1S. The resulting parameters are: $G=728$ MeV fm3 and $a=0.644$ fm Tian et al. (2009). 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 Afanasjev and Abdurazakov (2013); Agbemava et al. (2014); Dobaczewski et al. (2015).

III The effects of the basis truncation

Considering that the results published in Ref. Afanasjev et al. (2018) 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 $\beta_{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 $\beta_{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. 1c and d, respectively. One can see that the deformation of basis $\beta_{0}=0.5\beta_{2}$ typically leads to the lowest in energy solutions. Here $\beta_{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. Afanasjev et al. (2018) and the present manuscript.

Figs. 1a 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 $\beta_{0}=0.5\beta_{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 $\beta_{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 $\beta_{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<\beta_{2}<1.8$ (see Fig. 1b). The deformation ranges $-3.0<\beta_{2}<-1.8$ and $1.8<\beta_{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 $\beta_{2}\sim-4.0$ are seen in high-$Z$/high-$N$ hyperheavy nuclei (see Fig. 1b for the ${}^{466}156$ results and Fig. 1 in Ref. Afanasjev et al. (2018) 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. 2 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 $\beta_{2}=0.0$, the increase of prolate deformation leads to the emergence of hyperdeformed shapes at $\beta_{2}=1.0$, which evolve into the shapes consisting of two fragments connected by neck at higher $\beta_{2}$ values. The separation of the fragments and the size of the neck increases/decreases with increasing $\beta_{2}$ values. Figs. 1b, c and d clearly indicate increased dependence of the results on the parameters $N_{F}$ and $\beta_{0}$ of the basis for these shapes.

The evolution of the shapes for negative $\beta_{2}$ values is shown in Figs. 2a-e. Highly deformed oblate shape exists at $\beta_{2}=-0.5$ deformation which transforms into biconcave disk shape at $\beta_{2}=-1.0$. Further decrease of the $\beta_{2}$ values leads to toroidal shapes. Note that with the increase of absolute value of $\beta_{2}$ the radius of the toroid increases and the tube radius decreases. Total energies and equilibrium deformations of toroidal shapes with $\beta_{2}\sim-1.45$, corresponding to local minimum seen in deformation energy curves of Figs. 1b-d, are rather well described with $N_{F}=20$ and show almost no dependence on the deformation of basis $\beta_{0}$. On the contrary, toroidal shapes with larger (in absolute sense) $\beta_{2}$ values show substantial increase of the dependence of total energies and equilibrium deformations on $N_{F}$ and $\beta_{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 $\beta_{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 $\beta_{2}$ values and two-fragment ones connected by neck.

The equilibrium $\beta_{2}$ values of the $Z=122-138$ nuclei presented in Fig. 2 of Ref. Afanasjev et al. (2018) 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 $\beta_{2}<-1.4$ represent the lowest in energy solutions at axial symmetry in the $Z=140-180$ part of nuclear landscape (see Fig. 3). 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 $\beta_{2}$ values have to be considered as lower limits (in absolute sense). As illustrated in Fig. 1, further increase of the size of fermionic basis will lead to the increase (in absolute sense) of the $\beta_{2}$ values of LEMAS and to more energetically favored status of toroidal shapes as compared with ellipsoidal ones.

IV Density profiles, charge radii and neutron skins of

spherical hyperheavy nuclei

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 and $\alpha$-decay, have been predicted in Ref. Afanasjev et al. (2018). The largest region with the highest fission barriers is centered at $Z\sim 156,N\sim 310$; other two regions are smaller with smaller fission barriers (see Fig. 6 in Ref. Afanasjev et al. (2018)). 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 Afanasjev and Agbemava (2016). 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 ${}^{368}138$, ${}^{466}156$ and ${}^{584}174$ nuclei. As a benchmark, we are using the 208Pb and ${}^{292}120$ nuclei. The properties of latter nucleus were studied in details in Ref. Afanasjev and Frauendorf (2005).

Charge radii $r_{ch}$ and neutron skins $r_{skin}$ of these nuclei are presented in Table 1 and related density distributions are shown in Fig. 4. The predictions of different functionals for charge radii of 208Pb differ by less than 0.01 fm (see Table 1); this is also seen in proton density distributions (see Fig. 4a). 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. 4a). 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 ${}^{292}120$ nucleus (see Table 1 and Fig. 4b).

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 ${}^{368}138$, ${}^{466}156$ and ${}^{584}174$ 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 ${}^{292}120$ nucleus to higher $Z$/higher $N$ nuclei increase also for neutron skins (see Table 1). 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 ${}^{292}120$ 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. Agbemava et al. (2014)). 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. Agbemava et al. (2014) and in Ref. Abrahamyan et al. (2012)). 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\pm 0.17$ Abrahamyan et al. (2012). 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}\sim 0.3$ fm (see Table 1). 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 PRE . Thus, the results obtained in future PREX-2 experiment PRE would be quite useful in helping to discriminate the predictions.

General features of proton and neutron density distributions of hyperheavy nuclei seen in Fig. 4 do not depend on employed functional. Apart of neutron density in ${}^{368}138$ 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 $\omega=\rho^{cent}/\rho^{surf}$ of the density at the center $\rho^{cent}$ to the maximum density at the surface $\rho^{surf}$ averaged over the set of employed functionals to characterize this depression. The central density depressions in neutron subsystems of the ${}^{466}156$ and 584174 nuclei are rather modest with $\omega_{\nu}=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. Dechargé et al. (1999)). However, central depressions are significantly more pronounced in proton subsystems of hyperheavy nuclei with $\omega_{\pi}=0.753$, 0.651 and 0.534 for the ${}^{368}138$, ${}^{466}156$ and 584174 nuclei, respectively. Thus, they are close to the ones expected for semi-bubble nuclei (see Ref. Dechargé et al. (1999)). 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 ${}^{292}120$ superheavy nucleus (with $\omega_{\pi}\sim\omega_{\nu}\sim 0.65$) as compared with higher-$Z$ hyperheavy nuclei (see Fig. 4); the only exception is proton subsystem of the 584174 nucleus. The detailed analysis of Ref.  Afanasjev and Frauendorf (2005) strongly suggests that such central depression in the density distributions of the ${}^{292}120$ 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. Afanasjev and Frauendorf (2005)). 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 Schuetrumpf et al. (2017)). 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.

V 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. Agbemava et al. (2015). Fig. 5 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.  6a). This drop is also visible in calculated $\alpha$-decay half-lives (see Figs. 4 and 5 in supplemental material to Ref. Afanasjev et al. (2018)). However, the impact of the $Z=154$ proton shell closure on two-proton separation energies is substantially smaller (see Fig. 6b) 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 $\alpha$-decay half-lives presented in Fig. 6 of supplemental material to Ref. Afanasjev et al. (2018) 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. Afanasjev et al. (2018) do not reveal the impact of the $N=308$ shell gap. These observations suggest that in the center of the $Z\sim 156,N\sim 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. Afanasjev et al. (2018)] are larger for the DD-PC1 and DD-ME2 CEDFs as compared with the NL3* and PC-PK1 ones.

Figs.  7 and fig.  8 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. 6, 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. 5b), the $N=228$ shell gap between the $2h_{11/2}$ and $2h_{9/2}$ orbitals in the 366138 nucleus (Fig. 7b) and the $N=406$ shell gap between the $2j_{15/2}$ and $2j_{13/2}$ orbitals in the 580174 nucleus (Fig. 8b). 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. Bender et al. (1999)). 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. IV). Note that the gaps obtained in the DD-PC1/DD-ME2 functionals are larger than those seen in the NL3*/PC-PK1 ones.

VI The stability of spherical hyperheavy nuclei with respect of

octupole and triaxial distortions

Fig. 6 in Ref. Afanasjev et al. (2018) 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. 9 shows the potential energy surfaces of these nuclei in the $(\beta_{2},\beta_{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 $\beta_{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 $\beta_{2}=-0.196,\beta_{3}=0.0$ (as obtained in axial RHB calculations without octupole deformation) to $\beta_{2}=-0.198,\beta_{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.  9 do not suggest that cluster radioactivity from spherical 368138 and ${}^{466}156$ hyperheavy nuclei plays an important role. Super-asymmetric fission channel leading to a cluster radioactivity shows itself as a narrow fission path in the $(\beta_{2},\beta_{3})$ plane separate from main fission path (see example in Fig. 1 of Ref. Warda et al. (2018)). No such path is visible in Figs. 9a and b. However, because of convergence problems one cannot define whether cluster radioactivity is important in the 584174 nucleus.

Fig. 10 shows the potential energy surfaces obtained in the TRHB calculations. In the 584174 nucleus, the axial saddles are located at $\beta_{2}=0.17,\gamma=60^{\circ}$ (the barrier on the oblate side of spherical minimum) and at $\beta_{2}=0.197,\gamma=0^{\circ}$ (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 $\beta_{2}=0.16,\gamma=36^{\circ}$ 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 $\beta_{2}$. The saddles of axial fission barriers (on oblate and prolate sides of spherical minimum) are located at modest quadrupole deformation of $\beta_{2}\sim 0.2$. As a result, the distance between these two saddles in the $(\beta_{2},\gamma)$ plane is relatively small, so that large changes in binding energy due to triaxiality for nearly constant $\beta_{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.

VII Pairing interaction in hyperheavy nuclei

The magnitude of pairing interaction could be accessed via two calculated quantities: average pairing gap and pairing energy Ring and Schuck (1980). There are several definitions of average pairing gaps in literature (see discussion in Section IV of Ref. Agbemava et al. (2014)). The pairing gap

[TABLE]

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. Agbemava et al. (2014) showed that the $\Delta_{\rm 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

[TABLE]

where $\Delta$ and $\kappa$ are pairing field and pairing tensor, respectively Ring and Schuck (1980). 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. Afanasjev et al. (2000); Valor et al. (2000)).

Fig. 13 shows the evolution of neutron and proton pairing energies $E_{pairing}$ and pairing gaps $\Delta_{uv}$ as a function of $\beta_{2}$ for the lowest in energy solution in the 466156 nucleus. The oscillating behavior of these quantities as a function of $\beta_{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. Karatzikos et al. (2010)). 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 $\Delta_{uv}$ corresponds to low/high density of the single-particle states in the vicinity of the Fermi level (see Ref. Karatzikos et al. (2010)). Indeed, in the nucleus under consideration there is a substantial reduction of neutron $|E_{pairing}|$ and $\Delta_{uv}$ values near spherical shape which is attributable to the presence of large $N=308$ shell gap (see Fig. 5). Note also that there is a collapse of neutron pairing correlations at toroidal shapes in the $\beta_{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. 1d). The local minima in potential energy surfaces correspond to low density of the single-particle states in the vicinity of the Fermi level Strutinsky (1967, 1968); 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\sim 156,N\sim 310)$ island of stability of spherical hyperheavy nuclei are shown in Fig. 14. 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. Anguiano et al. (2002)) will restore the pairing correlations in these nuclei. However, in no way it will affect the conclusions of the present work and of Ref. Afanasjev et al. (2018). 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. 14b; 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\geq 168$ nuclei the $|E_{pairing}|$ values exceed 20 MeV.

Figure 14 shows the evolution of proton and neutron pairing energies as a function of the $\beta_{2}$ and $\gamma$ 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.

VIII 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 Bürvenich et al. (2004); Abusara et al. (2012). This is illustrated in Fig. 11 which shows the competition in energy of two local minima, namely, prolate superdeformed minimum with $\beta_{2}\sim 0.6$ and oblate one with $\beta_{2}\sim-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. Abusara et al. (2012)). 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. Afanasjev et al. (2018)).

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. 11), we have performed the calculations in axial octupole deformed code RHB-OCT developed in Ref. Agbemava et al. (2016). The results of such calculations are illustrated in Fig. 12. One can see that the inclusion of octupole deformation leads to the instability of the prolate superdeformed $\beta_{2}\sim 0.6$ minima with respect of octupole deformation. In the ${}^{316,320}124$ nuclei this instability will lead to asymmetric fission, while the competition of symmetric and asymmetric fissions is possible in the ${}^{296,300}122$ 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 ${}^{296}122$ and ${}^{316,320}124$ nuclei in Fig. 12) 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 ${}^{300}122$ nucleus in Fig. 12). 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. Abusara et al. (2012)). 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. Afanasjev et al. (2018)).

IX 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. 16. The minimum of deformation energy curve is located at $\beta_{2}\sim-4.0$ for proton-rich nuclei with $N=186-190$. This minimum appear at $\beta_{2}\sim-2.5$ for the $N=194-230$ nuclei. Further increase of neutron number leads to the growth (in absolute sense) of the $\beta_{2}$ values: the nuclei with $N=246-262$ have minima at $\beta_{2}\sim-3.0$. The minima of the deformation energy curves are located at $\beta_{2}\sim-3.5$ for neutron-rich nuclei with $N=266-326$. As discussed in details in Ref. Afanasjev et al. (2018), the nuclei have toroidal shapes in such minima. The competing local minima with $\beta_{2}\sim-0.8$ (which corresponds to oblate ellipsoidal shape, see Ref. Afanasjev et al. (2018)) 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\sim 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. Afanasjev et al. (2018), these minima in deformation energy curves with toroidal shapes are potentially unstable with respect of the transition to prolate shape via $\gamma$-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 $\beta_{2}$ part of deformation energy curve onto negative $\beta_{2}$ values] with solid ones in Figs. 1c and d of Ref. Afanasjev et al. (2018)). However, as discussed in Ref. Afanasjev et al. (2018) and in Sec. X below, some of these minima could be stable.

Proton and neutron chemical potentials for the solutions, displayed in Fig. 16, are shown in Fig. 17. They behave differently as a function of $\beta_{2}$. Neutron chemical potential on average is almost flat as a function of $\beta_{2}$. The magnitude of the fluctuations of the $\lambda_{n}(\beta_{2})$ values with respect of average values decreases on approaching two-neutron drip line. On the contrary, there is a pronounced slope in the $\lambda_{p}(\beta_{2})$ values: they on average decrease with increasing absolute value of $\beta_{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 $\beta_{2}$ values in which the proton chemical potential is positive (see top row of Fig. 17) and the shape of the nucleus in respective local minimum is oblate ellipsoidal. Even-even nuclei with $\lambda_{p}>0$ are typically expected to be unstable with respect of proton emission222 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 $\lambda_{p}=dE/dZ$ as a point (nucleus) of the transition from negative $\lambda_{p}$ (“bound” nuclei) to positive $\lambda_{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 Afanasjev et al. (2015). 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 $\lambda_{n}\sim 0.05$ MeV. These results were obtained in Ref. Afanasjev et al. (2015) from the calculations of Refs. Afanasjev et al. (2013); Agbemava et al. (2014) 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 $\beta_{2}\sim-1.5$ and it becomes more negative with the increase of absolute value of $\beta_{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 $\beta_{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. 24 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 $\beta_{2}$.

X Toroidal shapes: stability and shape evolution

along the fission path.

The investigations of Ref. Afanasjev et al. (2018) showed that some toroidal shapes could be stable with respect of triaxial distortions. Fig. 18 shows potential energy surfaces of the 354134 and 348138 nuclei around minima of such configurations located at $\beta_{2}\sim 2.3,\beta_{4}\sim+1.5,\gamma=60^{\circ}$. 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. Afanasjev et al. (2018).

To understand these reasons the evolution of toroidal shapes along the fission path of the configuration in the 354134 nucleus (shown in Fig. 18) is displayed in Fig. 19. The toroid and its tube are fully symmetric at the minimum (Fig. 19a-c). The deviations from axial symmetry lead to the distortions which are already seen at $\gamma=50^{\circ}$; the toroid is stretched out in the direction of the axis of symmetry and squeezed in perpendicular direction (Fig. 19f)). However, the tubes of the toroid still remain cylindrical (Fig. 19d,e)). Figs. 19g,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 $\gamma$- and $\beta_{2}$ deformations leads to increasing distortion of the geometry of toroid ((Fig. 19l)) the outer shape of which looks now similar to parallelogram and the shape of toroid hole comes closer to square. In addition, Figs. 19j,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.

XI 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. Afanasjev et al. (2018) and Sec. VIII 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. 20 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 $\gamma=0^{\circ}$ 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. 20) as compared with superheavy nuclei and the mountain centered around $\beta_{2}\sim 0$ is formed in this canyon. The slope of the mountain in the direction of the $\beta_{2}$-deformation at $\gamma=0^{\circ}$ 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. 20) or at very large $\gamma$-deformation (see example of the 432134 nucleus in Fig.  20). In addition, the fission paths from these minima proceed at larger $\gamma$-deformations as compared with superheavy nuclei. Not only the fission through the $\gamma$-plane gets more energetically favored, but also the fission path through $\gamma$-plane becomes much shorter than the one through the $\gamma=0^{\circ}$ axis (see also the discussion in Ref. Afanasjev et al. (2018)).

The general conclusion is that the barriers along the fission paths emerging from the oblate minima located within the $-1.0<\beta_{2}\leq 0.0$ range decrease with increasing proton number (see Ref. Afanasjev et al. (2018)). As a result (see discussion in Sect. XII 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\leq 120$ and $N\leq 184$ Abusara et al. (2010); Warda et al. (2002); Staszczak et al. (2009); Möller et al. (2009); Abusara et al. (2012); Schunck and Robledo (2016) . 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 Afanasjev et al. (2018).

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+\Delta N$, where $\Delta 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 2. Together with the results of the TRHB calculations presented in Refs. Agbemava et al. (2017); Afanasjev et al. (2018) 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 landscape333The 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. Afanasjev et al. (2018) and the discussion in Sec. XII). This fact has been ignored in many studies of nuclear landscape in neutron-rich actinides and superheavy nuclei (such as Refs.  Erler et al. (2012a); Agbemava et al. (2014)) 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. XII 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 Arnould et al. (2007); Martínez-Pinedo et al. (2007). However, so far the fission barriers for such nuclei were calculated only in theoretical frameworks restricted to axial symmetry Goriely et al. (2007); Erler et al. (2012b); Giuliani et al. (2018).

The results for inner fission barrier heights $E^{B}_{triax}$ and the decrease of the height of inner fission barrier due to triaxiality $\Delta E^{gain}$, presented in Table 2, are summarized in Fig. 21. The $E^{B}_{triax}$ values show oscillatory behavior as a function of neutron number $N$ with maxima seen at $N\sim 162$ and $N\sim 222$ (at $N=172$ and 232 in the $Z=118$ nuclei) and minima at $N=192$ and $N\sim 252$ (at $N=202$ and 272 in the $Z=118$ nuclei). More erratic behavior is seen for the $\Delta E^{gain}$ values; note that the triaxiality has no effect on the heights of inner fission barriers for $\Delta 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 $\Delta E^{gain}$ curves are the consequences of underlying shell structure and its evolution with proton and neutron numbers (see discussion in Ref. Abusara et al. (2010)). A general trend of the lowering of the height of inner fission barrier with the increase of proton number is clearly seen in Table 2 and Fig. 21a.

Table 2 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\sim 114$, $N\sim 192$) and ($Z\sim 118$, $N\sim 202$) which are shown in Fig. 22b. 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. Afanasjev et al. (2018)). These factors together lead to a substantial reduction of the region of potentially stable elliposoidal-like nuclei in the $N\geq 258$ region (compare panels (a) and (b) in Fig. 22 below).

XII 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. Möller et al. (1995); Delaroche et al. (2010); Erler et al. (2012a, 2013); Afanasjev et al. (2013); Agbemava et al. (2014) and, in particular, the compilation presented in Sec. VIII of Ref. Agbemava et al. (2014)). The studies of Refs. Erler et al. (2012a); Afanasjev et al. (2013); Agbemava et al. (2014) 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. Afanasjev et al. (2018). 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. Afanasjev et al. (2018) 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 $\gamma$-deformation (see Ref. Möller et al. (2008) for the results obtained in microscopic+macroscopic method) and octupole deformation shows up in the ground states of the nuclei in few localized regions Möller et al. (2008); Agbemava et al. (2016); Agbemava and Afanasjev (2017).

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. Afanasjev et al. (2018) and Table 2 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. Abusara et al. (2012)). 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. Afanasjev et al. (2018)), with increasing proton number beyond $Z=140$, the LEMAS always have toroidal shapes (see discussion in Sec. III and in Ref. Afanasjev et al. (2018)). 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.  24).

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. 23 shows the calculated Coulomb energies $E_{Coul}$ as a function of the $\beta_{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 $\beta_{2}$ this effect is especially pronounced for negative $\beta_{2}$ values. Moreover, the magnitude of $E_{Coul}$ and its slope with deformation for negative $\beta_{2}$ values drastically increases with increasing proton number (see Fig. 23). This is also quantified in Fig. 23 by the $\Delta E_{Coul}$ quantity which is the reduction of the Coulomb energy due to the transition from spherical shape to typical toroidal one. The $\Delta 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. Afanasjev et al. (2018)).

It is clear that nuclear landscape shown in Fig. 24 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. Afanasjev et al. (2018) 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. X. 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. III) 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. Wong (1973) 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 multifragmentation444There are also some experimental indications of the role of multifragmentation in toroidal nuclei, but they are restricted to a single 86Kr + 93Nb reaction Stone et al. (1997).. 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 ${}^{354}134$ and ${}^{348}138$ nuclei in Ref. Afanasjev et al. (2018) (see also Sec. X 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. Afanasjev et al. (2018)). These islands centered around ($Z\sim 138,N\sim 230$), ($Z\sim 156,N\sim 310$) and ($Z\sim 174,N\sim 410$) are shown in Fig. 24.

The analysis of Figs. 22 and 24 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 Dobaczewski et al. (2014) which are especially large for the position of two-neutron drip line Erler et al. (2012a); Afanasjev et al. (2013); Agbemava et al. (2014); Afanasjev et al. (2015); Afanasjev and Agbemava (2016) and fission barriers Agbemava et al. (2017); Jachimowicz et al. (2017). In the CDFT framework, systematic theoretical uncertainties due to the form of the CEDF are substantially larger than statistical errors Agbemava et al. (2019). 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. Afanasjev et al. (2018)). 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.

XIII 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\sim 138,N\sim 230$), ($Z\sim 156,N\sim 310$) and ($Z\sim 174,N\sim 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 PRE 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\sim 156,N\sim 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 $\beta_{2}\sim 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\sim 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 $\beta_{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\sim 138,N\sim 230$), ($Z\sim 156,N\sim 310$) and ($Z\sim 174,N\sim 410$) islands of stability of spherical hyperheavy nuclei have neutron to proton ratios of $N/Z\sim 1.67$, $N/Z\sim 1.99$ and $N/Z\sim 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 Abbott et al. (2017). In a similar fashion, the regions of neutron stars with nuclear pasta phases Caplan et al. (2015); Fattoyev et al. (2017); Kycia et al. (2017) 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\sim 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. Stone et al. (1997).

Any extrapolation beyond known regions of nuclear chart in which the functionals have been fitted is associated with theoretical uncertainties Erler et al. (2012a); Agbemava et al. (2014); Afanasjev et al. (2015). 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. Warda (2007) 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 Agbemava et al. (2017). 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. XI in the present paper and Refs. Brodziński and Skalski (2013); Afanasjev et al. (2018)). 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. Afanasjev et al. (2018)) and on the details of underlying shell structure (see Sect. V 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 Afanasjev et al. (2015). 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.

XIV 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.