Ground-state properties of finite nuclei in relativistic Hartree-Bogoliubov theory with an improved quark mass density-dependent model

Nuclear energy density functionals (EDFs), grounded in density functional theory and the mean-field approach, offer a universally applicable and powerful framework for describing the properties of finite nuclei and infinite nuclear matter. Over the past several decades, a variety of EDFs have been proposed, which can be classified into two primary categories: non-relativistic and relativistic approaches. Within the non-relativistic domain, the most successful EDFs are the Hartree-Fock models based on density-dependent forces, such as the Skyrme force for zero-range interactions [1, 2] and the Gogny force for finite-range interactions [3]. As a relativistic version of the EDFs, the relativistic mean field (RMF) model describes nucleon interactions through the exchange of various virtual mesons [4−15]. Naturally incorporating relativistic effects, the RMF model treats spin-orbit coupling in a fully self-consistent manner, without requiring additional adjustable parameters. By fitting a few free parameters to experimental data, both non-relativistic and relativistic mean-field models enable an insightful understanding of a wide range of nuclear properties[16−47], including phenomena such as the nuclear skin and halo [16−27], shell effects [28−31], shape coexistence [32−37], exotic nuclear structures [38, 39], and the proton and neutron drip lines [40−43]. These models also shed light on various properties of neutron stars [48−52].

The RMF framework has evolved into a wide range of theoretical formulations to date. In the original Walecka model [4], the nuclear interaction is described via RMF theory with nucleons interacting through the exchange of scalar (σ) and vector (ω) mesons. Subsequently, Boguta and Bodmer extended the Walecka model by incorporating cubic and quartic self-interactions of the σ-meson, thereby enhancing the model's description of the nuclear matter incompressibility and surface characteristics [53]. The self-coupling and cross-coupling terms of the vector fields have also been systematically investigated [8, 54]. Notably, the introduction of isoscalar-isovector (ω-ρ) coupling significantly influences the density dependence of the symmetry energy [13, 55−57]. Although the standard approach outlined above is widely used, several studies have extended it to explore modifications and alternatives. For instance, within the standard RMF framework, the derivative couplings of the scalar and vector fields have been proposed [58, 59]. Relative to the nonlinear meson self-interactions, relativistic models that explicitly incorporate density-dependent meson-nucleon couplings have seen significant advancements [60−64]. Besides these, the point-coupling model, which describes nucleon interactions via relativistic zero-range effective interactions, also marks a key progression within the RMF framework [65−68].

Alternatively, various effective models incorporating sub-nucleonic degrees of freedom have also been successfully developed. A prominent example is the quark-meson coupling (QMC) model, originally formulated by Guichon in 1988 [69]. This innovative approach couples the meson fields directly to the confined quarks within nucleons, establishing a theoretical framework that enables systematic investigation of nucleon modifications in the nuclear medium and provides critical insights into the fundamental properties of nuclear matter [70−74]. Other theoretical frameworks based on quark-meson degrees of freedom include the chiral SU(3) quark model [75−78], the quark mean field (QMF) model [79−81], and the improved quark mass density-dependent (IQMDD) model [82−84], among others. The IQMDD model is developed on the basis of the quark mass density-dependent (QMDD) model proposed by Fowler, Raha and Weiner [85]. According to the QMDD model, the masses of the up (u), down (d), and strange (s) quarks, as well as their corresponding antiquarks, are given by

where $ n_B $ denotes the baryon number density, $ m_{s0} $ represents the current mass of the strange quark, and B is the bag constant. Since perturbative QCD fails to provide a confinement solution for quarks, a confinement potential, typically proportional to r (or $ r^2 $), is introduced into the quark system to prevent quarks from reaching regions that are infinitely distant or excessively large [79, 80]. In the QMDD model, the confinement mechanism is manifested by requiring the mass of isolated quarks to approach infinity, rendering the vacuum incapable of supporting them [86, 87]. As the volume approaches infinity or the density tends to zero, the quark mass diverges, which is precisely the behavior characterized by equations (1) and (2). This confinement mechanism is analogous to that in the MIT bag model [88−90]. Building upon the QMDD framework, the IQMDD model extends the original formalism by systematically incorporating σ, ω, and ρ meson fields. This generalization establishes a self-consistent coupling between quark degrees of freedom and mesonic fields throughout the nuclear medium, thereby enabling rigorous nuclear many-body calculations within the mean-field approximation [91].

