Refined nuclear magnetic dipole moment of rhenium: $^{185}$Re and $^{187}$Re

The refined values of the magnetic dipole moments of $^{185}$Re and $^{187}$Re nuclei are obtained. For this, we perform a combined relativistic coupled cluster and density functional theory calculation of the shielding constant for the ReO$_4^-$ anion. In this calculation, we explicitly include the effect of the finite nuclear magnetization distribution in the single-particle nuclear model using the Woods-Saxon potential for the valence nucleon. By combining the obtained value of the shielding constant $\sigma=4069(389)$~ppm with the available experimental nuclear magnetic resonance data we obtain the values: $\mu(^{185}{\rm Re})=3.1567(3)(12) \mu_N, \mu(^{187}{\rm Re})=3.1891(3)(12) \mu_N$, where the first uncertainty is the experimental one and the second is due to theory. The refined values of magnetic moments are in disagreement with the tabulated values, $\mu(^{185}{\rm Re})=3.1871(3) \mu_N, \mu(^{187}{\rm Re})=3.2197(3) \mu_N$, which were obtained using the shielding constant value calculated for the atomic cation Re$^{7+}$ rather than the molecular anion. The updated values of the nuclear magnetic moments resolve the disagreement between theoretical predictions of the hyperfine structure of H-like rhenium ions which were based on the tabulated magnetic moment values and available experimental measurements. Using these experimental data we also extract the value of the parameter of nuclear magnetization distribution introduced in [J. Chem. Phys. \textbf{153}, 114114 (2020)], which is required to predict hyperfine structure constants for rhenium compounds.

The refined values of the magnetic dipole moments of 185 Re and 187 Re nuclei are obtained. For this, we perform a combined relativistic coupled cluster and density functional theory calculation of the shielding constant for the ReO − 4 anion. In this calculation, we explicitly include the effect of the finite nuclear magnetization distribution in the single-particle nuclear model using the Woods-Saxon potential for the valence nucleon. By combining the obtained value of the shielding constant σ = 4069(389) ppm with the available experimental nuclear magnetic resonance data we obtain the values: µ( 185 Re) = 3.1567(3)(12)µN , µ( 187 Re) = 3.1891(3)(12)µN , where the first uncertainty is the experimental one and the second is due to theory. The refined values of magnetic moments are in disagreement with the tabulated values, µ( 185 Re) = 3.1871(3)µN , µ( 187 Re) = 3.2197(3)µN , which were obtained using the shielding constant value calculated for the atomic cation Re 7+ rather than the molecular anion. The updated values of the nuclear magnetic moments resolve the disagreement between theoretical predictions of the hyperfine structure of H-like rhenium ions which were based on the tabulated magnetic moment values and available experimental measurements. Using these experimental data we also extract the value of the parameter of nuclear magnetization distribution introduced in [J. Chem. Phys. 153, 114114 (2020)], which is required to predict hyperfine structure constants for rhenium compounds.

