Charge radii, moments and masses of mercury isotopes across the N = 126 shell closure

Combining laser spectroscopy in a Versatile Arc Discharge and Laser Ion Source, with Penning-trap mass spectrometry at the CERN-ISOLDE facility, this work reports on mean-square charge radii of neutron-rich mercury isotopes across the $N = 126$ shell closure, the electromagnetic moments of $^{207}$Hg and more precise mass values of $^{206-208}$Hg. The odd-even staggering (OES) of the mean square charge radii and the kink at $N = 126$ are analyzed within the framework of covariant density functional theory (CDFT), with comparisons between different functionals to investigate the dependence of the results on the underlying single-particle structure. The observed features are defined predominantly in the particle-hole channel in CDFT, since both are present in the calculations without pairing. However, the magnitude of the kink is still affected by the occupation of the $1i_{11/2}$ and $2g_{9/2}$ orbitals with a dependence on the relative energies as well as pairing.


I. INTRODUCTION
The kink in the relative mean square charge radii (δ r 2 ) at the N = 126 shell closure has long been considered as a benchmark for testing theoretical calculations. Traditionally the lead isotopic chain was employed [1][2][3][4][5], but new experimental results in this region revealing the systemics of other isotopic chains [6][7][8][9], mass measurements, and odd-even staggering (OES) in charge radii offer the opportunity to broaden this benchmark. The droplet model is unable to reproduce this kink because of the absence of single-particle degrees of freedom [10]. Early non-relativistic mean field approaches were also incapable of reproducing a kink at 208 Pb [11], while conversely, relativistic mean field approaches were demonstrated to be successful in doing so [2].
Two alternative modifications were suggested to correct this arXiv:2111.10468v1 [nucl-ex] 19 Nov 2021 deficiency in non-relativistic models. The first relies on the modification of the spin-orbit interaction, either through a fitting procedure (see Refs. [3,12]) or via the introduction of a density dependence (see Refs. [13,14]). This leads to a reasonable reproduction of the experimental isotope shifts (see Refs. [3,15]). The second approach (employing so-called Fayans functionals) introduces gradient terms into the pairing and surface terms of the functional [4,16,17]. This significantly improves the general description of experimental data, however, discrepancies are still apparent in the lead and tin isotopic chains [18]. Moreover, pairing becomes a dominant contributor to the kink and OES [18], in contradiction with the results obtained in relativistic Hartree-Bogoliubov (RHB) calculations with the DD-ME2 functional and nonrelativistic Hartree-Fock-Bogoliubov (NR-HFB) calculations with the M3Py-P6a functional presented in Ref. [9]. This work is an in-depth follow-up to Ref. [9], which reported on the δ r 2 of mercury isotopes across N=126 and employed these results, together with existing lead data, to compare RHB and NR-HFB approaches. A new OES mechanism was additionally suggested, related to the staggering in the occupation of the different neutron orbitals in odd-and even-A nuclei and facilitated by particle-vibration coupling (PVC) in odd-A nuclei. Here we report on the magneticdipole and electric-quadrupole moments of 207 Hg, new and improved mass measurements of 206−208 Hg and a detailed theoretical study within the RHB framework to better understand the kinks and OES in lead and mercury isotopes. Multiple state-of-the-art covariant energy density functionals (CEDFs) are employed (NL3* [19], DD-PC1 [20], DD-ME2 [21] and DD-MEδ [22]) to assess the dependence of the theoretical results on the underlying single-particle structure. The global performance of these functionals in describing ground state properties such as masses and charge radii of even-even nuclei has been tested in Refs. [23] and [24].
This article is arranged as follows. The experimental techniques are presented in Sec. II. The radiogenic production of 207,208 Hg in a molten target is discussed in Sec. III. Experimental results on mean square charge radii, -dipole and electric-quadrupole moments and masses are summarized in Sec. IV. The discussion of experimental results is presented in Sec. V. Theoretical formalism and theoretical analysis of the kinks and OES in charge radii are presented in Sec. VI, together with their dependence on the underlying single-particle structure and pairing and a comparison of experimental and calculated binding energies. Finally, we give a brief summary in Sec. VII.

