Examination of Possible Proton Magic Number Z = 126 with the Deformed Relativistic Hartree-Bogoliubov Theory in Continuum

Abstract

Whether $Z=126$ is a proton magic number has been controversial in nuclear physics. The even-even ${}_{126}\mathrm{Ubh}$ isotopes are calculated based on the DRHBc calculations with PC-PK1. The evolutions of quadrupole deformation and pairing energies for neutron and proton are analyzed to study the possible nuclear magicity. Spherical shape occurs and neutron pairing energy vanishes at $N=258$ and 350, which are the results of possible neutron magicity, while the proton pairing energy never vanishes in Ubh isotopes, which does not support the proton magicity at $Z=126$. In the single-proton spectrum, there is no discernible gap at $Z=126$, while significant gaps appear at $Z=120$ and 138. Therefore, $Z=126$ is not supported as a proton magic number, while $Z=120$ and 138 are suggested as candidates of proton magic numbers.

1. Introduction

The exploration of the limit of nuclear existence has been a very fascinating topic in nuclear physics [1,2,3,4]. With neutron number N and proton number Z as horizontal and vertical axes, respectively, the southwest coast of the nuclear landscape corresponds the beginning of periodic table, H, the lightest element. The northwest and southeast coasts are the proton and neutron drip lines, respectively. The proton drip line has been experimentally determined up to Np with $Z=93$ [5], while the neutron drip line has only been determined up to Ne with $Z=10$ [6]. Near the drip lines, many interesting phenomena are discovered, including the nuclear halo [7], changes of nuclear magic numbers [8] and pygmy resonances [9], and have attracted worldwide attentions. The northeast coast of the nuclear landscape corresponds to the upper limit of mass number A, and is still unknown because data are extremely limited. The heaviest nuclides discovered so far are ${}^{294}\mathrm{Og}$ and ${}^{294}\mathrm{Ts}$ [10], with nucleon numbers $(Z,N)=(118,176)$ and $(117,177)$, respectively. These two nuclides have not reached the so-called “island of stability”, which is theoretically predicted based on the possible neutron magic number $N=184$ and proton magic number $Z=114$ [11,12,13,14].

The concept of nuclear magic numbers refers to the particular stability of nuclei with certain nucleon numbers. Such stability is reflected on the extra nuclear binding energy, and a more intuitive manifestation is a significantly larger separation energy than that of the next nuclide. The nuclear magicity is closely related to the nuclear shell effect, which is a hallmark characteristic in the atomic nucleus as a quantum system. The experimentally confirmed magic numbers are the neutron numbers $N=2,8,20,28,50,82,126$ and proton numbers $Z=2,8,20,28,50,82$ [15]. These magic numbers can be perfectly reproduced by assuming a harmonic oscillator potential plus spin-orbital coupling [16,17]. In a naïve consideration with this model, extra stability is suggested at nucleon numbers $2,8,20,28,50,82,126,184,258,350$, …, which not only explains all the magic numbers above, but also provides a reference for new possible magic numbers such as 184 and 258, as well as $Z=126$. The traditional concept of magic number has been implicitly associated with spherical symmetry. Recent development has extended this concept to exotic shape symmetries and fourfold octupole magic number N = 136 and 196 have been predicted with symmetry-induced increase in stability [18,19]. The nuclei with exotic shape are predicted to form islands in the superheavy region [20].

The possibility of a proton magic number at $Z=126$ has been discussed since as early as 1955 [21]. Following the establishment of the nuclear shell model, the potential existence of stable nuclides near $Z=126$ was studied and there were also controversies in the conclusions [22,23,24,25]. In addition to $Z=126$, subsequent efforts have also been made to study other possible magic numbers in the superheavy nuclear region. Based on phenomenological models such as finite-range droplet model (FRDM), $Z=114$ and $N=184$ are predicted to have large shell gaps, corresponding to new magic numbers [26]. In the FRDM, other large shell gaps are predicted to be located at $N=162$ and $Z=104,106,108,110$, instead of $Z=126$. In order to predict experimentally unknown areas, the prediction stability is a crucial issue to address. The predictive stability of micro-macro models has been investigated thoroughly in many literatures, such as Refs. [27,28] which use the inverse problem theory of applied mathematics and Monte Carlo simulations. For microscopic models, the theoretical uncertainties have also been investigated in a systematic manner [29,30].