On the other hand, the construction and operation of new-generation radioactive ion beam (RIB) facilities have produced an increasing number of nuclei far from β-stability valley, opening up new frontiers in nuclear physics that expand our understanding of nuclear phenomena from stable nuclei to exotic ones [92−94]. Exotic nuclei near the drip lines are weakly-bound systems, where the Fermi surface of protons (or neutrons) is close to the particle continuum. This leads to an increased scattering of Cooper pairs into the continuum states as a result of pairing correlations, and furthermore, the conventional BCS model with the monopole pairing force can only provide a poor approximation in this case [95]. For weakly bound nuclei, a unified and self-consistent treatment of the mean field and pairing correlations becomes essential. The development of the RHB model has provided an effective approach to address this issue [33, 96−99]. In the RHB model, a unified description of particle-hole (ph) and particle-particle (pp) correlations is given, and the continuum effects are also effectively taken into account in the coordinate representation [100−105]. In this paper, RHB calculations incorporating quark-meson coupling are presented within the IQMDD framework, serving as a natural extension of the IQMDD model. With this new model, we systematically calculated the ground-state properties of even-even nuclei with proton numbers ranging from Z=8 to Z=118 where experimental data are available. This paper is organized as follows. The main formulas of the IQMDD+RHB model are presented in Section 2. In Section 3, the model parameters and numerical results are presented and discussed. Finally, a summary is given in Section 4.

The effective Lagrangian density of the IQMDD model can be expressed as

The strength tensors for the vector mesons and the electromagnetic field are defined as follows: $ \Omega_{\mu\nu}=\partial_\mu\omega_\nu- \partial_\nu\omega_\mu $, $ G_{\mu\nu}=\partial_{\mu}\rho_{\nu}-\partial_{\nu}\rho_{\mu} $ and $ F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu} $. The quark mass $ m_q $ ($ q = u $, d), characterized by its density dependence, is given by equation (1). The self-interaction potential of the σ field is expressed in the following form

The bag constant B is introduced such that $ U(\sigma_{\nu})=0 $, where $ \sigma_{\nu} $ denotes the value of the sigma field at which the potential $ U(\sigma) $ attains its absolute minimum. The effective quark mass $ m_q^* $ is given by

In nuclear matter, three quarks form a soliton bag, and the effective nucleon mass is derived from the bag energy and is expressed as

where $ K^q_F $ is the Fermi momentum of the quarks, $ \Gamma_q $ denotes the quark degeneracy. $ dN_q/dk $ is the density of states for various types of quarks in a confined space. The expressions for $ dN_q/dk $ and $ K^q_F $ used in this work can be found in Ref. [106]. The bag radius R can be obtained from the equilibrium condition for the nucleon bag

Our focus now shifts to examining finite nuclei within the IQMDD model. Analogous to the QMC model [74], one can formulate a relativistic Lagrangian density at the hadronic level:

Starting from the Lagrangian density (8), the RHB equation for the nucleons, obtained through a self-consistent treatment of the mean field and pairing correlations, is presented as [107]

Here, λ represents the chemical potential, and $ E_k $ denotes the quasiparticle energy. The single-nucleon Dirac Hamiltonian, $ h_D $, is given by

The vector potential reads

where $ \boldsymbol{I} $ represents the identity matrix. Formally, the equations of motion for the mesons and photon are given by

Here, $ \rho_s $, $ \rho_\upsilon $, $ \rho_0^T $, $ \rho_3 $ and $ \rho_p $ represent the densities of scalar, vector, tensor, the third component of isovector, and proton, respectively. The relativistic pairing field Δ in Eq. (9) can be expressed as

The pairing tensor is written as

In the two-body pairing channel, the finite range force is used, and its form in coordinate space is given by [108],