II. EXPERIMENTAL TECHNIQUE
The neutron-rich mercury isotopes were studied at the ISOLDE facility [25] as part of a wider experimental campaign which investigated both ends of the isotopic chain [9,26,27]. The mercury nuclei were produced using the Isotope Separator On-Line (ISOL) method [28,29] and studied via insource resonance ionization spectroscopy [30,31] as depicted in Fig. 1(a). of mercury isotopes. b) Schematic of a VADLIS coupled to a molten lead target via a temperature controlled chimney. c) The variation in the RILIS-mode and VADIS-mode ion currents on mass A=197 for differing anode bias voltages. See text for additional details and the definition of acronyms.
A molten lead target (thickness 170 g/cm 2 ) was bombarded with 1.4-GeV protons, resulting in a cocktail of reaction products which effused via a temperature controlled chimney [32] into the anode volume of a Versatile Arc Discharge and Laser Ion Source (VADLIS) [33]. The target and ion source were biased at 30 kV and laser light from the ISOLDE Resonance Ionization Laser Ion Source (RILIS) [34] was directed into the anode volume for multi-step resonance ionization [35] of mercury isotopes. A {λ 1 |λ 2 | λ 3 } = {254 nm|313 nm|532 nm} ionization scheme [36] was applied, with the first (254 nm) resonant transition used to investigate the hyperfine structure (hfs) and isotope shifts in the 5d 10 6s 2 1 S 0 → 5d 10 6s6p 3 P • 1 atomic transition.
The ions were accelerated by the electric field resulting from the grounded extraction electrode depicted in Fig. 1(a) and Fig. 1(b) to form a 30-keV radioactive ion beam (RIB). The ISOLDE General Purpose Separator [25] was employed for mass separation before the RIB was directed to either a Faraday cup for direct ion current measurement or to the ISOLTRAP radio frequency quadrupole cooler-buncher (RFQ C-B) [37]. Downstream of the RFQ C-B, the RIB was injected into the Multi-Reflection Time-of-Flight Mass Spectrometer (MR-ToF MS) [38] for either isobaric separation and subsequent detection [39] or for mass measurements, either by measuring the time-of-flight of the ions [40] or by utilizing the downstream Penning traps [41].
The lead target-VADLIS combination was required to avoid the overwhelming isobaric francium contamination present on masses A=207, 208 when employing a standard UC x target with a hot cavity surface ion source for the laser light-atom interaction region [42]. Alternative approaches have struggled to suppress certain isotopes of francium, and would additionally be expected to result in a factor of ≈20 reduction of the signal of interest [43,44].
A schematic of the VADLIS is presented in Fig. 1(b) together with the relative bias of the components. In the standard Versatile Arc Discharge Ion Source (VADIS)-mode of operation, atoms and molecules are ionized by electrons that are emitted from the ≈2000 • C cathode and accelerated into the anode volume by a relative anode voltage of 100-200 V [45]. The selective RILIS-mode of operation was employed for this experiment, where the anode voltage is optimized for laser-ion extraction while maintaining it below what is required for significant electron impact ionization [33]. Figure 1(c) presents the on-line optimization of the RILIS-mode with radiogenically produced 197 Hg. The RILIS-mode and VADIS-mode related signals were separated by blocking and unblocking the laser light exciting the 254-nm transition. A clear maximum is visible with a near-negligible background with the anode voltage set to ≈8 V. The alternative, applying the RILIS lasers with a 100-200 V anode bias (termed RILIS+VADIS-mode) would have resulted in significant isobaric contamination and a reduced signal-to-noise ratio as a result of the competing ionization processes.
The benefits of combining the RILIS-mode of operation with the ISOLTRAP MR-ToF MS are highlighted in Fig. 2(a). Operating in RILIS-mode reduced the isobaric 208 Pb background by seven orders of magnitude compared  for a given frequency tripled laser frequency. The RILIS+VADIS-mode is dominated by VADIS ionized 208 Pb, the count rate was estimated based on Faraday cup measurements. Further discussion in the text. b) and c) time-of-flight spectra on mass A=208, measured downstream of the MR-ToF MS with the VADLIS in RILIS-mode and the lasers on resonance and off-resonance, respectively. The y-axis scales of b) and c) are identical.
with RILIS+VADIS-mode. This enabled the MR-ToF MS to be employed for selective detection and for determining the mass of 208 Hg. Time-of-flight spectra recorded on and off resonance with the MR-ToF MS are shown in Fig. 2(b) and 2(c), respectively. By applying time gates in the ToF spectra, it was possible to separate the 208 Hg signal from the remaining isobaric contamination. Considering the natural lead ( 206−208 Pb) target material used for this experiment, the creation of mercury isotopes with N ≤ 126 is comparatively well understood as the result of spallation reactions induced by the incident 1.4 GeV proton beam. However, when going beyond N = 126 ( 206 Hg) the production mechanism changes, and a range of other processes may become relevant [46,47] including secondary reactions induced by the light and energetic products of the primary spallation reactions. The production of 207 Hg via 208 Pb(n, 2p) 207 Hg is a good example of such a process, and was first reported at an ISOL facility in Ref. [48]. The mechanism for producing 208 Hg is significantly more exotic, as evidenced by a factor of 2400 decrease in the measured ion rate between 207 Hg and 208 Hg. There are a number of potential production channels including 208 Pb(t,3p) 208 Hg, 208 Pb(α,4p) 208 Hg, or reactions with radiogenically produced 209 Pb (t 1/2 ≈ 3 h) or 210 Pb (t 1/2 ≈ 22 y) which build up within the target during the experiment.
The in-target production of mercury isotopes was calculated via ABRABLA [49,50], FLUKA [51,52] and Geant4 [53][54][55] simulations. The results are assessed by considering the isotope specific extraction and ionization efficiencies (ε), determined by dividing the measured yield by the cal- culated in-target production. Figure 3 presents the relationship between ε and half-life for the mercury isotopes measured with the MR-ToF during this experimental campaign. For 202 Hg (stable) and 203 Hg (t 1/2 ≈47 days) the half-lives are set at 1×10 5 s to facilitate their inclusion. The data is fitted using Eq. (4) from Ref. [42], with the hollow data points omitted from the fits to enable them to converge. As expected, ε generally increases with increasing half-life, and stabilizes at a point where the half-life is sufficiently long compared to the release time.
All of the results broadly agree with an extraction and ionization efficiency of approximately 1% for sufficiently long-lived isotopes. ABRABLA [49] is not capable of reproducing the secondary reactions required for the production of 207,208 Hg, however, as it is commonly used for calculating in-target production, it is useful for benchmarking the other codes for this application. FLUKA [52] was found to reproduce some 207 Hg production, though based on Fig. 3 the rate appears to be underestimated. The Geant4 (geant4−10−07) simulations employing the Liege (QGSP INCLXX HP) model [54] combined with the native de-excitation code were the most successful in reproducing both 207 Hg and 208 Hg production, though with an apparent overestimation of the latter. Discrepancies with the Geant4 results may be a consequence of the necessity to scale from the simulation of a reduced density target, which was required to enable a feasible simulation time.
While 208 Hg production was not present in the FLUKA results, the simulations presented the possibility to investigate the 208 Pb(t,3p) 208 Hg channel. The results for tritium production (using FLUKA2021.0) were convoluted with cross-section data from TENDL17 [56] over a  MeV interval. This resulted in an in-target production rate of 110 atoms/µC for 208 Hg, significantly below the estimated rate of ≈56,000 atoms/µC calculated considering a 1% extraction and ionization efficiency. This suggests that 208 Pb(t,3p) 208 Hg reactions only contribute to a fraction of the observed 208 Hg yields. Based on this, we tentatively conclude that the observed 208 Hg production is the result of a combination of multiple reaction channels.