The relativistic density functional theory has been proven to be a powerful tool in nuclear physics, due to its successful descriptions on many nuclear phenomena [4]. For the exotic nuclei far away from stability line, the occupation of single-nucleon spectrum is very close to the continuum threshold, and the pairing interaction can scatter nucleons from bound states to resonant states in the continuum, leading to a more diffuse density and the dripline locations might be influenced, called the continuum effects [31,32,33,34]. A proper treatment for continuum is solving the nucleon system described by Bogoliubov transformation in coordinate space, where wave functions are approximated on a spatial lattice, and the continuum is discretized by suitably large box boundary conditions [35,36,37,38]. Based on the relativistic density functional theory and solving the problem in coordinate space, the spherical relativistic continuum Hartree-Bogoliubov (RCHB) theory was developed [38,39], which can properly take into account the effects of pairing correlations and continuum for the nuclei near the limit of nuclear landscape. The RCHB theory has been applied to many studies for both stable and exotic nuclei, including describing the halo in ${}^{11}\mathrm{Li}$ [39], predicting giant halos [40], interpreting the pseudospin symmetry in exotic nuclei [41], reproducing the interaction cross section and charge-changing cross sections in light nuclei [42,43], etc. Based on the RCHB theory, the shell structures for superheavy nuclei were studied, and $Z=120,132,138$, and $N=172,184,198,228,238,258$ were suggested to be the magic numbers [44]. In Ref. [45], the first nuclear mass table for the nuclei with $8⩽Z⩽120$ that incorporates continuum effects was constructed based on the RCHB theory, where the evolution of shell structures and magic numbers were studied.

Considering that most nuclei deviate from the spherical shape, the deformed relativistic Hartree-Bogoliubov theory in continuum (DRHBc) was developed [46,47]. The DRHBc theory takes into account the effects of deformation, pairing correlation and continuum, and can provide proper descriptions for both the stable nuclei and the unstable exotic nuclei near the boundary of nuclear landscape. As the advantages of the RCHB theory are inherited and the deformation degrees of freedom are further included, the DRHBc theory has been applied in many studies on exotic nuclei, including the halo structures [48,49,50,51,52,53,54,55,56], dripline locations [57,58,59], evolution of deformation and shape coexistence [60,61,62,63], etc. Recently, a nuclear mass table for the nuclei with $8⩽Z⩽120$ is in progress [64,65], and the even-Z part has been established [66,67]. The mass description for superheavy nuclei reach accuracy of several hundred keV [57] in reproducing available experimental data as well as empirical data in AME2020 [68]. Based on the DRHBc theory, $Z=120$ is suggested as a candidate of proton magic number [66], while $Z=126$ has not been considered yet.

In this work, the DRHBc theory is employed to examine the possible proton magic number $Z=126$. The $Z=126$ element has a temporary systematic IUPAC name as Unbiohexium (Ubh) [69]. The theoretical framework is briefly introduced in Section 2. The numerical details are given in Section 3. The results and discussions are presented in Section 4. Finally, a summary is given in Section 5.

2. Theoretical Framework

The starting point of the relativistic density functional theory is an effective Lagrangian density with either the meson-exchange or point-coupling interactions [4]. For the point-coupling interaction, the Lagrangian density reads