where $ \boldsymbol{r} $ and $ \boldsymbol{R} $ stand for the relative and center-of-mass coordinates of the two-nucleon system, respectively. The values of parameters G and a are adopted from Ref. [108].

The coupled equations mentioned above can be solved through a self-consistent iterative scheme, where the wave functions are expanded in terms of the harmonic oscillator basis functions. In our calculation, the number of oscillator shells is chosen as $ N_f=12 $ and $ N_b=20 $. From this solution, we can calculate the total binding energy, which has the following form:

Here, $ E_{\mathrm{pair}} $ denotes the nuclear ground-state pairing energy, which is given by

The microscopic center-of-mass correction, which is employed in our calculations, is given as:

where $ P_{\mathrm{c.m.}} $ denotes the total momentum as viewed from the center-of-mass frame.

Prior to conducting the numerical calculations, we first examine the parameters of this model. The parameter set IQMDD-2*, incorporating both the quark-ω tensor coupling and the non-linear ω-ρ interaction term, has been implemented to investigate the properties of finite nuclear systems [109]. This parameter set was initially proposed a decade ago, and subsequent experimental advancements have yielded substantial new data on nuclear mass and neutron skin thickness measurements. To explore potential improvements in finite nuclear property predictions, we conduct a global fitting procedure that combines the simulated annealing algorithm with the Levenberg-Marquardt method, resulting in the parameterization referred to as IQMDD3. The fitted nuclear properties include the binding energies and charge radii of twenty even-even nuclei: ¹⁶O, ⁴⁰Ca, ⁴⁸Ca, ⁵⁸Ni, ⁸⁶Kr, ⁸⁸Sr, ⁹⁰Zr, ⁹²Mo, ¹¹⁶Sn, ¹²⁴Sn, ¹³²Sn, ¹³⁴Te, ¹³⁸Ba, ¹⁴²Nd, ¹⁴⁴Sm, ¹⁴⁸Gd, ¹⁸⁶Pb, ²⁰⁸Pb, ²¹⁰Pb and ²¹⁰Po. Spherical RHB calculations are employed in the fitting process. The experimental binding energy values are adopted from the atomic mass evaluation [110], and charge radii from Ref. [111].

The interaction strengths between the nucleon and the ω- and ρ-mesons are selected as: $ g_{\omega}=3g_{\omega}^q $ and $ g_{\rho}=g_{\rho}^q $ [71]. The bag constant is set to $ B=174.0 \; \mathrm{MeV}\cdot\mathrm{fm}^{-3} $. The optimized parameters obtained for the IQMDD3 Lagrangian parameterization are presented in Table 1, with comparative results from the previous IQMDD-2*. To achieve better agreement with experimental neutron skin thickness measurements of ²⁰⁸Pb, the ω-ρ coupling constant (η) is set to zero in the IQMDD3 parameterization. Recent parity-violating electron scattering experiments at JLab have achieved a precision measurement of the neutron skin thickness in ²⁰⁸Pb, yielding a result of $ R_n-R_p=0.283\pm0.071 $ fm [113]. Using dispersive optical model (DOM) analysis, Ref. [114] determined the neutron skin thickness of ²⁰⁸Pb to be $ 0.12-0.25 $ fm. Additional experimental results for this quantity are presented in Fig. 12 of Ref. [115]. The calculated value for ²⁰⁸Pb from the IQMDD3 model is $ R_n-R_p=0.24 $ fm, which falls within the experimental uncertainty range reported in Refs.[113, 114]. In Table 2, we present the IQMDD3 predictions for the total binding energies, charge radii, and neutron-skin thickness of the nuclei included in the fit, along with the calculated results for other isotopic chains. All calculations are performed within a spherical RHB approach. For the twenty nuclei used in the fitting, the rms deviations from IQMDD3 are 0.023 MeV for the binding energy per nucleon, 2.33 MeV for the total binding energy, and 0.01 fm for the charge radius.

	$ m_\sigma $	$ m_\omega $	$ m_\rho $	$ g^q_\sigma $	$ {g^q_\omega} $	$ {g^q_\rho} $	b	$ f_\omega $	η