IV. EXPERIMENTAL RESULTS
A. Laser spectroscopy of mercury isotopes across N=126 Mean square charge radii, magnetic-dipole and electricquadrupole moments were studied via the measurement of isotope shifts and hfs in the 254-nm transition. Sample spectra are presented in Fig. 4(a), with the (substate weighted) centroids indicated with solid black lines.
Reference scans of 198 Hg were taken periodically to monitor the stability of the experimental setup, with a 10 h interval. Multiple measurements of the hfs of each isotope were taken (2× 202,203 Hg, 3× 206,207 Hg, 5× 208 Hg) and the fitting of the resulting spectra was cross-checked using multiple software packages: Origin 2016 with a chi-squared minimization [57] performed with a Levenberg-Marquardt algorithm [58,59], the SATLAS open source Python package [60] and a similar program written in ROOT [61]. The results are presented in Table I, together with literature data ( 202,203,206 Hg) for comparison. Relative mean square charge radii, electric dipole and magnetic quadrupole moments were extracted from the spectra by applying standard methods, these are summarized in Appendix A.
The nuclear spin of 207 Hg could not be determined unambiguously because the spectroscopic transition is between atomic states with electronic spins J=0 and J=1. I π = 9/2 + was assumed for the analysis of the 207 Hg measurements based on Refs. [48,64]. The extracted isotope shifts for 202,203 Hg are in good agreement with literature. The same is true for the neutron deficient isotopes that were re-measured during this experimental campaign [26,27]. The 500-MHz discrepancy between the δν 206,198 value of [65] and this work is discussed in [66]. The general agreement of the extracted δv A,198 and hyperfine a and b factors with the previously published literature values further validates the method of insource resonance ionization spectroscopy with a VADLIS ion source.