$$\begin{array}{cc}\hfill \mathcal{L}=& \overline{\psi }(i{\gamma }_{\mu }{\partial }^{\mu }-M)\psi \hfill \\ \hfill & -{\displaystyle \frac{1}{2}}{\alpha }_S\left(\overline{\psi }\psi \right)\left(\overline{\psi }\psi \right)-{\displaystyle \frac{1}{2}}{\alpha }_V\left(\overline{\psi }{\gamma }_{\mu }\psi \right)\left(\overline{\psi }{\gamma }^{\mu }\psi \right)-{\displaystyle \frac{1}{2}}{\alpha }_{TV}\left(\overline{\psi }\overrightarrow{\tau }{\gamma }_{\mu }\psi \right)\left(\overline{\psi }\overrightarrow{\tau }{\gamma }^{\mu }\psi \right)\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}{\alpha }_{TS}\left(\overline{\psi }\overrightarrow{\tau }\psi \right)\left(\overline{\psi }\overrightarrow{\tau }\psi \right)-{\displaystyle \frac{1}{3}}{\beta }_S{\left(\overline{\psi }\psi \right)}^3-{\displaystyle \frac{1}{4}}{\gamma }_S{\left(\overline{\psi }\psi \right)}^4-{\displaystyle \frac{1}{4}}{\gamma }_V{\left[\left(\overline{\psi }{\gamma }_{\mu }\psi \right)\left(\overline{\psi }{\gamma }^{\mu }\psi \right)\right]}^2\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}{\delta }_S{\partial }_{\nu }\left(\overline{\psi }\psi \right){\partial }^{\nu }\left(\overline{\psi }\psi \right)-{\displaystyle \frac{1}{2}}{\delta }_V{\partial }_{\nu }\left(\overline{\psi }{\gamma }_{\mu }\psi \right){\partial }^{\nu }\left(\overline{\psi }{\gamma }^{\mu }\psi \right)-{\displaystyle \frac{1}{2}}{\delta }_{TV}{\partial }_{\nu }\left(\overline{\psi }\overrightarrow{\tau }{\gamma }_{\mu }\psi \right){\partial }^{\nu }\left(\overline{\psi }\overrightarrow{\tau }{\gamma }^{\mu }\psi \right)\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}{\delta }_{TS}{\partial }_{\nu }\left(\overline{\psi }\overrightarrow{\tau }\psi \right){\partial }^{\nu }\left(\overline{\psi }\overrightarrow{\tau }\psi \right)-{\displaystyle \frac{1}{4}}{F}^{\mu \nu }{F}_{\mu \nu }-e\overline{\psi }{\gamma }^{\mu }{\displaystyle \frac{1-{\tau }_3}{2}}{A}_{\mu }\psi ,\hfill \end{array}$$

(1)

where M is the nucleon mass, e is the charge unit, and $\psi$ is the field operator for nucleon. The constants $\alpha ,\beta ,\gamma ,\delta$ terms correspond to four-fermion, third-order, fourth-order, and derivative couplings terms, respectively [70,71]. The subscripts S, V, and T mean scalar, vector, and isovector, respectively. The last two terms describe the electromagnetic interaction in the standard quantum electrodynamics Lagrangian, and $A_{\mu }$ and $F_{\mu \nu }$ are the four-vector potential and field strength tensor of the electromagnetic field, respectively.

Starting from Equation (1) and performing Legendre transformation, the Hamiltonian is obtained, which, by using the mean-field approximation, is then utilized to derive the energy density functional. By simultaneously treating the mean fields and pairing correlations in a self-consistent manner, the nucleons are described by the the relativistic Hartree-Bogoliubov (RHB) equation [72],

$$\left(\begin{array}{cc}{\widehat{h}}_D-{\lambda }_{\tau }& \widehat{\Delta }\\ -{\widehat{\Delta }}^{*}& -{\widehat{h}}_D^{*}+{\lambda }_{\tau }\end{array}\right)\left(\begin{array}{c}{U}_k\\ V_k\end{array}\right)=E_k\left(\begin{array}{c}{U}_k\\ V_k\end{array}\right),$$

(2)

where ${\lambda }_{\tau }$ is the Fermi energy of neutron or proton $(\tau =n,p)$, $E_k$ is the quasiparticle energy, $U_k$ and $V_k$ are the quasiparticle wave functions, ${\widehat{h}}_D$ is the Dirac Hamiltonian, and $\widehat{\Delta }$ is the pairing potential.

In the coordinate space, the Dirac Hamiltonian can be written as

$$h_D\left(\mathit{r}\right)=\mathit{\alpha }·\mathit{p}+V\left(\mathit{r}\right)+\beta [M+S\left(\mathit{r}\right)],$$

(3)

where $S\left(\mathit{r}\right)$ and $V\left(\mathit{r}\right)$ are scalar and vector potentials, respectively,

$$\begin{array}{cc}\hfill S\left(\mathit{r}\right)& ={\alpha }_S{\rho }_S+{\beta }_S{\rho }_S^2+{\gamma }_S{\rho }_S^3+{\delta }_S▵{\rho }_S,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill V\left(\mathit{r}\right)& ={\alpha }_V{\rho }_V+{\gamma }_V{\rho }_V^3+{\delta }_V▵{\rho }_V+e{A}^0+{\alpha }_{TV}{\tau }_3{\rho }_3+{\delta }_{TV}{\tau }_3▵{\rho }_3.\hfill \end{array}$$

(5)