IQMDD-2*	508.0	782.5	763.0	5.09	2.9796	9.7001	–3400.0	2.1	0.03
IQMDD3	498.5998	782.5	763.0	4.9667	2.9404	9.2314	–3114.1540	2.0993	0.0

Table 1.  Model parameters of the effective interactions IQMDD3 and IQMDD-2*. The parameter b and the masses of the mesons $ m_\sigma $, $ m_\omega $ and $ m_\rho $ are measured in MeV.

Nucleus	B.E. (MeV)	$ R_c $ (fm)	$ R_{np} $ (fm)
16O	128.328 (127.619)	2.722 (2.699)	−0.03
40Ca	341.070 (342.052)	3.464 (3.478)	−0.05
48Ca	416.619 (416.001)	3.493 (3.477)	0.20
50Ti	438.112 (437.786)	3.574 (3.570)	0.13
52Cr	455.211 (456.352)	3.656 (3.645)	0.06
54Fe	469.653 (471.765)	3.704 (3.693)	0.00
58Ni	504.782 (506.460)	3.782 (3.776)	0.00
68Ni	590.678 (590.408)	3.899 (3.892$ ^a $)	0.21
86Kr	751.827 (749.235)	4.189 (4.184)	0.17
88Sr	770.741 (768.468)	4.235 (4.224)	0.13
90Zr	786.180 (783.897)	4.274 (4.269)	0.10
92Mo	799.110 (796.511)	4.319 (4.315)	0.06
112Sn	953.695 (953.525)	4.583 (4.595)	0.09
116Sn	989.254 (988.682)	4.616 (4.625)	0.14
124Sn	1051.822 (1049.958)	4.671 (4.674)	0.22
132Sn	1104.904 (1102.843)	4.728 (4.732)	0.30
134Te	1126.579 (1123.408)	4.765 (4.757)	0.26
136Xe	1145.636 (1141.882)	4.805 (4.796)	0.23
138Ba	1162.421 (1158.292)	4.843 (4.838)	0.20
140Ce	1176.638 (1172.683)	4.880 (4.877)	0.17
142Nd	1188.969 (1185.136)	4.915 (4.912)	0.14
144Sm	1199.256 (1195.730)	4.953 (4.952)	0.12
148Gd	1222.693 (1220.754)	5.003 (5.008)	0.11
184Pb	1435.552 (1432.022)	5.386 (5.393)	0.06
186Pb	1454.809 (1451.793)	5.397 (5.403)	0.08
204Pb	1608.613 (1607.506)	5.491 (5.480)	0.20
208Pb	1636.544 (1636.430)	5.512 (5.501)	0.24
210Pb	1645.453 (1645.553)	5.528 (5.521)	0.25
214Pb	1661.461 (1663.293)	5.562 (5.558)	0.27
210Po	1647.514 (1645.213)	5.548 (5.570)	0.21
aAvailable experimental data taken from [112];the remaining experimental charge radii are from [111].

Table 2.  Total binding energies B.E., charge radii $ R_c $, and neutron-skin thickness $ R_{np} $($ R_{np} $=$ R_{n}-R_{p} $) calculated with IQMDD3, compared to experimental data [110−112] (values in parentheses).

Table 3 presents the spin-orbit splittings for ¹⁶O, ⁴⁰Ca, ⁴⁸Ca, ⁵⁶Ni, ⁷⁸Ni, ¹⁰⁰Sn, ¹³²Sn and ²⁰⁸Pb calculated with the DD-ME2, IQMDD-2*, and IQMDD3 parameter sets. The inclusion of tensor coupling significantly impacts the IQMDD model's description of nuclear spin-orbit splittings. As the tensor coupling strengths ($ g^q_\omega $ and $ f_\omega $) of the IQMDD3 parameter set are slightly lower than those of IQMDD-2*, the calculated spin-orbit splittings for IQMDD3 are overall smaller than those for IQMDD-2*, as indicated in Table 3.