B. Mass spectrometry of 206−208 Hg
The masses of 206−208 Hg were measured, employing different techniques with respect to earlier experiments that are referenced in the AME2020 [67]. A short outline of the time-offlight methods that were used in our experiment can be found in Appendix B. A summary of the measured values is presented in Table II and a comparison with the AME2020 is shown in Figure 5.   [62] and [63]. The spin assignment of 207 Hg is discussed in the text. Statistical uncertainties are listed in parentheses and the systematic uncertainties related to F 254 and M are listed in curly brackets. The atomic mass of 206 Hg was previously deduced from α-decay measurements [68]. Our measurements of 206 Hg, using the Time-of-Flight Ion-Cyclotron-Resonance (ToF-ICR) technique, represent the first direct determination of the mass of this nucleus. For the ToF-ICR measurements in a Penning trap, two excitation schemes were employed: a single pulse excitation of 400 ms, as well as a Ramsey-type scheme [69], presented in Figure 4b, where two excitation pulses of 60 ms were applied, separated by a waiting time of 480 ms.
The masses of 207 Hg and 208 Hg have been determined previously by storage-ring measurements at GSI using Schottky mass spectrometry [70][71][72]. In the present experiment, ToF-ICR measurements were performed for 207 Hg with a singlepulse excitation of 1.2 s (see Figure 4c). For 208 Hg, the hyperfine structure was measured in five different laser scans, using ISOLTRAP's MR-ToF MS as mass separator and ion counter. Mass data were extracted from the scans by summing individual binned data per scan step into a single histogram, resulting in one histogram per scan. The ToF-distributions corresponding to singly charged ions of the isotope of interest, 208 Hg + , and to an on-line reference ion, 208 Pb + , were aligned as outlined in [73] and fitted by employing an unbinned maximumlikelihood estimation where the fit function was constructed as an Exponential-Gaussian-Hybrid (EGH) to account for the tails of the ToF-distributions towards longer flight times. The mass was extracted by calculating the average C ToF of the five scans, using the 208 Pb + present in the RIB and 133 Cs + from an off-line ion source as reference ions. The results we report are in general agreement with the literature and improve upon the precision.  Table I. The g-factors of the νg 9/2 iso- tones 209 Pb, 211 Po are plotted in Fig. 6 together with the Schmidt value, the energies of the first excited 2 + states (E(2 + ) −2 ) of the N=126 cores and g( 210 Bi) calculated from the measured magnetic moments of the 5 -, 7and 9isomeric states [74,75] using the additivity relation and assuming a pure [πh 9/2 ⊗ νg 9/2 ] configuration for these states. The inverse square of E(2 + ) is of particular relevance due to its approximate proportionality to the second-order perturbation theory correction [76,77]. The experimentally determined g-factors differ significantly from the Schmidt value g S chmidt =−0.425 and there is a noticeable Z-dependence. The deviation of the g-factor of the magic 209 Pb isotope was explained in non-relativistic [81,82] and relativistic [83] approaches by taking into account corrections for the meson exchange current and first and second order core polarization (CP). The configuration admixture contributing to the magnetic dipole moment in the first order of perturbation theory corresponds to a particle-hole excitation from an orbit j = l + 1/2 to its spin-orbit partner j = l-1/2 (CP1 correction [84]). In the second order of perturbation theory, the most important magnetic moment correction stems from the odd-particle coupling with the lowest 2 + excitation of the core (CP2 correction [76,77]). In the vicinity of the doubly magic 208 Pb, the most important Z-dependent CP1 correction corresponds to the proton (h −1 11/2 h 9/2 ) particle-hole excitation. A corresponding increase in the occupancy of the πh 9/2 orbital with the increase of Z, reduces the probability of proton (h −1 11/2 h 9/2 ) core excitations, thereby decreasing the magnitude of the CP1 correction. The opposite is apparent in Fig. 6, where the deviation from the Schmidt value increases between 209 82 Pb and 211 84 Po, thus CP1 corrections do not appear to be the dominant driver for the g-factor Z-dependence. Additionally, meson exchange corrections have been shown to TABLE II. Mass-measurement results for the mercury isotopes, given either as the ratio R of cyclotron frequencies from ToF-ICR or as the C ToF from the MR-ToF MS. The computed mass excess values M exc from this work are compared to the literature values found in [67]. Additionally, half-lives T 1/2 from [67] are given as well as the reference ions that were used to extract the mass values. The mass excesses are given in terms of energy divided by the square of the speed of light c 0 . have a weak A dependence in the vicinity of 208 Pb [77,85]. It could then be suggested that CP2 corrections are primarily responsible for the Z-dependence of the discrepancy with the Schmidt value.
The same mechanism (particle-quadrupole-vibration coupling) was used to explain the magnetic moment evolution in the vicinity of the magic numbers [76,77]. This explanation is supported by the apparent correspondence of the energies of the first excited 2 + states of the N=126 cores and the g-factors of the N=127 isotones in Fig. 6. Considering that E(2 + , 208 Pb)≈4.1 MeV, E(2 + , 206 Hg)≈1.1 MeV and E(2 + , 210 Po)≈1.2 MeV [80], the admixture of the (2 + ,νi 13/2 ) 9/2+ should be larger for 207 Hg and 211 Po than for 209 Pb, resulting in an increased CP2 correction. It is worth noting that this interpretation indicates the importance of particle-vibration coupling in the description of the groundstate properties of the odd-A nuclei in the vicinity of shell closure. This same mechanism proves to be decisive for explanation of the charge radii behavior in this region of the chart of the nuclides (see Sec. VI).

B. Quadrupole moment of 207 Hg
Quadrupole moments near the closed proton and neutron shells have a predominantly single-particle nature and are usually well described by the shell-model formula: where j is the spin of the odd particle, r 2 j is its mean square radius and e e f f is an effective charge. Taking r 2 j from Ref. [86], one obtains: e e f f ( 207 Hg; ν2g 9/2 )=2.4 (13) e. This neutron effective charge is noticeably larger than the value of the "universal" neutron e e f f = 0.95 e which describes fairly well the measured quadrupole moments for all closed-shell ±1 nuclear states in the vicinity of 208 Pb [87]. This points to the rapid increase of a quadrupole core polarization when moving away from the magic Z=82, similar to that observed for the proton e e f f when moving away from the magic N=126 [87].

VI. THEORETICAL INTERPRETATION
A. Theoretical framework and the details of the calculations A theoretical interpretation of experimental data is performed within the framework of covariant density functional theory [88] using the RHB computer code for spherical nuclei first employed in Ref. [9]. This code enables the blocking of selected single-particle orbitals and allows for fully self-consistent calculations of the ground and excited states in even-even and odd-A spherical nuclei. In the pairing channel, the RHB code employs a separable version of the Gogny pairing [89] with the pairing strength defined in Ref. [23].
Several restrictions/constraints are employed in the present paper. First, we consider only spherical nuclei (i.e. for which β 2 2 1/2 < 0.1 where β 2 2 1/2 is the mean-square deformation deduced from experimental δ r 2 using the droplet model (see Refs. [90,91]).) This restriction corresponds to N ≥ 116 and N ≥ 121 for lead and mercury isotopes, respectively, and it was already used in our earlier paper [9].
Second, two different procedures labeled as "LES" and "EGS" are used for the blocking in odd-A nuclei, and the results of the respective calculations are labeled by these abbreviations. In the LES (lowest in energy solution) procedure, the lowest in energy configuration is used, which is similar to all earlier calculations of OES in non-relativistic DFTs [4,17]. In the EGS (experimental ground state) procedure, the configuration with the spin and parity of the blocked state corresponding to those of the experimental ground state is employed, although it is not necessarily the lowest in energy. The need for this procedure is due to the following considerations. First, the nodal structure of the wavefunctions and neutron radius of the single-neutron orbital depend on its quantum numbers such as total and orbital angular momenta (see Refs. [3,5,92]). Thus, the occupation of different neutron single-particle states impacts the resulting charge radii (see Refs. [3,5] and the discussion of Fig. 7(a) below). Second, the structure of the experimental ground states in odd-A nuclei is reproduced globally only in approximately 40% of the nuclei in the non-relativistic and relativistic DFTs [93,94] If there is a mismatch in the structure of experimental and calculated ground states the impact of the blocked orbital on the physical observable of interest (charge radius, deformations, binding etc.) in odd-A is expected to deviate from the value observed in experiment. The consequences of this mismatch exist for the deformations of one-quasiparticle states (see Ref. [94]), odd-even staggerings in charge radii (see Ref. [9]) and binding energies (see Ref. [97]). Note that the latter is used to define pairing indicators. In such a situation it is safer to use the blocked solution with the structure corresponding to experimental ground state (even if it leads to an excited solution) for the description of charge radii since moderate shift in the energy of single-particle state has negligible effect on its neutron single-particle rms radius r 2 sp .