${\rho }_S$, ${\rho }_V$, and ${\rho }_3$ are the local densities, and are defined as

$$\begin{array}{cc}\hfill {\rho }_S\left(\mathit{r}\right)& =\sum _{k>0}{V}_k^{†}\left(\mathit{r}\right){\gamma }_0{V}_k\left(\mathit{r}\right),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\rho }_V\left(\mathit{r}\right)& =\sum _{k>0}{V}_k^{†}\left(\mathit{r}\right)V_k\left(\mathit{r}\right),\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\rho }_3\left(\mathit{r}\right)& =\sum _{k>0}{V}_k^{†}\left(\mathit{r}\right){\tau }_3{V}_k\left(\mathit{r}\right).\hfill \end{array}$$

(8)

$k>0$ means that the summation runs over the quasiparticle states in the Fermi sea, which corresponds to the so-called no-sea approximation.

The pairing potential $\widehat{\Delta }$ reads

$$\Delta ({\mathit{r}}_1,{\mathit{r}}_2)=V^{pp}({\mathit{r}}_1,{\mathit{r}}_2)\kappa ({\mathit{r}}_1,{\mathit{r}}_2),$$

(9)

where $V^{pp}$ is the pairing force, and $\kappa$ is the pairing tensor. In this work, the density-dependent zero-range pairing force

$$V^{pp}({\mathbf{1}}_1,{\mathit{r}}_2)=V_0{\displaystyle \frac{1}{2}}(1-P^{\sigma })\delta ({\mathit{r}}_1-{\mathit{r}}_2)\left(1-{\displaystyle \frac{\rho \left({\mathit{r}}_1\right)}{{\rho }_{\mathrm{sat}}}}\right)$$

(10)

is adopted. The pairing tensor is defined by using quasiparticle wavefunctions as [73]

$$\kappa =U\ast V^T.$$

(11)

In the DRHBc theory, since axial deformation and spatial reflection symmetry are assumed, the potentials and densities are expanded in terms of the Legendre polynomials,

$$f\left(\mathit{r}\right)=\sum _{\lambda }{f}_{\lambda }\left(r\right)P_{\lambda }(cos\theta ),\lambda =0,2,4,\dots ,{\lambda }_{max}.$$

(12)

It is worth mentioning that this work focuses on axial and spatial reflection symmetries, while higher-order deformation, spatial-reflection asymmetry as well as triaxial shapes, which might also play important roles in heavy nuclei [19,20,28,30], are not included here. The conclusion is limited due to such an assumption on symmetry, and future works are expected by taking into account the full degrees of freedom that nuclei can take.

In order to correctly take into account the continuum effect, the deformed RHB Equation (2) are solved in a spherical Dirac Woods-Saxon basis [74], which can properly describe the asymptotic behavior of the density distribution at a large r.

After self-consistently solving the RHB equations, the expectation values can be calculated. The total energy is calculated as [47,64]

$$\begin{array}{cc}\hfill E_{\mathrm{tot}}=& \sum _{k>0}({\lambda }_{\tau }-E_k)v_k^2-E_{\mathrm{pair}}\hfill \\ \hfill & -\int {\mathrm{d}}^3\mathit{r}\left({\displaystyle \frac{1}{2}}{\alpha }_S{\rho }_S^2+{\displaystyle \frac{1}{2}}{\alpha }_V{\rho }_V^2+{\displaystyle \frac{1}{2}}{\alpha }_{TV}{\rho }_3^2\right.+{\displaystyle \frac{2}{3}}{\beta }_S{\rho }_S^3+{\displaystyle \frac{3}{4}}{\gamma }_S{\rho }_S^4+{\displaystyle \frac{3}{4}}{\gamma }_V{\rho }_V^4\hfill \\ \hfill & +\left.{\displaystyle \frac{1}{2}}{\delta }_S{\rho }_S▵{\rho }_S+{\displaystyle \frac{1}{2}}{\delta }_V{\rho }_V▵{\rho }_V+{\displaystyle \frac{1}{2}}{\delta }_{TV}{\rho }_3▵{\rho }_3+{\displaystyle \frac{1}{2}}{\rho }_{p}e{A}^0\right)+E_{\mathrm{c}.\mathrm{m}.},\hfill \end{array}$$

(13)

where

$$v_k^2=\int d^3\mathit{r}{V}_k^{†}\left(\mathit{r}\right)V_k\left(\mathit{r}\right),$$