Additionally, we present the findings of the bulk characteristics for infinite nuclear matter in Table 4. It presents corrected predictions for nuclear matter properties based on the IQMDD-2* parameter set as reported in Ref. [109]. Compared with IQMDD-2*, the IQMDD3 parameterization, which lacks isoscalar-isovector coupling, yields a larger symmetry energy at saturation density. For the six parametrizations in Table 4, the neutron-skin thickness of ²⁰⁸Pb increases monotonically from 0.160 fm to 0.285 fm as the symmetry energy slope increases from 47.2 MeV to 122.6 MeV. This positive correlation, previously highlighted by the literature [117], indicates that a larger symmetry energy slope L implies higher neutron matter pressure and a thicker neutron-skin thickness in ²⁰⁸Pb.

Model	$ \rho_0 $	$ \varepsilon_0 $	S	L	K	$ R_{np} $
IU-FSU	0.155	−16.40	31.30	47.2	231.3	0.160
FSU	0.148	−16.30	32.59	60.5	230.0	0.205
IQMDD-2*	0.155	−16.38	34.16	72.8	283.3	0.212
IQMDD3	0.153	−16.41	36.00	100.1	292.8	0.236
NL3	0.148	−16.24	37.40	118.5	271.8	0.277
NL3*	0.150	−16.31	38.68	122.6	258.1	0.285

Table 4.  Bulk properties for infinite nuclear matter at saturation density $ \rho_0 $ (in fm$ ^{-3} $) across different models [10, 13−15]. The quantities $ \varepsilon_0 $, S, L and K (all in MeV) represent the binding energy per nucleon, the symmetry energy, the slope of the symmetry energy and incompressibility coefficient of symmetric nuclear matter at saturation density. The neutron-skin thickness of ²⁰⁸Pb, represented by $ R_{np} $ (in fm), is also included.

We conduct systematic calculations for even-even isotopes from Oxygen (Z=8) to Oganesson (Z=118), focusing on nuclei with available experimental data. The differences in the total binding energy between calculated and experimental values for even-even nuclei spanning Z=8–118 are systematically presented in Fig. 1. The relevant experimental data are obtained from AME2020 atomic mass evaluation [110]. For the specified mass region, our calculations encompass all 868 even-even nuclei in AME2020, containing 646 experimentally measured and 222 estimated mass values. Note that the binding energies herein and all subsequent results are obtained from axially deformed RHB calculations. As shown in Fig. 1, IQMDD3 and IQMDD-2* exhibit similar overall performance for total binding energy. In the fitting process, the binding energy per nucleon was employed as a constraint instead of the total binding energy. This adjustment effectively increased the fitting weight for light nuclei but led to a moderately larger deviation in the total binding energy of heavy nuclei compared to experimental values. This is reflected in the relatively large deviation in the total binding energy for heavy nuclei shown in Fig. 1. Statistical analysis reveals that the rms deviation of binding energy per nucleon is $ \sigma_{B/A}=0.034 $ MeV for IQMDD3 and $ \sigma_{B/A}=0.036 $ MeV for IQMDD-2*. For total binding energy, the corresponding values are $ \sigma_{B}=2.893 $ MeV for IQMDD3 and $ \sigma_{B}=2.984 $ MeV for IQMDD-2*. Additionally, statistical analysis shows that both models perform better in the medium and heavy nucleus regions compared to the light nucleus region in terms of the average binding energy per nucleon. When the analysis is confined to the range from Calcium (Z=20) to Oganesson (Z=118), the rms deviation decreases significantly to $ \sigma_{B/A} =0.020 $ MeV for IQMDD3 and $ \sigma_{B/A} =0.022 $ MeV for IQMDD-2*.

Figure 1.  (color online) Differences in total binding energy between the calculations and experimental data for even-even nuclei using the IQMDD3 and IQMDD-2*. To guide the eye, the dashed frames correspond to the positions of the magic numbers 8, 20, 28, 50, 82, and 126.

For comparison, the typical binding energy rms deviations of established EDFs are 1.518 MeV for PC-PK1 by including the rotational correction energies [118], 1.74 MeV for QMCπ-Ⅲ-T [119], 2.96 MeV for NL3* [120], 2.39 MeV for DD-ME2 [120], 2.29 MeV for DD-MEδ [120], 2.01 MeV for DD-PC1 [120], and 1.339 MeV for 91 spherical nuclei for PC-L3R [121]. In terms of predicting nuclear binding energies, the performance of IQMDD3 is overall comparable to that of NL3*.