B. Charge radii and related indicators
The charge radii were calculated from the corresponding point proton radii as where r 2 stands for mean square radius of proton density distribution and the factor 0.64 accounts for the finite-size of the proton. Three indicators are commonly used to facilitate the quantitative comparison of the experimental results with those from theoretical calculations. The first is differential mean-square charge radius 2 δ r 2 N,N = r 2 (N) − r 2 (N ) = r 2 ch (N) − r 2 ch (N ) (3) where N is the neutron number of the reference nucleus. The second one is the ξ even indicator ξ even = δ r 2 128,126 δ r 2 126,124 , introduced in Ref. [8,66]. It provides a quantitative measure of the change of the slope (i.e. δ r 2 N,126 /δN) of differential charge radii as a function of neutron number N at N = 126. The applicability of ξ even is restricted by the limited availability of the data for the N=128 isotones because for 84 ≤ Z ≤ 88 they have half-lives (t 1/2 ) <300 µs [98] which limits the potential for laser spectroscopy measurements. Thus, in these cases the ξ even indicator is replaced by where N 0 is the lowest even neutron number at N > 126 with measured isotopic shift: N 0 = 132 for 84 88 Ra [105], and N 0 = 138 for 89 Ac [106]. For the 82 Pb and 83 Bi isotopes, the data for the nuclei with N = 124, 126, 128 were taken from [78] and [8], respectively. Note that for the Bi and Pb isotopes, for which the N = 128 data are available, 2 N 0 −126 δ r 2 N 0 ,126 ≈ δ r 2 128,126 (N 0 = 130 for Bi and N 0 = 132 for Pb) since differential radii increase nearly linearly for the N > 126 nuclei under study.
The third is the three-point indicator which quantifies OES in charge radii.
C. The kink in charge radii and its relation to underlying single-particle structure and pairing The differential charge radii of the Pb and Hg isotopes are shown in Fig. 7. One can see that all of the employed CEDFs generate a kink at N = 126 and that the best description is provided by the CEDF DD-ME2. Thus, it is important to understand which physical features determine the differences between the functionals. For that we look at the energies of the neutron single-particle states and their occupation probabilities.
The energies of neutron single-particle states obtained in the 208 Pb nucleus with the employed CEDFs are shown in Fig. 8. One can see close similarities in the predictions of the energies and relative positions of the 2 f 5/2 , 3p 3/2 and 3p 1/2 states occupied in the N ≤ 126 nuclei. In contrast, the differences are more pronounced for the 1i 11/2 and 2g 9/2 orbitals located above the N = 126 shell closure. In all functionals, the 1i 11/2 orbital is the lowest in energy 3 but the energy gap between these two orbitals strongly depends on the functional. It is smallest in DD-ME2, gets larger in NL3* and DD-PC1 and become extremely large in DD-MEδ. Because of this feature, and the fact that the kink in charge radii at N = 126 is defined by the interplay of the occupation of these two orbitals, we focus our discussion on the N > 126 nuclei.
Let us analyze the slope of differential radii defined as δ r 2 N,N /δN. Note that in such an analysis we consider only even-even nuclei. The results of the calculations without pairing indicate that this slope is almost the same for N < 126 4 and N > 126 nuclei when only the 2g 9/2 states are occupied above N = 126 (see Fig. 7(a) for Pb isotopes and Fig. 7 for Hg isotopes). As a consequence, there would be either no kink or a very small kink in the charge radii at N = 126. However, the situation drastically changes when only 1i 11/2 states are occupied above N = 126. This leads to a substantial increase of the δ r 2 N,N /δN slope and as a result to a creation of large kink in the charge radii at N = 126. This is related to the fact that the neutron 1i 11/2 orbital is the n = 1 orbital which overlaps more strongly with the majority of the proton orbitals than the n = 2 neutron 2g 9/2 orbital [5]. As a consequence, it provides a larger pull of the 1i 11/2 neutron states on proton orbitals via the symmetry energy. This is despite the fact that in 208 Pb the rms radius of the neutron 1i 11/2 orbital (for example, r ch = 6.4131 fm in DD-ME2) is smaller than that of the neutron 2g 9/2 one (for example, r ch = 7.0227 fm in DD-ME2) by ≈ 0.6 fm in all employed CEDFs.
The inclusion of pairing modifies the situation in the N > 126 nuclei in such a way that both of these orbitals become partially occupied (see Fig. 9(a)). Note that we consider a cumulative occupation probability v 2 state which provides information on the filling of a given j-subshell: v 2 state could take any value between 0 (unoccupied subshell) and 2 j + 1 (fully occupied j-subshell). Pairing also leads to a partial occupation of the single-particle states located above the N = 148 gap (see Fig. 8), but their occupation probabilities are relatively small because of the presence of this gap. Thus, for the sake of simplicity we focus our discussion on the interplay of the occupation of the 2g 9/2 and 1i 11/2 states and the consequences of this interplay on the δ r 2 N,N /δN slope.
The large energy gap of 1.27 MeV between the 2g 9/2 and 1i 11/2 states in the DD-MEδ functional (see Fig. 8) is responsible for a predominant occupation of the lower lying 1i 11/2 states (see Fig. 9(a)). Considering the N = 134 nucleus as an example, the eight neutrons outside of the N = 126 shell closure are located almost entirely in the 1i 11/2 subshell (v 2 1i 11/2 ≈ 7.1) with only a small portion occupying 2g 9/2 (v 2 2g 9/2 ≈ 0.5). This leads to a large δ r 2 N,N /δN slope (see Fig. 7). Note that this slope is the largest among the considered functionals and it is not far away from the one obtained in the calculations without pairing, when only the 1i 11/2 states are occupied in the nuclei with N > 126 (see Fig. 7). The reduction of the gap between the 2g 9/2 and 1i 11/2 orbitals in the DD-PC1, NL3* and, especially, DD-ME2 functionals (see Figs. 8 and 9(b)) is responsible for a decrease of the difference in the occupation of these orbitals (see Fig. 9). For these three functionals, the eight neutrons outside of the N = 126 shell closure in the N = 134 nucleus are still located predominantly in the 1i 11/2 subshell (v 2 1i 11/2 ≈ 4.7) but with a significant portion also found in 2g 9/2 (v 2 2g 9/2 ≈ 2.5). This leads a reduction of the δ r 2 N,N /δN slope as compared with the case of the DD-MEδ functional. One can see in Fig. 7 that the experimental δ r 2 N,N /δN slope in the N > 126 nuclei is reproduced with comparable accuracy by the DD-ME2, NL3* and DD- PC1 functionals. Note that with increasing neutron number N the single-particle states become more bound, but for a given functional, the relative energies of the 2g 9/2 and 1i 11/2 states change only slightly (see Fig. 9(b)). As a consequence, the occupation probabilities of the single-particle states behave as almost linear functions of neutron number (see Fig. 9(a)). The magnitude of the kink in charge radii at N = 126 is better quantified by the ξ even and ξ * even indicators defined in Eqs. (4) and (5), respectively. Experimentally available values of these indicators (both for even-even and odd nuclei) and calculated (only for even-even nuclei) are compared in Fig. 10. Mercury is the first element below Z=82 and only the second even-Z element for which ξ even is experimentally determined; the kink in the mercury charge radii is of a similar magnitude to that of lead. Only experimental ξ * even indicators are available for Z > 83 nuclei. Note that experimental ξ even and ξ * even indicators form a smooth curve which indicates that neither addition nor subtraction of proton(s) from the Z = 82 nuclei affects drastically a kink in charge radii at N = 126. Figure 10 clearly shows that the best reproduction of this trend is achieved by the CEDF DD-ME2. Other functionals (including the Fayans Fy(∆r) functional) deviate somewhat from experimental data.
This analysis clearly indicates that the evolution of charge radii with neutron and proton numbers is sensitive to the details of underlying singe-particle structure and the occupation probabilities of these states. Note that the latter depends on the type of pairing interaction employed in the calculations and its strength. It also indicates that the magnitude of the kink in charge radii at N = 126 depends on the relative balance of the occupation of the 2g 9/2 and 1i 11/2 orbitals. In general, the substantial occupation of the 1i 11/2 orbital and the kink can be obtained even when the 2g 9/2 orbital is located lower in energy than the 1i 11/2 orbital (but still within its vicinity).This, for example, takes place in the non-relativistic HFB (NR-HFB) calculations with the semirealistic M3Y-P6a interaction, the spin-orbit properties of which were modified [13] to improve the description of the charge radii of proton-magic nuclei [13][14][15]. However, these calculations underestimate the kink. A similar situation exists in the Skyrme DFT calculations with the SLy4mod functional presented in Ref. [5].