(14)

and $E_{\mathrm{pair}}$ and $E_{\mathrm{c}.\mathrm{m}.}$ are the pairing energy and center-of-mass correction energy [47], respectively. It should be mentioned that the pairing energy $E_{\mathrm{pair}}$ is a measure for the size of pairing correlations in theoretical calculations, and under the zero-range pairing force it is calculated as

$$E_{\mathrm{pair}}=-{\displaystyle \frac{1}{2}}\int d^3\mathit{r}\kappa \left(\mathit{r}\right)\Delta \left(\mathit{r}\right).$$

(15)

The quadrupole deformation is calculated using

$${\beta }_2={\displaystyle \frac{\sqrt{5\pi }{Q}_2}{3N\langle r^2\rangle }}={\displaystyle \frac{2\sqrt{5\pi }\int {\mathrm{d}}^3\mathit{r}{\rho }_v\left(\mathit{r}\right)r^2{P}_2(cos\theta )}{3\left[\int {\mathrm{d}}^3\mathit{r}{\rho }_v\left(\mathit{r}\right)\right]\left[\int {\mathrm{d}}^3\mathit{r}{\rho }_v\left(\mathit{r}\right)r^2\right]}},$$

(16)

where N is the corresponding particle number, $\langle r^2\rangle$ is the root-mean-square (rms) radius, $Q_2$ is the intrinsic quadrupole moment, and ${\rho }_v$ is the vector density distribution.

3. Numerical Details

In this work, the relativistic density functional PC-PK1 [71] is adopted, which leads to the root-mean-square deviation of 1.14 MeV in reproducing available mass data [75,76], among one of highest-accuracy density-functional descriptions for mass. The accuracy for mass descriptions by PC-PK1 for medium-heavy to superheavy regions reach several hundred keV, as illustrated in Refs. [57,77]. In the coordinate space, the box size $R_{\mathrm{box}}=20$ fm and the mesh size $\Delta r=0.1$ fm. The cutoff parameters for the Dirac Woods-Saxon basis include the angular momentum cutoff $J_{max}=23/2\hslash$ and the energy cutoff for positive-energy states $E_{\mathrm{cut}}^{+}=300$ MeV, and the number of negative-energy states is taken the same as that of positive-energy states [74]. For the Legendre expansion in Equation (12), a cutoff ${\lambda }_{max}=10$ is taken. In the pairing channel, the pairing strength for Equation (10), $V_0=-325\mathrm{MeV}{\mathrm{fm}}^3$, and a pairing window of 100 MeV are taken. These pairing parameters are determined by reproducing the odd-even mass differences for Ca and Pb isotopes [64]. The examinations for the above numerical cutoffs and the pairing parameters can be found in Refs. [64,66].

4. Results and Discussions

The even-even nuclei in the ${}_{126}\mathrm{Ubh}$ isotopic chain are calculated with the DRHBc theory. Bulk properties for the ground states, such as binding energy, separation energies, quadrupole deformation, root-mean-square radii and pairing energies, as well as the corresponding density distributions and single-particle levels, are obtained. In order to discuss the possible magicity, in the following we will focus on the evolution of two-proton separation energy and two-proton gap, deformation and pairing energies, as well as single-particle levels.

The separation energy belongs to the first rang observables for magic numbers, and corresponds to the characteristic jumps in separation energies [78]. As has been shown in Figures 2 and 3 of Ref. [66], in the DRHBc calculations for $8⩽Z⩽120$, abrupt changes of separation energies are exhibited at magic numbers. In order to explore the possible magicity for superheavy nuclei, Figure 1a shows the two-proton separation energy $S_{2p}$ in the DRHBc calculations by taking the even-even $N=258$ isotopes with $Z>100$ as examples. It is shown that $S_{2p}$ decreases monotonously with the proton number, with the slopes larger near $Z=120$ and smaller near $Z=130$. For a more intuitive view, Figure 1b shows the corresponding decrease of $S_{2p}$ at each Z, i.e., the two-proton gap, ${\delta }_{2p}(Z,N)=S_{2p}(Z,N)-S_{2p}(Z+2,N)$. The drastic changes of the two-nucleon separation energies are implicated by the peaks of the two-nucleon gaps, as has been shown in Figures 4 and 5 of Ref. [66]. Here a peak at $Z=120$ is noted, which is consistent with several other density functional predictions that $Z=120$ is a proton magic number. On the contrary, the ${\delta }_{2p}$ of $Z=126$ corresponds to a decrease, instead of a peak, which has also been depicted in the slope of $S_{2p}$. Therefore, from the two-proton separation energy, $Z=126$ does not possess extra stability, and no characteristic of magicity is shown.