Fig. 2 illustrates the theoretical predictions of ground-state quadrupole deformation parameters $ \beta_2 $ for all the isotopic chains examined in this work. The deformation parameter $ \beta_2 $ is derived from the calculated nucleon quadrupole moments Q via [7]:

Figure 2.  (color online) Quadrupole deformation parameters for even-even nuclei with proton numbers $ 8\leq Z\leq118 $, calculated using the RHB approach with the IQMDD3 interaction.

with $ R_0=1.2A^{1/3} $(fm). The $ \beta_2 $ parameter, a dimensionless quantity with an absolute value less than 1, is characterized by relatively large empirical uncertainties. Given that our primary objective is to investigate the systematic trends and global characteristics of nuclear deformation patterns predicted by the IQMDD3 model, a comparison with experimental values is not included in the present analysis. As illustrated in Fig. 2, nuclei with spherical configurations are predominantly located in the region around the magic numbers. Nuclei with proton and neutron numbers situated between two magic numbers tend to exhibit larger deformations. Additionally, for nuclei located in the middle of the major shells, deformations in the light and medium-mass regions are generally more pronounced than those in heavy and superheavy nuclei.

For the neutron magic numbers N=8, 20, 50, 82, and 126, the corresponding isotonic chains are predominantly composed of spherical nuclei. However, the isotonic chain with neutron number N=28 exhibits significant shape shifts. Isotopic chains characterized by proton magic numbers (Z=8, 20, 50, 82) generally maintain spherical configurations, whereas the Z=28 isotopic chain exhibits moderate quadrupole deformation, which is not as pronounced as that in the isotonic chain with N=28.

For nuclei with neutron numbers Z>50, those located in the central region between major nuclear shells predominantly exhibit prolate deformation. A characteristic shape transition occurs as the neutron number approaches magic numbers, where nuclear deformation evolves from prolate to oblate configurations. Within the proton number range of $ 28<Z<50 $, oblate nuclear configurations are predominantly found near the central region of the major shells. Similarly, the range of $ 8<Z<20 $ also contains oblate nuclei. However, according to the IQMDD3 model, no significantly deformed oblate nuclei are predicted to exist within the $ 20<Z<28 $ range.

Besides nuclear deformations, the charge radius is also a crucial quantity for characterizing finite nuclei, and is expressed as:

where the factor 0.64 corrects for the finite size of the proton charge distribution. The deviations between the theoretical values derived from the IQMDD3 model and the experimental data cited in Ref. [111, 112] are depicted in Fig. 3. It shows that the IQMDD3 predictions are in reasonable agreement with the experimental charge radii, with calculated deviations typically within the $ \pm0.04 $ fm range. For the 368 even-even nuclei with known charge radii, the rms deviation between the theoretical and experimental values is $ \sigma = 0.022 $ fm, compared to $ \sigma = 0.021 $ fm for IQMDD-2*.

Figure 3.  (color online) Discrepancies between IQMDD3+RHB calculations and experimental data for even-even nuclei with measured charge radii.

The most significant discrepancies occur in the Mg(Z=12), Si(Z=14), and Cm(Z=96) isotopic chains, where theoretical values exceed experimental values by over 0.05 fm. For ²⁶Mg, ²⁸Si, and ³⁰Si, the theoretical overestimation of the charge radii relative to experimental values can be attributed to the significant oblate deformation ($ |\beta_2|>0.2 $) predicted in their ground states. Experiments reveal that the Cm (Z=96) isotopic chain exhibits anomalous charge radii, which are smaller than those of Pu (Z=94) [111]. However, IQMDD3 calculations indicate that charge radii typically increase with proton number in an isotonic chain. This results in a noticeable discrepancy between theoretical predictions and experimental observations for the charge radii of the Cm isotopic chain. Notably, this anomalous charge radius behavior exhibited by U-Pu-Cm isotopes has already been extensively discussed in Refs. [118, 120, 122].