D. Odd-even staggering in charge radii
It was shown in Ref. [9] both in the RHB calculations with DD-ME2 and in the NR-HFB studies with semi-realistic M3Y-P6a interaction that OES in charge radii is best reproduced when the EGS procedure is applied in odd-A nuclei. In contrast, the experimental OES is significantly underestimated when the LES procedure is used in odd-A nuclei in the RHB framework for all nuclei under study and for N < 126 nuclei in the NR-HFB approach. The same behavior is observed also in the results of RHB calculations with CEDFs DD-MEδ, DD-PC1 and NL3*, the results of which are shown in Fig. 11. This figure shows also that there is some dependence of the magnitude of the ∆ r 2 (3) (N) values on the employed functional which comes from the differences in the energies of the single-particle states and their occupations (see Figs. 8 and 9).
Particle-vibration coupling (PVC) plays a critical role in the emergence of such significant OES in charge radii because it leads to the nucleonic configuration with the blocked state corresponding to the experimental ground state (see Ref. [9] for details). Let us illustrate that with the case of the N > 126 Pb nuclei. In the odd-A isotopes, the PVC coupling lowers the 2g 9/2 state below the 1i 11/2 state making it the ground state despite the fact that at the mean field level (as well as in eveneven nuclei) the energy of the 2g 9/2 state is higher than that of the 1i 11/2 state (see Fig. 5 of Ref. [95]). Although these results were obtained with NL3*, a comparable effect of PVC on relative energies of these states is expected in other functionals. However, the feasibility of such a scenario also crucially depends on the relative energies of the single-particle orbitals of interest at the mean field level. For example, if the energy gap between the 2g 9/2 and 1i 11/2 states is too large (like in the case of the DD-MEδ functional, see Figs. 8 and 9(b)), the impact of the PVC would most likely not be enough to make the 2g 9/2 state as a ground state. This functional however suffers from some significant deficiencies in the Z > 82 region which probably are the consequences of the overly large energy gap between these two states. For example, it does not predict octupole deformed actinides [109] and predicts fission barriers in superheavy nuclei which are too small to make them relatively stable [110]. One has to keep in mind, however, that  Eq. (5), using data from [8,9,62,78] and the references listed in the explanation of Eq. (5), while the calculated ones are from the charge radii defined in the present paper and in Ref. [23] which is publicly available at Ref. [107]. The ξ even (Pb) value for Fayans Fy(∆r) functional was extracted from the data presented in Ref. [18]).
the interaction DD-MEδ is different from the other interactions discussed here. It is the most microscopic one among considered functionals: only four parameters at the saturation density are fitted to finite nuclei and the full density dependence of the parameters is derived from ab-initio calculations.
On the contrary, the other interactions contain an additional 2 (NL3*), 4 (DD-ME2) or 6 (DD-PC1) phenomenological parameters for the fine-tuning of the density dependence.