When the nucleon number equals to a traditional magic number, one result of RHB calculation is that at the lowest minimum of the potential energy curve, i.e., the ground state, usually has a spherical shape [66]. In order to further discuss the possibility of the proton magicity at $Z=126$, Figure 2 shows the quadrupole deformation ${\beta }_2$ as a function of the neutron number in the DRHBc calculations for even-even ${}_{126}\mathrm{Ubh}$ isotopes. From the proton drip line and by increasing the neutron number, the deformation ${\beta }_2$ first generally decreases from a large positive value to zero near $N=258$. After that, ${\beta }_2$ increases to a large value again and decreases to zero near $N=350$. The ${\beta }_2$ at $324⩽N⩽334$ are missing because the corresponding isotopes are predicted to be unbound. According to the evolution of two-neutron separation energy and neutron Fermi energy, the isotopes at $336⩽N⩽350$ are bound again, which shows a similar behavior with the “re-entrant binding” in Refs. [2,57,58,79], and forms a “peninsula” beyond the primary two-neutron drip line at $N=322$.

From the evolution of deformation, the information on nuclear structure can be extracted. Before focusing on the possible proton magicity, we first discuss the neutron structure, which is more obvious in the evolution with the increase of neutron number and can serve as a comparison for our discussion on the proton structure in the next paragraph. It is noted that except for a few spherical ones near $N=258$ and 350, all remaining Ubh isotopes are deformed. The spherical Ubh isotopes near $N=258$ and 350 correspond to new possible superheavy magic numbers. Sudden changes of ${\beta }_2$ occur near $N=208,240,258,268,292$ and 338, which are also related to the evolution of shell structure and may lead to shape coexistence [60,61]. Taking $N=266$ to 268 as an example, in Figure 2 a sudden change of deformation occurs between them. Figure 3 shows their potential energy curves from deformation-constrained DRHBc calculations by the solid curves, and the black squares represent the corresponding deformation minima from unconstrained DRHBc calculations. Both of them have two major minima at ${\beta }_2\approx -0.2$ and 0.4, and a minor minimum occurs at ${\beta }_2\approx 0.06$. These two major minima correspond to similar energies, supporting the possible shape coexistence. At $N=266$ the oblate minimum is lower, and at $N=268$ prolate minimum is lower, leading to the sudden change of ${\beta }_2$ at $N=268$. In addition, considering that most Ubh isotopes are prolate rather than oblate, the prolate-shape dominance still works in the superheavy region [63]. Therefore, the evolution of neutron structure for Ubh isotopes is reflected on the quadrupole deformation, and the sphericity near $N=258$ and 350 are significant signals corresponding to possible neutron magicity [80]. It should be also mentioned that Figure 2 focuses on the axial quadrupole deformation. When one wishes to study general stability properties in nuclear physics, especially of heavy nuclei, more quantities and effects, such as the deformation parameters with orders higher than quadrupole one, as well as the degrees of freedom from spatial-reflection asymmetry and triaxial shapes, also play important roles and should be taken into account [19,20,28,30].

Then we come to the proton structure. As has been discussed above, a na$\ddot{\mathrm{v}}$e analysis suggests $Z=126$ to be a possible traditional proton magic number. The deformation shown in Figure 2 is significantly different with the behavior in the proton-magic isotopes with $Z⩽120$. As shown in Ref. [66], in the isotopic chains with proton magic numbers $Z=8,20,28,50$, almost all nuclei are spherical, and in the isotopic chain with $Z=82$, about half of the isotopes are spherical. In comparison, for $Z=126$ here, only five isotopes near predicted neutron magic numbers $N=258$ and 350 are spherical. In conclusion, not only the evolutions of $S_{2p}$ and ${\delta }_{2p}$ depicted in Figure 1, but also the behavior of ${\beta }_2$ in Figure 2, do not support the proton magic number $Z=126$. It should also be noted that the conclusion drawn here is limited by the choice of one single axial deformation parameter ${\beta }_2$. Higher orders such as ${\beta }_4$ and ${\beta }_6$, as well as the reflection asymmetry degrees of freedom such like ${\beta }_3$ and ${\beta }_5$ which has not been incorporated in the DRHBc theory yet, also impact the calculations [81,82], and are expected to be studied in future works.