Two-nucleon separation energies ($ S_{2n} $ and $ S_{2p} $) offer essential insights into the investigation of neutron and proton drip lines, nuclear shell evolution, and nucleon radioactivity phenomena [123]. They can be expressed as:

where $ B(A,Z) $ represents the binding energy of a nucleus with mass number A and proton number Z, calculated from equation (16). In Fig. 4, the predicted two-neutron separation energies, $ S_{2n} $, for even-even nuclei using the IQMDD3 parameters are depicted, along with the discrepancies between the theoretical predictions and experimental data from Ref. [110]. The left panel of Fig. 4 shows that, for a given isotopic chain, $ S_{2n} $ decreases with increasing neutron number, while for a given isotonic chain, $ S_{2n} $ increases with the increase in proton number. Meanwhile, as the proton number increases, the two-neutron separation energy $ S_{2n} $ on the proton-rich side will gradually decrease. The theoretical values of $ S_{2n} $ calculated by IQMDD3 are compared with the experimental values. As shown in the right panel of Fig. 4, most of the differences between theory and experiment remain within $ \pm2 $ MeV, and for the studied nuclei, the rms of the differences is determined to be 0.984 MeV.

Figure 4.  (color online) Left: The two-neutron separation energies of even-even nuclei calculated by IQMDD3. Right: The deviation between the theoretical and experimental values of the two-neutron separation energies.

It is quite interesting to note that, as seen from the left panel of Fig. 4, there are pronounced sudden drops in $ S_{2n} $ at the positions of N=50, 82 and 126. Similar sudden drops can also be observed at N=20 and 28. These abrupt decreases in $ S_{2n} $ indicate the emergence of neutron shells, which are well-reproduced by IQMDD3.

Fig. 5 presents the systematic trends of two-proton separation energy ($ S_{2p} $) for even-even nuclei, accompanied by the differences between theoretical predictions and experimental measurements. In contrast to the two-neutron separation energy ($ S_{2n} $), for a given isotopic chain, $ S_{2p} $ exhibits an increasing trend with neutron number, while in isotonic chains, it demonstrates a decreasing trend with proton number. Furthermore, the data indicate a gradual reduction in $ S_{2p} $ values within neutron-rich regions as the proton number increases. The right panel of Fig. 5 demonstrates that the majority of theoretical-experimental deviations fall within a range of $ \pm2 $ MeV. For the nuclei under investigation, the rms deviation of $ S_{2p} $ between theory and experiment is calculated to be 0.848 MeV, indicating acceptable agreement between theoretical predictions and experimental observations.

Figure 5.  (color online) Left: Theoretical predictions of two-proton separation energies for even-even nuclei obtained from IQMDD3 calculations. Right: The discrepancy between theoretical calculations and experimental measurements of the two-proton separation energies.

As observed in the left panel of Fig. 5, there are notable sudden drops in $ S_{2p} $ at the proton numbers Z=20, 28 and 50. These abrupt decreases in $ S_{2p} $ signify the appearance of proton shells. For the isotope chain with Z=82, the variation in $ S_{2p} $ is relatively subtle. Nevertheless, the two-proton shell gaps, which will be discussed in the following subsection, provide a clearer manifestation of this change.

The two-nucleon energy gaps, essential for examining the shell effects in finite nuclei, are defined as follows:

Fig. 6 displays the two-neutron gaps ($ \delta_{2n} $) for even-even nuclei in the region of $ 8\leq Z\leq114 $. It is shown that $ \delta_{2n} $ exhibits significant sharp changes at neutron numbers N=8, 20, 28, 50, 82, and 126, which can be regarded as one of the signatures of magic numbers. For the magic number N=28, $ \delta_{2n} $ becomes relatively small in both proton-rich and neutron-rich regions, suggesting that the N=28 shell closure may be quenched in these areas.

Figure 6.  (color online) Two-neutron shell gaps $ \delta_{2n} $ for even-even nuclei, as calculated with the IQMDD3 parameterization.