E. Binding energies and ∆ (3) E indicators
Empirical approaches considering binding energies can help with understanding the interplay between nucleons. With exception of the DD-MEδ functional, there is a good agreement between the calculated and the experimental binding energies for both lead and mercury isotopes. E are reproduced reasonably well (see Fig. 12), especially when the LES are used in odd-A nuclei. The level of agreement is comparable with that provided by the Fayans Fy(∆r) functionals in Ref. [17]. Figure 12 also illustrates the challenges faced by all existing theories. Polarization effects in deformation/radial density distributions and pairing depend on blocked state in odd-A nuclei. Thus, on the one hand, the ∆ (3) E indicators are distorted by incorrect polarizations effects when a wrong (as compared with experiment) state is used for the ground state of odd-A nucleus. The related uncertainties in binding energies due to polarization effects in the pairing channel are on the level of 150 keV and those due to deformation/radial density distribution polarization effects are more difficult to estimate but are expected typically to be less than 100 keV. On the other hand, large theoretical uncertainties in the predictions of single-particle energies (see Refs. [94][95][96] and Fig. 8) reveal themselves when the ∆ (3) E indicators are defined using the EGS procedure in odd-A nuclei. This is especially pronounced in the calculations for the N > 126 nuclei with the DD-MEδ functional (see Fig. 12) which is characterized by a large energy gap between the 2g 9/2 and 1i 11/2 subshells (see Fig. 8). Comparing these uncertainties, it is safer to use the LES procedure in the definition of the ∆ (3) E indicators and their association with pairing indicators (see detailed discussion in Ref. [97]). Note that the calculations somewhat overestimate experimental ∆ (3) E indicators. The PVC provides additional binding (on the level of few 100 keV) to the ground states of odd-A nuclei. Thus, its inclusion into the calculations is expected to decrease the calculated ∆ (3) E and, as a consequence, to improve the description of experimental data.

VII. CONCLUSION
The kink in the δ r 2 systematics of the mercury isotopic chain has been analyzed considering a variety of dimensionless parameters and employed in the comparison of the results of a range of covariant energy density functionals. The results of mass measurements of 206−208 Hg are presented which improve upon the precision of previously measured values. The observed Z dependence of the g-factors of I=9/2, N=127 isotones is interpreted as the result of CP2 corrections. The magnitude of the extracted neutron effective charge for 207 Hg suggests a rapidly increasing quadrupole core polarization when moving away from the Z=82 proton shell.
The theoretical analysis of charge radii and related indicators in the Pb and Hg isotopic chain has been performed within the RHB framework using several CEDFs characterized by different single-particle properties. This analysis supports the conclusion that the kink at N = 126 in δ r 2 N,N originates from the occupation of the ν1i 11/2 orbital located above the N = 126 shell gap. The pairing effect does not play a critical role here since the kink is present also in the calculations without pairing. This is in contrast to non-relativistic Skyrme and Fayans functionals, in which the pairing becomes a dominant contributor to the kink and OES [18]. However, the pairing plays an important role in defining the magnitude of the kink which depends on the balance of the occupation of the 1i 11/2 and 2g 9/2 orbitals. This balance sensitively depends on the relative energies of these two orbitals. The DD-ME2 functional provides the best description of kinks at N = 126 not only in the Pb and Hg isotopes but also in all isotopic chains for which experimental data is available. A reasonable description of OES in charge radii has been achieved with all employed functionals. This confirms a new mechanism of OES suggested in Ref. [9] which is related to the staggering in the occupation of neutron orbitals between odd and even isotopes facilitated by PVC in odd-mass nuclei.