For magic nuclei, a notable result from RHB calculations is that a large gap exists in the single-particle spectrum, which further implies that the occupied nucleons cannot be scattered into higher orbits, resulting in the vanishing of pairing energy. Conversely, the sudden disappearance of pairing energy may serve as a theoretical signal to aid in the search for possible magic numbers. Figure 4 shows the neutron and proton pairing energies as functions of the neutron number in the DRHBc calculations for even-even ${}_{126}\mathrm{Ubh}$ isotopes. The evolution of neutron pairing energy is smooth at $268⩽N⩽320$, except a few turning points, while at other neutron numbers, a significant staggering pattern occurs between zero and nonzero values, such as the results near $N=240$ and 340. Since we only consider even-even Ubh nuclei in this work, the staggerings in neutron pairing energy are irrelevant to odd-even effect. The zero values at $N=258$ and 350 correspond to the predicted neutron magic numbers, while those at $N=232,240,248,320,336,344$ reflect the evolution of neutron shell structure in deformed nuclei.

The evolution of proton pairing energy is generally smooth for most Ubh isotopes, except for a few turning points. Its behavior shows no notable staggering, and never vanishes, which is remarkably different from the neutron one. For the spherical nuclei near $N=258$ and 350, the corresponding proton pairing energies are significantly larger than their neighbors. Therefore, the analysis on pairing energies does not support the proton magic number $Z=126$.

In order to further discuss the possibility of proton magicity at $Z=126$, taking the spherical isotope ${}_{126}{}^{384}\mathrm{Ubh}$ as an example, Figure 5 shows the single-proton levels around the Fermi energy ${\lambda }_p$ in the canonical basis. If we neglect the pairing correlation and let the proton occupies the orbits one by one, $Z=126$ would fully occupy the $i_{11/2}$ orbit at ${ϵ}_p=-7.847$. This orbit is very close to the $p_{1/2}$ above and the $p_{3/2}$ below, resulting in no discernible gap formed in the single-proton spectrum. Although the neutron magic number $N=126$ has been confirmed and successfully reproduced by the DRHBc calculations [58], $Z=126$ is not supported to be a magic number here.

It should also be noticed that two significant gaps appear at $Z=120$ and 138, consistent with the prediction of proton magic number $Z=120$ from the evolution of $S_{2p}$ and ${\delta }_{2p}$ in Figure 1. As has been discussed above, the DRHBc calculations have self-consistently reproduced all traditional magic numbers, and the prediction on $Z=120$ and 138 is also consistent with existing literature, such as Refs. [44,83]. We also expect further DRHBc studies on other superheavy isotopic chains to confirm the prediction of these new magic numbers.

5. Summary

The even-even ${}_{126}\mathrm{Ubh}$ isotopes are calculated based on the DRHBc calculations with PC-PK1. The evolutions of quadrupole deformation and pairing energies for neutron and proton are analyzed to study the possible nuclear magicity. Spherical shape occurs and neutron pairing energy vanishes at $N=258$ and 350, supporting that they are neutron magic numbers, while the proton pairing energy never vanishes in Ubh isotopes. In the single-proton spectrum, there is no discernible gap at $Z=126$, while significant gaps appear at $Z=120$ and 138. Therefore, $Z=126$ is not supported as a proton magic number, while $Z=120$ and 138 are suggested as candidates of proton magic numbers.

We expect further DRHBc studies on other superheavy isotopic chains to provide the corresponding ground-state properties, so that we can extract the two-proton gaps, which can be used as one of the signatures for magic numbers [44,83], to confirm the conclusions on proton magic numbers in this work. In addition, in the present studies we have only considered the ground-state static properties related to energies and shapes, while the $\alpha$, $\beta$ decays or fissions are not discussed yet. It should also be noted that the conclusion drawn here is limited by the choice of one single axial deformation parameter ${\beta }_2$, while higher orders such as ${\beta }_4$ and ${\beta }_6$ are not discussed in this work, and the reflection asymmetry degrees of freedom such like ${\beta }_3$ and ${\beta }_5$ is not incorporated in the present DRHBc theory. They may significantly influence the stability in superheavy nuclei, and their effects are expected to be investigated in future works.