Fig. 7 presents the two-proton gaps ($ \delta_{2p} $) for even-even nuclei in the region of $ 10\leq Z\leq116 $. It can be observed that $ \delta_{2p} $ exhibits distinct peaks at proton numbers Z=20, 28, and 50, which is consistent with the conclusions regarding shell closures shown in Fig. 5. For Z=82, $ \delta_{2p} $ also shows a significant sharp change, whereas the change in $ S_{2p} $ is less evident in Fig. 5. This indicates that, compared to the two-proton separation energy, the two-proton gap provides a more intuitive indication of potential magic numbers. Additionally, it is noteworthy that a relatively large $ \delta_{2p} $ is also observed at Z = 14, suggesting that Z = 14 may correspond to a subshell closure.

Figure 7.  (color online) Two-proton shell gaps $ \delta_{2p} $ in even-even nuclei from IQMDD3 calculations.

α-decay usually occurs in heavy and superheavy radioactive nuclides. The α-decay energy, denoted as $ Q_{\alpha} $, is expressed as:

Fig. 8 presents the discrepancies between the $ Q_{\alpha} $ values calculated using the IQMDD3 model and the experimental data from Ref. [110]. Among the 405 studied even-even nuclei with positive $ Q_{\alpha} $ values, relatively larger deviations are observed in specific nuclei of the Z=52 and Z=94 isotopic chains, as well as in some nuclei along the N=84 isotonic chain. Overall, most of the discrepancies are within $ \pm1.5 $ MeV, and the rms of the discrepancies is 0.733 MeV. Recently, machine learning approaches have been successfully employed in the description of α-decay energies, achieving a remarkably high level of precision [124, 125].

Figure 8.  (color online) The discrepancies between the IQMDD3-calculated and experimental values of α decay energy $ Q_{\alpha} $.

In this work, we have developed a relativistic Hartree-Bogoliubov (RHB) model incorporating quark-meson couplings and established a new parameter set named IQMDD3 by fitting to the binding energies and charge radii of twenty selected even-even nuclei. Using this parameter set, we systematically studied the ground-state properties of 868 even-even nuclei, with proton numbers from Z=8 to Z=118. For the twenty nuclei employed for fitting, the rms deviation of the binding energy per nucleon obtained with the IQMDD3 model is 0.023 MeV, that of the total binding energy is 2.33 MeV, and that of the charge radius is 0.01 fm. The predicted neutron skin thickness for ²⁰⁸Pb with this parameter set is 0.24 fm, with the symmetry energy and its slope at saturation density being 36.00 MeV and 100.1 MeV, respectively.

For the 868 even-even nuclei studied, the rms deviation of the binding energy per nucleon is 0.034 MeV, and the rms deviation of the total binding energy is 2.893 MeV. We investigated the systematic trends and general features of quadrupole deformations in these even-even nuclei. Results show that nuclei in the magic number regions and those in their vicinity are mostly spherical, except for the N=28 isotonic chain, which exhibits significant shape shifts. In contrast, nuclei far from the major shells exhibit larger deformations, with the largest deformations often occurring in the central regions between major shells. The analysis of charge radii shows that, for the 368 studied even-even nuclei with available experimental data, the rms deviation of the charge radii is 0.022 fm.

We systematically examined the two-nucleon separation energies, as well as the two-nucleon shell gaps. For the same isotopic chain, the two-neutron separation energy is larger in proton-rich regions than in neutron-rich regions, while the two-proton separation energy is the opposite, being larger in neutron-rich regions than in proton-rich regions. The conventional magic numbers were successfully reconstructed based on the systematic trends observed in two-nucleon separation energies and two-nucleon shell gaps. Additionally, the two-proton shell gaps indicate that Z = 14 might be a subshell closure.

Finally, for nuclei with positive α-decay energies, the extracted $ Q_{\alpha} $ values were compared with the 405 available experimental data, yielding a rms deviation of $ \sigma = 0.733 $ MeV.

To further examine the validity and reliability of the IQMDD3 parameter set in describing finite nuclei, future work will employ this model to study the properties of drip-line nuclei far from the β-stability line.