VIII. ACKNOWLEDGEMENTS
This project has received funding through the European Union's Seventh Framework Programme for Research Nuclear observables were extracted from the experimentally measured hfs through the application of standard methods [90]. The substate weighted centroid ν 0 and the hyperfine a and b parameters were extracted from the fitted spectra, where the shift ∆ν F of the state F of the hyperfine multiplet from ν 0 is given as where K=(F(F + 1) − I(I + 1) − J(J − 1)), I is the nuclear spin and J is the atomic angular momentum. The extraction of ν 0 enabled the calculation of δν A,A , the isotope shift between A the isotope under investigation and A a reference isotope. δ r 2 A,A , the change in the mean square charge radius of A with respect A was extracted as where for the spectroscopic transition F 254nm = −55.36 GHz fm −2 [62], K(Z = 80) = 0.931 (taking into account [111,112]) is a calculable correction factor and M is the mass shift factor representing the sum of the normal mass shift (M N MS ) and the specific mass shift (M(Z = 80) = (1 ± 0.5) · M N MS [62]).
The magnetic moments µ A were calculated as where the isomer 199m Hg (I = 13/2) was used as a reference nucleus with µ( 199m Hg) = −1.0147 (8) µ N [113] and a( 199m Hg) = −2298.3(2) MHz [113]. An additional correction for the hyperfine anomaly (HFA) is included for the µ values presented in Table I. Moskowitz and Lombardi [114] demonstrated that for mercury isotopes the "Bohr- [116] can be ignored. The following relation between the magnetic moments and the HFA, the so called Moskowitz-Lombardi (ML) rule was used, determining A 1 ∆ A 2 BW [115] as where α = 1 × 10 −2 and l is the orbital momentum of the last neutron. The ML rule was verified later by the microscopic theory [117]. The application of the ML rule to estimate the Bohr-Weisskopf correction of the magnetic moment of 207 Hg is justified based on the successful reproduction of the experimental HFA by this rule for neutron single particle states in mercury nuclei across a rather large range of masses [114]. For previously measured isotopes and isomers the maximum deviation of the experimental A 1 ∆ 199 BW from the ML calculation is equal to 2×10 −3 . Correspondingly, we conservatively estimated the error of the ML prediction for 207 ∆ 199 BW as 5×10 −3 . The uncertainty of this correction was estimated based on the omitted A 1 ∆ A 2 BR correction. It was shown in [118] that A 1 ∆ A 2 BR is proportional to δ r 2 A 1 ,A 2 . Thus 207 ∆ 199 BR for 207 Hg can be estimated by scaling the calculated 201 ∆ 199 BR [119]. It should be noted that 205 ∆ 203 BR for thallium isotopes, calculated in [119] by solving the one-electron Dirac equation, practically coincides with that calculated in [118] by the Dirac-Fock approach. Taking into account the independence of the A 1 ∆ A 2 BR correction on the details of the atomic calculations and the uncertainty of δ r 2 , we estimate the uncertainty of this correction as 10%.
The spectroscopic quadrupole moments Q A s were calculated using the relation where Q 201 Hg s =0.387(6) b (from [120]) and b201 Hg =−280.107(5) MHz [121]. The results are presented in Table I.

Appendix B: Mass spectrometry methods
The mass measurements were performed with the ISOLTRAP mass-spectrometer [122]. As described in parts above, the setup consists of four ion traps for beam preparation and mass measurements. First, a linear radio-frequency quadrupole trap was used to accumulate, cool and bunch the quasi-continuous radioactive ion beam delivered by ISOLDE [37]. Using a pulsed drift tube, the energy of the bunched beam is then reduced from the initial 30 keV to 3.2 keV.
Ions are captured by the Multi-Reflection Time-of-Flight Mass Spectrometer/Mass Separator (MR-ToF MS) [38] using the in-trap lift electrode [123], following this, they undergo a certain number of revolutions in-between the mirror electrodes of the device. Here, the time-of-flight is given as t = a √ m + b, where a and b are device-dependent parameters. Due to the mass-dependence of the trapped ions moving at the same kinetic energy, the different isobaric species separate in time-of-flight. In mass-spectrometry mode, the MR-ToF MS can be used to determine the mass of the ion of interest by using well known reference masses to account for the calibration parameters a and b. This is done by expressing the mass m using the so-called C T oF -value [40]: where m 1 , m 2 and t 1 , t 2 are the masses and time-of-flights of the two reference species, respectively, and C T oF = (2t − t 1 − t 2 )/ [2(t 1 − t 2 )]. An example of a time-of-flight spectrum is presented in Figs. 2(b) and 2(c).
In mass-separation mode, the MR-ToF MS selectively ejects the ions of interest towards the Penning traps located downstream of the setup. In the so-called preparation Penning trap, the ions are captured and cooled using a mass-selective buffer-gas method to improve beam emittance [124] and to further reduce contamination. Subsequently, the ions are ejected and recaptured in the precision Penning trap in which the high precision mass measurement is performed. In the present measurement, this is achieved by employing the Time-of-Flight Ion Cyclotron Resonance technique (ToF-ICR), which determines the cyclotron frequency of the trapped ions by scanning the frequency of an applied quadrupolar electric field [125]. This frequency can be written as ω c = qB m , where q is the electric charge of the ion, B is the magnetic field, and m is the ion mass. By performing cross-reference measurements with known mass calibrants, the mass m of the ion of interest can be extracted by comparing its cyclotron frequency with that of a well-known mass. Expressed as a frequency ratio the magnetic field and the charge cancel out. For the atomic mass, the electron mass is added to the measured ion mass. The ionization energy is negligible.
Expressing the ion masses in MR-ToF MS and ToF-ICR measurements via the C T oF -value and the frequency ratio R, respectively, facilitates an easy recalculation of the mass of interest in case one of the employed reference masses are measured more precisely in the future.