I. INTRODUCTION
Nuclear magnetic dipole moments are of wide interest for many physical problems. They can be used to test predictions of the nuclear theory. They are required as external parameters to predict the hyperfine structure (HFS) of neutral atoms, and molecules. Such data are required to probe the accuracy of calculated electronic wave functions, which are used for calculation of characteristics of symmetry-violation interactions in atoms [1][2][3][4][5][6] and molecules [7][8][9][10][11][12][13][14][15][16]. Such characteristics cannot be directly measured, but they are required to extract the value of the T,P-violating nuclear Schiff and magnetic quadrupole moments, the electron electric dipole moment and other similar effects from the experimental data [2,6,[17][18][19][20]. Magnetic dipole moments of stable isotopes can be combined with the experimental and theoretical data on hyperfine structure for stable and short-lived isotopes to obtain magnetic moments of short-lived isotopes [21][22][23][24][25][26][27][28]. Magnetic moments are used to predict hyperfine splittings in highly-charged ions, which can be used to test predictions of the bound state quantum electrodynamics [29].
Magnetic dipole moments of stable nuclei can be obtained from nuclear magnetic resonance (NMR) experiments on molecules, though there are suggestions to extract them from precise g-factor experiments on highly charged ions [30][31][32]. In molecular NMR experiments, one usually obtains so-called uncorrected values of the nuclear magnetic moment. To obtain the intrinsic value of the nuclear magnetic moment, one has to apply a correction on the shielding effect. It is induced by electrons surrounding the nuclei of interest in a given atom or a molecule. Accurate calculation of the shielding constant in molecules containing heavy atoms is rather complicated. Therefore, one often uses shielding corrections calculated for the corresponding atomic ions. Such approach can lead to serious errors [33][34][35].
In the present paper we study nuclear magnetic moments for two stable isotopes of rhenium, 185 Re and 187 Re, both having nuclear spin I = 2.5. Nuclear magnetic resonance experiments with the aqueous solution of the NaReO 4 molecule were carried out in 1951 [36]. The tabulated values of the nuclear magnetic moments of 185 Re and 187 Re [37] are based on those experimental data combined with the shielding constant calculated for the Re 7+ atomic ion [38][39][40]. In Ref. [41] it has been noted that such interpretation is not free from possible errors due to neglect of the chemical shift effect, i.e. contribution of molecular environment. The authors of Ref. [34] have used a combination of the nonrelativistic coupled cluster theory and the relativistic density function theory to calculate the molecular shielding constant. Here we perform a precise study of the shielding effect within the relativistic coupled cluster and relativistic density functional theories and show that the completely relativistic treatment allows one to significantly reduce the uncertainty of the shielding constant. In the present paper, we also explore the influence of the finite nuclear magnetization distribution effect in the singleparticle approximation on the shielding constant value. Using the refined value of the shielding constant we obtain the updated values of the nuclear magnetic moments of 185 Re and 187 Re.
The structure of the article is as follows. In Section II, we give a brief overview of the theory for calculating the shielding constant and the nuclear magnetization distribution effect. In Section III A, we discuss the shielding constant calculation scheme. Section III B provides an analysis of possible uncertainties. These two sections contain some technical details. In Section III C, we derive the nuclear magnetic moments from experimental NMR data and the theoretical value of the shielding constant obtained in sections III A and III B. Section IV analyzes the available experimental data on hyperfine splitting in H-like rhenium ions in various aspects.
The relativistic units (m = = c = 1) and the charge units α = e 2 /(4π) are used in this paper.

II. THEORY
One can use the following definition of the shielding tensor corresponding to the nucleus j in a given molecule: Here E is the energy of the system, µ j,a is the a'th component of the nuclear magnetic moment vector µ j of j'th nucleus, B b is b'th component of the uniform external magnetic field vector B. For the interpretation of the molecular NMR experiments, performed in a solution, we need the isotropic part σ of the shielding tensor, σ = 1/3 a σ a,a . From the nuclear magnetic resonance experiment, it is possible to obtain the uncorrected value µ uncorr. of the nuclear magnetic dipole moment value, i.e. the value, which is not corrected for the magnetic shielding. The intrinsic value of the magnetic moment µ can be obtained as: The interaction of electrons in a molecule with an external uniform magnetic field B can be described by the following term included in the Dirac-Coulomb Hamiltonian: where α are the Dirac matrices and r G = r − R G , R G is the gauge origin [42], i.e. the origin for the coordinate system that describes the electron radius-vector in this equation. In principle, the choice of R G can influence the results obtained in the modest basis sets (see below). In the point magnetic dipole approximation, the hyperfine interaction of an electron with the magnetic moment µ j of the jth nucleus can be written in the following way: where r j = r − R j , R j is the position of the nucleus j. Note, that the interaction (4) does not take into account the finite nuclear magnetization distribution effect.
In the theory of atomic hyperfine structure this effect is called the Bohr-Weisskopf (BW) effect [43][44][45]. One can use the following substitution to consider this effect [46][47][48]: Function F (r) takes into account the nuclear magnetization distribution inside the finite nucleus. In the point magnetic dipole moment approximation F (r) = 1. In the finite distribution case F (r) can significantly differ from 1 inside the nucleus. Expressions for different models can be found in Refs. [46,[48][49][50][51]. In the simplest uniformly magnetized ball model function F (r) = (r/r n ) 3 for r inside the sphere of radius r n = 5/3r c (r c is the root-mean-square charge radius) and is equal to 1 outside [49]. In studies of neutral atoms, this model is most widely used to calculate the BW correction [25,[52][53][54][55].
In the present paper, we mainly use the model which implies that magnetization can be ascribed to the singleparticle structure of the nucleus. In this model function F (r) is given by [46]: for I = L + 1/2, and for I = L − 1/2. Here µ N is the nuclear magneton, m p is the proton mass, I is the nuclear spin, |u(r)| 2 is the density of the valence nucleon, φ SO is the radial part of the spin-orbit interaction V SO = φ SO σ σ σ · l l l, l l l is the angular moment operator and σ σ σ is the vector of Pauli matrices. In the Woods-Saxon (WS) model of the nucleus, the wave function of the valence nucleon is determined as a solution of the Schrödinger equation with the WS potential. A detailed description of the implementation and parameters of the potential can be found in Ref. [28] and references therein. For the valence proton we set g L = 1, for the valence neutron g L = 0. Parameter g S is obtained from the following equations: for I = L + 1/2, and µ µ N = − I 2(I + 1) g S + I(2I + 3) 2(I + 1) − 2I + 1 4(I + 1) m p φ SO r 2 g L (9) for I = L − 1/2. In the simple single-particle model with the uniform distribution of the valence nucleon, the density of the valence nucleon |u(r)| 2 is a constant inside the nucleus volume and there is no spin-orbit term in this model [56].
In the one-electron case, the tensor (1) can be calculated using the sum-over-states method corresponding to the second-order perturbation theory with perturbations (3) and (4): where |0 is the unperturbed one-particle state of interest, |n is the unoccupied n'th unperturbed state (orbital) and h.c. is the Hermitian conjugate. For the case of the four-component Dirac theory, the summation should include both positive energy and negative energy states |n [57]. The part of the sum associated with the positive energy states is called "paramagnetic" term. The part associated with the negative energy states is called "diamagnetic term" [57]. For the cases of the Dirac-Hartree-Fock (DHF) and density functional theory (DFT) manyelectron methods, one can use the response technique to calculate both terms [57][58][59][60]. The result of the application of this technique is equivalent to (and derived from) the analytical calculation of the DHF/DFT energy derivative (1). In the present paper we have used the implementation of the method within the dirac [60,61] code.
In calculations of the shielding constant, we have used the following Gaussian-type basis sets to describe electronic wave functions. The first one corresponds to the uncontracted Dyall's AE4Z [62,63] basis set for all atoms and will be called QZQZ below. This basis set contains [34s 30p 19d 14f 10g 5h 1i] primitive Gaussian functions for Re and [18s 10p 5d 3f 1g] functions for each oxygen. The second one, TZTZ, corresponds to the uncontracted AE3Z [62,63] basis set on rhenium and contracted aug-cc-pVTZ [64,65] on oxygen. This basis set contains [30s 24p 15d 11f 5g 1h] functions for Re and [5s 4p 3d 2f ]/([11s 6p 3d 2f ) for each oxygen, where in the [...]-brackets the numbers of contracted functions are given and in the (...)-brackets the corresponding numbers of primitive functions are given (e.g. each of five contracted s−type functions of oxygen is a linear combination of 11 primitive functions). We have also used the DZDZ basis set which corresponds to the uncontracted Dyall's AE2Z [62,63] basis set on rhenium and aug-cc-pVDZ [64,65]  Formally, the interaction of a molecule with an external uniform magnetic field should not depend on the choice of the origin R G in Eq. (3). But for finite-size basis sets there may be some dependence [42,58,59] which can affect the shielding constant value. To minimize such a dependence one can use the London atomic orbitals (LAOs) method, developed at the four-component DFT level in Refs. [58,59]. In this approach, basis functions are replaced by the so-called London atomic orbitals which are obtained from the original basis functions by applying a magnetic field-dependent factor [42,58,59]. This corresponds to the transformation of the wave function due to the gauge transformation of the vector potential in Eq. (3). The use of London orbitals guarantees the gauge-origin invariance of results in a finite basis approximation [42,58,59]. Even for usual basis set, the gauge-origin problem should decrease with the basis set size increase. In the present case, we are interested in the shielding constant for the rhenium nucleus. Therefore, it is natural to place the origin at this nucleus. According to our DFT estimates, the values of the shielding constant calculated for the QZQZ basis set (i) with such choice of the origin and employing usual basis functions or (ii) within the LAOs technique coincide within 7 ppm. This value is negligible in comparison with the total uncertainty of the present calculation (see below).
Geometry structure parameters of the ReO − 4 anion have been optimized using the four-component density functional theory with the Perdew-Burke-Ernzerhof, PBE0, functional [66] and using the TZTZ basis set. No solvent effects were considered at this stage. The optimized value of the Re-O bond length in the ReO − 4 anion with the regular tetrahedral symmetry was found to be 1.723Å. This value is in good agreement (within 0.003Å) with the study [34].
Relativistic four-component calculations were performed within the locally modified dirac15 [60,61] code. High-order correlation effects have been calculated using the mrcc code [67]. The code for calculating the BW matrix elements in the WS model has been developed in Ref. [28] for atoms and generalized to the molecular case in the present paper. Taking into account that the action of the corresponding operator is localized inside the nucleus we have neglected the contribution of basis functions centered on oxygen atoms in the present implementation. In the molecular electronic structure calculations, the Gaussian nuclear charge distribution model [68] has been used. We have used the following scheme to calculate the shielding constant for ReO − 4 and its contributions (see Table I). The diamagnetic part has been calculated at the four-component PBE0 method [66]. As in previous studies [33,69] we have found that this contribution is almost independent of the choice of the functional or method used. For example, the values calculated within the Dirac-Hartree-Fock (7630.8 ppm) and PBE0 (7633.3 ppm) response theories coincide within a few ppm. Moreover, the same value within a few ppm can be obtained even using the uncoupled Dirac-Hartree-Fock approach (7633.2 ppm), i.e. simple orbital perturbation theory. The latter corresponds to calculation using Eq. (10) with an additional summation over all occupied molecular orbitals, |0 , which are included in the Slater determinant, and the sum over n in Eq. (10) is limited to the negative energy orbitals.
The most challenging part of the problem is the calculation of the paramagnetic contribution to the shielding constant which is strongly affected by both correlation and relativistic effects. The authors of Ref. [34] have used the nonrelativistic coupled cluster theory combined with the relativistic correction calculated within the DFT approach to calculate the shielding constant for ReO − 4 . As it has been analysed in Ref. [34] the dominant source of the uncertainty of such an approach is "the systematic error of correlation and relativistic effects nonadditivity". To avoid such an error, we avoided the use of non-relativistic theory at all and used the four-component relativistic coupled clusters theory as the main approach for calculating the paramagnetic contribution to the shielding constant. We have directly calculated the mixed derivative (1) within the standard numerical finite-difference technique. For this we needed the (numerical) dependence of the energy E on the magnetic moment and external magnetic field values. To calculate this dependence we have added perturbations (3) and (4) with required (small) values of the nuclear magnetic moment and the external uniform magnetic field amplitude to the molecular relativistic Hamiltonian and solved coupled cluster equations with these perturbed Hamiltonian to obtain perturbed values of the energy E. This electronic correlation calculation has been performed within the TZTZ basis set using the relativistic coupled cluster with single, double and perturbative triple cluster amplitudes method, CCSD(T) [70]. All 108 electrons of ReO − 4 were included in the correlation calculation and no virtual energy cutoff has been applied. In Table I we separate the CCSD(T) value into the CCSD value (line "TZTZ/108e-CCSD") and a contribution of perturbative triple cluster amplitudes (line "TZTZ/108e-CCSD(T) -108e-CCSD"). Thus the "TZTZ/108e-CCSD(T) -108e-CCSD" line gives the difference of the shielding constants calculated within the CCSD(T) and CCSD methods (in both cases all electrons were correlated and the TZTZ basis set has been used).
To explore even higher-order correlation effects, we have performed correlation calculations within the relativistic coupled cluster with single, double and iterative triple cluster amplitudes method, CCSDT, and compared it with the CCSD(T) one. Due to extremely high complexity of the CCSDT approach (e.g. in the present calculation our cluster operator included 2.3 billions cluster amplitudes), we have included 24 valence electrons of ReO − 4 in these two calculations and employed the DZDZ basis set. One should note that the contribution of perturbative triple cluster amplitudes in the DZDZ/24e-CCSD(T) calculation reproduces such contribution obtained in the main calculation (350 ppm, see line "TZTZ/108e-CCSD(T) -108e-CCSD" in Table I) within 80%. As one can see, the difference between CCSDT and CCSD(T) results (-81 ppm) is rather small (see the "DZDZ/24e-CCSDT -24e-CCSD(T)" line in Table I).
To take into account the effect of the extended basis set with respect to the main TZTZ one, we have calculated basis set correction within the relativistic PBE0 approach. In this calculation, we have increased the basis set up to the QZQZ one, i.e. we have calculated the difference of the paramagnetic contributions to shielding constants calculated within the QZQZ and TZTZ basis set using the PBE0 approach. One can see from Table I that this correction is small (-10 ppm). Note, that there is no guarantee, that DFT can reasonably take into account basis set correction accurately. Therefore, we have estimated the influence of the basis set size increase from the DZDZ basis set to the TZTZ one on the paramagnetic part of the shielding constant. Here we studied how DFT (PBE0) can reproduce this effect, calculated within the wave function-based relativistic CCSD(T) method. Obtained corrections are: -30 ppm within DFT vs. -71 ppm within CCSD(T). As expected, DFT underestimated the effect of the basis set size increase (by a factor of 2.4). Thus, the mentioned correction (-10 ppm) in Table I can be underestimated. We take this fact into account in the uncertainty estimation below.
The effect of solvent has been extensively analyzed in Ref. [34]. It seems that the "δ 4 " scheme used in Ref. [34] is the most elaborate study of this effect at present. It explicitly takes into account the effect of the first solvation shell and approximately takes into account the influence of the solution on the shielding constant under consideration (within the polarizable continuum model) at the DFT level [34]. Therefore, we include this contribution in our final value.
The Gaunt interaction contribution has been calculated at the relativistic DFT (PBE0) level using the TZTZ basis set, i.e. it has been calculated as a difference between the shielding constant values obtained within the relativistic DFT method with inclusion and without inclusion of the Gaunt interaction into the electronic Hamiltonian.
In the present paper, we have studied the influence of the finite nuclear magnetization distribution on the shielding constant. For this, we have used the substitution given by Eq. (5) in the hyperfine interaction operator (4) and have used the single-particle WS model, described above. Calculation of the paramagnetic part has been performed at the relativistic CCSD(T) level using the DZDZ basis set, while the diamagnetic part has been calculated at the uncoupled Dirac-Hartree-Fock level (this method is described above). The latter contribution, 8 ppm, to the considered correction was found to be much smaller than the paramagnetic part, −81 ppm, of this correction. As one can see from Table I the finite nuclear magnetization distribution effect (-73 ppm or 1.8% of the total σ) is more important than the solvent effect for the system under consideration. For both isotopes considered, nuclear magnetization distribution effects were found to be almost identical and they are not distinguished in Table I. We did not find previous attempts to take into account the influence of the finite nuclear magnetization distribution effect on the shielding constants in many-electron molecules, although such studies have been performed for H-like ions [71,72].

B. Shielding constant uncertainty estimation
The total uncertainty (δ) of the present calculation of shielding constant can be estimated as a square root of the sum of squares of electronic correlation calculation uncertainty (δ el.corr. ), basis set uncertainty (δ BS ), Gaunt interaction effects inclusion uncertainty (δ G ), finite nuclear magnetization distribution uncertainty (δ FMD ), quantum electrodynamics (δ QED ), uncertainty due to geometry structure uncertainty (δ geom ) and solvent effects uncertainty (δ sol ): δ = δ 2 el.corr. + δ 2 FMD + δ 2 G + δ 2 BS + δ 2 sol + δ 2 QED + δ 2 geom . (11) Uncertainty of the calculated electronic correlation effect δ el.corr. can be estimated as the contribution of the perturbative triple cluster amplitudes given in Table I in the line "TZTZ/108e-CCSD(T) -108e-CCSD", i.e. δ el.corr. =350 ppm. For a more detailed analysis of this contribution, we also calculated it within the smaller DZDZ basis set, i.e. calculated the difference between shielding constants calculated using the DZDZ/108e-CCSD(T) and DZDZ/108e-CCSD approaches. The calculated value, 355 ppm, reproduces with high accuracy the contribution of perturbative triple clusters amplitudes within the TZTZ basis set given above, 350 ppm. This suggests a small contribution from the interference effect between the high-order correlation effects defined by perturbative triple amplitudes and the size of the basis set (the uncertainty associated with the size of the basis set itself is discussed below). Note that the indicated uncertainty δ el.corr. seems to be quite conservative, given that the estimate of the contribution of higher-order correlation effects (outside the considered CCSD(T) approximation) calculated within the framework of the CCSDT approach (see "DZDZ/24e-CCSDT -24e-CCSD(T)" line in Table I and description in section I) is several times smaller than δ el.corr. .
The uncertainty δ FMD of the calculated value of the finite nuclear magnetization distribution contribution (-73 ppm) can be estimated as about 30% of the value of this contribution (see also the analysis of the HFS data for Hlike Re below): δ FMD = 23 ppm. This value is obtained by comparing contributions of the finite nuclear magnetization distribution to the shielding constant calculated in the WS model, -73 ppm, and in the uniformly magnetized ball model, -96 ppm. The latter value has been calculated using the same approach as that employed for the first one (described in the previous section) with only replacement of the nuclear magnetization distribution function F (r).
As a measure of the uncertainties of the Gaunt interaction and solvent effects calculated within one method (DFT), we have used corresponding values of these effects given in Table I, i.e. δ G = 15 ppm, δ sol = 25 ppm. This means that the (conservative) uncertainties of these corrections are suggested to be 100%. A more accurate calculation of these effects and their uncertainties will be required when the remaining uncertainties are made smaller (e.g. the uncertainty δ el.corr. considered above is more than an order of magnitude bigger than the effects under consideration).
Let us estimate the uncertainty due to the basis set incompleteness δ BS . For this we can compare shielding constant calculated within the best employed QZQZ basis set and a smaller one, TZTZ basis set. These differences for the diamagnetic contributions, paramagnetic contribution and the total shielding constant values calculated within the QZQZ and TZTZ basis sets using the DFT (PBE0) method are +38 ppm, −10 ppm and +28 ppm, respectively. Thus, we can suggest that δ BS is about 28 ppm. However, as we have mentioned above, the DFT approach can underestimate the effect of the basis set size increase for the paramagnetic contribution by a factor of 2.4. This will lead to the estimation for the uncertainty of the basis set incompleteness as 38 − 10 · 2.4 =14 ppm. For the conservative estimate we choose the largest of these two possible values, i.e. we set δ BS = 28 ppm.
The geometry parameters, i.e. the Re-O distances in the ReO − 4 anion have been optimized at the relativistic PBE0 level (see above). To check the uncertainty of the optimized geometry we have also performed geometry optimization within another popular functional -B3LYP [73]. Within this approach the optimized Re-O distance was found to be 1.736Å, i.e. the estimation for the geometry parameters uncertainty of ReO − 4 is about 0.013Å. According to our DFT-based estimation, this uncertainty in the geometry structure parameters leads to contribution to the uncertainty of the shielding constant of about δ geom = 158 ppm. This is about six times bigger than the solvent effect contribution (and is about an order of magnitude bigger than the Gaunt interaction effect).
In a recent paper [74] it was estimated that contribution of quantum electrodynamics effects to the shielding constant is about 0.5% for such many-electon atoms as astatine (atomic number 85). For H-like ions, ab-initio calculations are available [71,72]. According to these Refs., the QED contribution to the shielding constant of H-like Bi (atomic number 83) is 0.7% and generally increases with atomic number. In the present work, we do not take into account QED effects but according to the notes above suggest that their contribution can be about 1% of the total value of the shielding constant and include this value in the uncertainty, i.e δ QED =41 ppm.
The final value of the uncertainty of the shielding constant is dominated by the correlation contribution uncertainty δ el.corr. . Substituting all estimated values of the uncertainties in Eq. (11) we obtain the total uncertainty value: δ=389 ppm, i.e. the final value of the shielding constant is σ =4069(389) ppm. It is in reasonable agreement with the previous calculation 3698(927) ppm [34] but has reduced uncertainty due to the use of the relativistic coupled cluster approach for the most challenging part of the calculation. Note also, that in the present paper we consider more sources of the uncertainties. To obtain the nuclear magnetic dipole moments of 185 Re and 187 Re we need corresponding uncorrected values of these moments, i.e. magnetic dipole moments which are not corrected for the magnetic shielding, see Eq.
(2). The nuclear magnetic resonance experiment on the ReO − 4 anion has been performed in the aqueous solution of NaReO 4 with the magnetic dipole moment of 23 Na as the reference [36]. In the experiment the two resonances of 185 Re and 187 Re were located near a frequency of 6.4 Mc in an external field of 6700 gauss. The following results for the resonance frequencies ν were obtained in the experiment [36]: ν( 185 Re)/ν( 23 Na) = 0.85114(9), (12) ν( 187 Re)/ν( 23 Na) = 0.85987 (9).
As in Ref. [75] we use the uncorrected NMR value µ uncorr. ( 23 Na) from Ref. [76] to obtain the uncorrected values of the nuclear magnetic moments of 185 Re and 187 Re [75] according to the relation [75] µ uncorr. (Re) = µ uncorr. (Na)(ν(Re)/ν(Na))(I Re /I Na ) (where the 23 Na nuclear spin I23 Na =1.5): Using these uncorrected values and our theoretical shielding constant value σ we obtain the final values of the nuclear magnetic moments according to Eq. (2): Here the first uncertainty is due to the experiment and the second is due to the present theory. The obtained value of the shielding constant for ReO − 4 molecular anion is about three times smaller than the shielding constant for the Re 7+ atomic cation [38][39][40] which was used in some of the previous interpretations of the molecular NMR data [37,75,77]. It means that in previous studies the uncertainty of the shielding correction used for interpretation of the molecular NMR experiment has been substantially underestimated [75].

A. Nuclear magnetic moments from HFS data for H-like Re
Rhenium is one of several elements for which measurements of hyperfine splitting for H-like (Re 74+ ) ion was carried out [41]. In principle, it is possible to extract the magnetic moment value from these HFS data if the values of BW and QED contributions are known [78]. For this one can use the following expression for the hyperfine structure constant A [79]: where A (0) is the HFS constant calculated in the point magnetic dipole approximation, A BW is the Bohr-Weisskopf contribution to the HFS constant, ε is the relative Bohr-Weisskopf correction, and A QED is the QED contribution. Constants A (0) , A BW and A QED are proportional to the nuclear magnetic moment µ = g I Iµ N , where g I is the g-factor of nucleus with spin I. A (0) and A QED can be accurately calculated [47,80] and A BW can be estimated within some nuclear magnetization distribution model. In Ref. [28] we have calculated the relative BW correction ε for H-like 185 Re. Here we have also estimated BW effect for the 187 Re isotope using four different nuclear magnetization distribution models (see Table II): uniformly magnetized ball model (Ball), singleparticle model with the uniform distribution of the valence nucleon (UD), WS model with and without spinorbit interaction in Schrödinger equation for valence nucleon (see the Theory section). The obtained values for both considered isotopes are almost identical, since these nuclei have similar single-particle structures and charge radii. Using the experimental values of A [41], the values of BW corrections ε calculated within different nuclear magnetization distribution models and given in Table II, the values of the ratio A (0) /g I I = 0.2926(3) eV calculated in [47,56] [81] and QED effect calculated in Refs. [47,80] A . The values of µ deduced in such a way using different models of nuclear magnetization distribution are given in Table III. In such approach the main uncertainty is due to the BW effect, as it is hard to reliably treat many-body nuclear structure effects. We suppose that this uncertainty can be estimated by comparing different nuclear magnetization distribution models given in Table II. Using this approach, the relative Bohr-Weisskopf correction can be estimated as ε =1.32(37)% for both isotopes [82]. It corresponds to the following value of the magnetic moments derived from the experimental HFS data [41] for H-like ions: Here the first uncertainty is due to the experimental determination of A [41], the second one is due to uncertainties of A (0) (the uncertainty of A QED is negligible) and the third one corresponds to the uncertainty of the calculated BW effect. These values are in agreement with Ref. [78] and in agreement with the values (16) and (17) derived from the NMR data above but have an order  Tables I and II). However, the influence of the BW effect on the final value of the magnetic moment, extracted from the NMR data is much smaller than in the case of H-like HFS data. In the first case, the finite nuclear magnetization distribution effect is a correction to the shielding effect, which is 4069 ppm, i.e. only 0.4% itself (see above). The uncertainty of the molecular electronic structure shielding constant calculation can be controlled, see section III B. In the second case, the BW effect directly contributes to the magnetic moment, and an estimation of its uncertainty is complicated due to the absence of direct many-body calculations of this effect for the rhenium nucleus. It is possible to employ the obtained results also as follows. If one uses the tabulated values [37] [47,83,84] and experimental values [41].
B. BW effect from HFS data for H-like ions and magnetic moments from molecular NMR data As mentioned in the Introduction, accurate theoretical prediction of the hyperfine structure of atoms and molecules can be used to probe the accuracy of the electronic wave function. However, such predictions depend on the nuclear magnetic dipole moment value and the function of the nuclear magnetization distribution F (r) in Eq. (5). If both of these components are accurately known, then one can predict the hyperfine structure of the compound or ion under consideration. It was shown in Ref. [15] that for many-electron heavy-atom molecules and heavy atoms to a good approximation it is possible to factorize the Bohr-Weisskopf contribution to the hyperfine structure constant into a pure electronic part and just one universal numerical parameter, which depends on the nuclear magnetization distribution, see Eq. (29) in Ref. [15] (recently a related approach has been considered for atoms in s and p 1/2 states in Ref. [85]). The latter parameter B s [15] is proportional to the BW contribution A BW to the hyperfine structure constant in Eq. (18) for the H-like ion in the ground electronic state and can be calculated as B s = A BW /2g I in this case; actually, the constant of interest is the product B S g I = A BW /2 [15]. Using the experimental values of the H-like rhenium HFS constants [41], QED corrections from Refs. [47,80] and the values of the nuclear magnetic moments (16) and (17) refined in the present paper above, we obtain the following values of the BW contribution to the HFS constants, A BW("exp ′′ ) , according to Eq. (18): A BW("exp ′′ ) ( 185 Re) = 0.0124 (6)(8)(4) A BW("exp ′′ ) ( 187 Re) = 0.0130(6)(9)(4) eV.
Here the first uncertainty is due to the experimental HFS data for H-like ions [41], the second one is due to uncertainty of A (0) and the third one corresponds to nuclear magnetic moment uncertainties in (16) and (17). From these values one can obtain for the product B S g I = A BW /2: B S g I ( 185 Re) = 0.0062(3)(4)(2) eV, B S g I ( 187 Re) = 0.0065(3)(5)(2) eV. Here the uncertainties correspond to the uncertainties in Eqs. (21) and (22) above.

V. CONCLUSION
In the present paper, we have obtained refined values of the magnetic moments of 185 Re and 187 Re nuclei. For this, we have calculated the shielding constant for the ReO − 4 anion using a combination of the relativistic coupled cluster and relativistic density functional theories. We have studied the influence of the finite nuclear magnetization distribution effect on the shielding constant value. Such effect is usually omitted in molecular calculations. However, according to our study, this effect can be more important than the solvent effect which is often estimated. Updated values of the nuclear magnetic moments resolve the disagreement between theoretical predictions [47,83,84] and experimental values [41] for the hyperfine splittings of H-like rhenium ions. In addition to the nuclear magnetic moment values, we have also used H-like data for rhenium HFS constants to extract the universal parameter [15] of the nuclear magnetization distribution. The values of the nuclear magnetic moment and this parameter are necessary ingredients for the theoretical prediction of HFS constants in different rhenium ions and compounds.

ACKNOWLEDGMENTS
We are grateful to V.M. Shabaev for useful discussions. Electronic structure calculations have been carried out using computing resources of the federal collective usage center Complex for Simulation and Data Processing for Mega-science Facilities at National Research Centre "Kurchatov Institute", http://ckp.nrcki.ru/, and partly using the computing resources of the quantum chemistry laboratory.
Molecular electronic structure calculations performed at NRC "Kurchatov Institute" -PNPI have been supported by the Russian Science Foundation Grant No. 19-72-10019. Calculations of nucleon wave function performed at SPbSU were supported by the foundation for the advancement of theoretical physics and mathematics "BASIS" grant according to Project No. 21-1-2-47-1.