Stefan Knecht
email: knecht@ifk.sdu.dk

Trond Saue
email: trond.saue@irsamc.ups-tlse.fr

Nuclear size effects in rotational spectra: a tale with a twist

We report a 4-component relativistic benchmark study of the isotopic field shift in the rotational spectrum of three diatomic molecules: TlI, PbTe and PtSi. A central quantity in the theory is the derivative with respect to internuclear distance of an effective electron density associated with a given nucleus, calculated at the equilibrium distance. The effective density, which is related to the mean electron density within the nuclear volume, is usually replaced with the contact density, that is, the electron density at the origin of the nucleus. Our computational study shows that for the chosen systems this induces errors on the order of 10%, which is not acceptable for high-precision work. On the other hand, the systematic nature of the error suggests that it can be handled by an atom-specific correction factor. Our calibration study reveals that relativistic effects increase the contact density gradient by about an order of magnitude, and that the proper transformation of the associated property operator is mandatory in 1-and 2-component relativistic calculations.

Our results show very good agreement with the experimental data presented by Schlembach and Tiemann [Chem. Phys. 68 (1982) 21], but disagree completely with the revised results given by the same group in a later paper [Chem. Phys., 93 (1985) 349]. We have carefully re-derived the relevant formulas and can not see that the rescaling of results is justified. Curiously previous DFT calculations agree quite well with the revised results for TlI and PbTe, but we demonstrate that this is because the authors inadvertently employed a non-relativistic Hamiltonian, which by chance induces an error of the same magnitude as the suggested scaling. For the PtSi molecule our results for the correction term due to nuclear volume disagree with experiment by a factor five, and we recommend a re-examination of the experimental data.

I. INTRODUCTION

The incompatibility of the α-scattering experiments of Geiger and Marsden [START_REF] Geiger | On a diffuse reflection of the α-particles[END_REF] with the multiple scattering predicted by Thomson's "plum pudding" atomic model [START_REF] Thomson | On the structure of the atom: an investigation of the stability and periods of oscillation of a number of corpuscles arranged at equal intervals around the circumference of a circle; with application of the results to the theory of atomic structure[END_REF] led [START_REF] Rutherford | The Scattering of α and β particles by matter and the structure of the atom[END_REF] to propose that all positive charge was concentrated in a nucleus at the center of a homogeneously negatively charged sphere [START_REF] Rutherford | The Scattering of α and β particles by matter and the structure of the atom[END_REF][START_REF] Heilbron | The scattering of α and β particles and Rutherford's Atom[END_REF]. Rutherford concluded that the size of the nucleus had to be very much smaller than the atomic radius, on the order of 10 -15 m or less.

Since the nuclear size is indeed negligible for the description of most phenomena on an atomic scale, an overwhelming majority of theoretical studies of molecules and solids, starting with the celebrated scattering formula of Rutherford, treat nuclei as point charges. A notable exception are relativistic molecular electronic structure calculations based on the finite basis approximation, where the point nucleus assumption leads to basis set convergence problems due to the weak singularities displayed by the solutions of the Dirac Hamiltonian at the nuclei. To overcome this problem it was early realized [START_REF] Ishikawa | On the use of an extended nucleus in Dirac-Fock Gaussian basis set calculations[END_REF][START_REF] Visser | The use of Gaussian nuclear charge distributions for the calculation of relativistic electronic wavefunctions using basis set expansions[END_REF] that the introduction of a finite (Gaussian) nuclear charge distribution not only curbs the singularity but also leads to Gaussian type solutions near the origin, which greatly alleviates the integral evaluation in relativistic wave function calculations using basis set expansions.

However, many phenomena known from atomic and molecular spectroscopy as well as isotope chemistry even lack a reasonable explanation if the finite size of a nucleus is not taken into account. A notable example is the isotope shift in the electronic spectra of atoms and molecules [START_REF] Van Vleck | On the isotope corrections in molecular spectra[END_REF][START_REF] Heilig | Changes in mean-square nuclear charge radii from optical isotope shifts[END_REF][START_REF] King | Isotope Shifts in Atomic Spectra[END_REF] which is well understood as an interplay between i) the mass shift, scaling approximately with ∆M/M 2 (M being the mass), and ii) the nuclear field shift, arising from the difference in size and shape of the extended nucleus of each isotope. The mass shift term typically becomes less significant with increasing nuclear charge values Z due to its scaling properties whereas the nuclear field shift is known to be the main contributor to the total isotope shift for Z ≥ 58 [START_REF] Heilig | Changes in mean-square nuclear charge radii from optical isotope shifts[END_REF]. Following similar lines of thoughts Bigeleisen [START_REF] Bigeleisen | Temperature dependence of the isotope chemistry of the heavy elements[END_REF][START_REF] Bigeleisen | Nuclear size and shape effects in chemical reactions. Isotope chemistry of the heavy elements[END_REF] developed in the mid 1990s a general theoretical model which was able to explain isotopic anomalies in chemical reactions, as for example observed in U(IV)-U(VI) exchange reaction experiments [START_REF] Fujii | An anomalous isotope effect of U-235 in U(IV)-U(VI) chemical-exchange[END_REF] which yielded a 235 U isotopic enrichment over 238 U in the U(VI) species of an initial 235 U/ 238 U mixture. As a result, Bigeleisen could unambiguously attribute the experimentally observed excess separation factor (beyond the expected regular mass-dependent fractionation) to a nuclear field shift effect rather than a nuclear spin effect [START_REF] Fujii | An anomalous isotope effect of U-235 in U(IV)-U(VI) chemical-exchange[END_REF]. The theory proved also successful to the description of related findings in 233 U/ 238 U isotope studies [START_REF] Nomura | Mass Dependence of Uranium Isotope Effects in the U(IV)-U(VI) Exchange Reaction[END_REF] and ab initio atomic and molecular studies of uranium isotope fractionation have been reported by Abe and co-workers [START_REF] Abe | An ab initio study based on a finite nucleus model for isotope fractionation in the U(III)-U(IV) exchange reaction system[END_REF][START_REF] Abe | An ab initio molecular orbital study of the nuclear volume effects in uranium isotope fractionations[END_REF]. The effect of nuclear volume has also been considered as a possible source of explanation for isotopic anomalies observed in the isotope fractionation of other heavy atoms such as mercury [START_REF] Schauble | Mass-dependent and -independent fractionation of Hg isotopes by photoreduction in aquatic systems[END_REF][START_REF] Bergquist | Mass-dependent and -independent fractionation of Hg isotopes by photoreduction in aquatic systems[END_REF][START_REF] Estrade | Mercury isotope fractionation during liquidvapor evaporation experiments[END_REF][START_REF] Wiederhold | Equilibrium mercury isotope fractionation between dissolved Hg(II) species and thiol-bound Hg[END_REF], as well as observed in the early Solar System [START_REF] Fujii | Nuclear field vs. nucleosynthetic effects as cause of isotopic anomalies in the early Solar System[END_REF].

In the present work we consider the effect of nuclear volume no rotational spectra. In the early 1980s Tiemann and co-workers [START_REF] Tiemann | Observed adiabatic corrections to the Born-Oppenheimer approxmation for diatomic molecules with ten valence electrons[END_REF] conducted a series of high-resolution rotation spectroscopy experiments for isovalent, closed-shell diatomics AB, which aimed at an accurate assessment of spectroscopic constants as for example, the equilibrium bond length R e or the rotational constant B e . In the particular case of thallium and lead compounds, they were not able to bring their measured isotope data in agreement using the existing theory of adiabatic and non-adiabatic correction terms to the Born-Oppenheimer ("clamped nuclei") approximation [START_REF] Watson | The isotope dependence of the equilibrium rotational constants in 1 Σ states of diatomic molecules[END_REF][START_REF] Bunker | The nuclear mass dependence of the Dunham coefficients and the breakdown of the Born-Oppenheimer approximation[END_REF][START_REF] Watson | The isotope dependence of diatomic Dunham coefficients[END_REF], considering their order of magnitude required to derive consistent Dunham coefficients Y kl [START_REF] Dunham | The energy levels of a rotating vibrator[END_REF]. In a follow-up publication Knöckel and Tiemann [START_REF] Knöckel | Isotopic field shift in the transition A 0 + -X 1 Σ + of PbS[END_REF] therefore first unambiguously identified the necessity of, beyond the known mass-dependent corrections, introducing an additional correction factor arising from the finite size of the nuclei. Taking into consideration this finite extension in the initial molecular Hamiltonian, Schlembach and Tiemann [START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF] built a sound theoretical foundation of the molecular field shift in rotational spectra. Their final expression for the total correction term to the first Dunham coefficient Y 01 , which is of particular interest in rotational spectroscopy, thus not only contains mass-dependent contributions but also a further term that is proportional to the change in the mean-square nuclear charge radius δ r 2 n times the derivative dρ(0) dR Re of the electron density ρ(o) at a given nucleus with respect to the internuclear distance R, taken at R e . The so-called contact density ρ(0) is a central quantity in the theory of the Mössbauer isomer shift [START_REF] Mössbauer | Kernresonanzabsorption von Gammastrahlung in 191[END_REF][START_REF] Mössbauer | Kernresonanzfluoreszenz von Gammastrahlung in 191[END_REF][START_REF] Greenwood | Mössbauer Spectroscopy[END_REF][START_REF] Gibb | Principles of Mössbauer Spectroscopy[END_REF][START_REF] Gütlich | Mössbauer Spectroscopy and Transition Metal Chemistry[END_REF][START_REF] Gütlich | [END_REF], but appears as a consequence of approximations [START_REF] Fricke | Calculation of isomer shift in Mössbauer spectroscopy[END_REF][START_REF] Knecht | Mössbauer spectroscopy for heavy elements: a relativistic benchmark study of mercury[END_REF]. The same holds true in the case of the molecular field shift, as will be discussed in section II.

A careful review of the original spectroscopic data for the heavy-element Pb-and Tl-compounds [START_REF] Tiemann | Observed adiabatic corrections to the Born-Oppenheimer approxmation for diatomic molecules with ten valence electrons[END_REF] by means of the expressions derived in Ref. [START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF] revealed then the molecular field shift term as main contributor to the isotope shift corrections to Y 01 . Interestingly, though, in 1985 [START_REF] Knöckel | Molecular-beam-laser studies of the states X 1 Σ + and A0 + of PbS[END_REF] the same experimental group announced a revision of their entire set of earlier results (for a review see for example Ref. 37) based on a computational error where they recommend a scaling of the data by a factor 10. Their revised findings have later been corroborated by Cooke and co-workers [START_REF] Cooke | The calculation of field shift effects in the rotational spectra of heavy metal-containing diatomic molecules using density functional theory: comparison with experiment for the Tl-halides and Pb-chalcogenides[END_REF][START_REF] Cooke | The rotational spectrum, nuclear field shift effects, 195 Pt nuclear spin-rotation constant, and electric dipole moment of PtSi[END_REF][START_REF] Etchison | Born-Oppenheimer breakdown effects and hyperfine structure in the rotational spectrum of strontium monosulfide[END_REF][START_REF] Serafin | Concerning the electron density at the Pb nucleus in PbO as a function of bond length[END_REF] who performed a series of (scalar-relativistic) ab initio calculations based on density-functional theory (DFT).

They furthermore demonstrated the existence of a molecular field shift effect in the rotational spectrum of platinum silicide [START_REF] Cooke | The rotational spectrum, nuclear field shift effects, 195 Pt nuclear spin-rotation constant, and electric dipole moment of PtSi[END_REF] by carrying out high-resolution experiments accompanied by DFT calculations.

The purpose of the present study is to provide an independent check on both theoretical and experimental studies of the isotopic field shift in rotational spectra. On the theoretical side it is known that 1-and 2-component relativistic calculations of molecular properties probing the electron density near nuclei are highly sensitive to picture change errors [START_REF] Baerends | Relativistic atomic orbital contractions and expansions -magnitudes and explanations[END_REF][START_REF] Kellö | Picture change and calculations of expectation values in approximative relativistic theories[END_REF][START_REF] Dyall | Relativistic electric and magnetic property operators for two-component transformed Hamiltonians[END_REF][START_REF] Pernpointer | Accurate nuclear quadrupole moments of the gallium isotopes 69 Ga and 71 Ga within the PCNQM model[END_REF] and so we wanted to calibrate previous calculations against 4-component relativistic highly correlated calculations. Theoretical studies can also shed light on the major physical contributions to the molecular field shift.

The paper is therefore organized as follows: In section II we carefully re-derive an expression for the first-order modification of the rotational constant B e due to nuclear volume changes between isotopes. We have next carried out non-relativistic, scalar-relativistic two-component as well as relativistic two-and four-component ab initio DFT and high-level wave function Coupled Cluster calculations for three representative molecules, namely TlI, PbTe and PtSi.

Computational details are provided in section III, and in section IV we compare our results to previous theoretical and experimental results based on the original as well as revised expressions derived by Tiemann and co-workers. By means of our extensive reference data we shall reveal a curious twist in the tale before concluding in section V.

II. THEORY

Within the Born-Oppenheimer ("clamped nuclei") approximation the rovibrational energy levels of a closed-shell diatomic molecule AB may be determined from an effective radial Schrödinger equation of the form

- 2 2µ d dR 2 + E el (R) + 2 J(J + 1) 2µR 2 ψ ν,J (R) = E ν,J ψ ν,J (R) (1)
where R is the internuclear coordinate and µ the reduced mass. Solutions to this problem of rotating vibrator were provided by Dunham in 1932 [START_REF] Dunham | The energy levels of a rotating vibrator[END_REF] in the framework of the Jeffreys-Wentzel-Kramers-Brillouin (JWKB) approximation [START_REF] Dunham | The energy levels of a rotating vibrator[END_REF][START_REF] Jeffreys | On certain approximate solutions of linear differential equations of the second order[END_REF][START_REF] Wentzel | Eine Verallgemeinerung der Quantenbedingungen für die Zwecke der Wellenmechanik[END_REF][START_REF] Kramers | Wellenmechanik und halbzählige Quantisierung[END_REF][START_REF] Brillouin | La mécanique ondulatoire de Schrödinger: une méthode générale de resolution par approximations successives[END_REF] and expressed in the form nowadays known as the Dunham expansion

E ν,J = h k=0 l=0 Y kl (ν + 1/2) k [J (J + 1)] l . ( 2
)
where ν and J are vibrational and rotational quantum numbers, respectively. E el is the potential obtained from solving the electronic problem

H el ψ el r el ; R = E el (R)ψ el r el ; R (3)
where r el designates all electronic coordinates. The electronic Hamiltonian, relativistic or not, has the generic form

H el = i h(i) + 1 2 i =j g(i, j) + V AB ; h = h 0 + V eN (4)
where V AB is the classical repulsion of nuclei A and B and g(i, j) the two-electron operator. The one-electron operator h splits into the free-particle Hamiltonian h 0 and a term V eN describing the interaction with the nuclei.

Following Schlembach and Tiemann [START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF] we give the Dunham coefficients Y kl in units of frequency. As emphasized by Ogilvie [START_REF] Ogilvie | Vibrational and Rotational Spectrometry of Diatomic Molecules[END_REF][START_REF] Ogilvie | Dunham's formalism applied in reduction of spectral data of diatomic molecules and the development of computational spectrometry[END_REF] they are not freely adjustable fitting coefficients since they are interrelated. To lowest order in ν and J one has [START_REF] Dunham | The energy levels of a rotating vibrator[END_REF] [START_REF] Ishikawa | On the use of an extended nucleus in Dirac-Fock Gaussian basis set calculations[END_REF] where R e is the equilibrium bond distance. For a specific isotopomer α of the molecule AB the Dunham coefficients are given by

Y 10 = ωe 2π 1 + B 2 e /4ν
Y α kl = µ -(k/2+l) α U kl (6)
where U kl are isotope-independent coefficients. However, the above relation supposes i) that all isotopomers experience the same internuclear potential and ii) that the first-order JWKB approximation (the Bohr-Sommerfeld quantization condition) is exact [START_REF] Le Roy | Improved parametrization for combined isotopomer analysis of diatomic spectra and its application to HF and DF[END_REF]. A more general expression, taking into account the breakdown of the Born-Oppenheimer approximation, was proposed by Ross et al. [START_REF] Ross | Heterodyne measurements of 12 C 18 O, 13 C 16 O, and 13 C 18 O laser frequencies; mass dependence of Dunham coefficients[END_REF] in an experimental study of isotopomers of CO

Y α kl = µ -(k/2+l) α U kl 1 + m e M A ∆ A kl + m e M B ∆ B kl (7)
and later developed theoretically by Bunker [START_REF] Bunker | The nuclear mass dependence of the Dunham coefficients and the breakdown of the Born-Oppenheimer approximation[END_REF] and Watson [START_REF] Watson | The isotope dependence of diatomic Dunham coefficients[END_REF]. In the above expression ∆ A kl and ∆ B kl are massindependent atom-specific correction factors.

In the following we shall focus on the modification of the rotational constant by modification of the size of one nucleus, say nucleus A. Experimentally only discrete values of the nuclear radius are available through isotopic substitution A → A , but for the purpose of derivation it is more useful to consider the electronic energy E el as a continuous function of both internuclear distance R and the nucleus radius ξ, that is E el ≡ E el (R, ξ). We shall let ξ A correspond to the nuclear radius of a particular reference isotope of atom A. The equilibrium internuclear distance R e for any value of nuclear radius ξ is found by minimizing the electronic energy with respect to internuclear distance,

dE el dR ξ = 0, (8)
and thus becomes a function of the nuclear radius parameter, which allow us to attack the above problem by variational perturbation theory (see for instance Ref. [START_REF] Helgaker | Calculation of geometrical derivatives in molecular electronic structure theory[END_REF]). Since the variational condition Eq. ( 8) is valid for any nuclear radius ξ we find the first-order shift in the equilibrium distance by derivation with respect to nuclear radius ξ

d 2 E el dξdR ξ A = d 2 E el dR 2 dR dξ + ∂ ∂ξ dE el dR ξ A = 0 (9)
The expression for the first-order equilibrium distance, Eq. ( 9), contains the explicit derivative of the electronic energy E el with respect to nuclear radius. Such derivatives also appear in recent theoretical studies of the isomer shift in Mössbauer spectroscopy [START_REF] Knecht | Mössbauer spectroscopy for heavy elements: a relativistic benchmark study of mercury[END_REF][START_REF] Filatov | On the calculation of Mössbauer isomer shift[END_REF][START_REF] Filatov | First principles calculation of Mössbauer isomer shift[END_REF]. We can safely ignore contributions from the classical repulsion V AB of nuclei such that the only non-zero contribution comes from the electrostatic interaction between the electrons and nucleus A. Its modification upon the isotope substitution A → A is given by

δE el A A = ρ e (r e ; R) φ A (r e ; ξ A ) -φ A (r e ; ξ A ) dτ e . (10)
It should be noted that in the above expression we ignore the implicit dependence on nuclear radius of the electronic charge distribution ρ e . We express the scalar potential of nucleus A in terms of the normalized nuclear charge

distribution ρ A φ A (r e ; ξ) = Ze 4πε 0 ρ A (r n ; ξ) r ne dτ n . (11)
The perhaps simplest model for a finite nuclear charge distribution is the homogeneous charged sphere of radius ξ

ρ H n (r) =      ρ 0 ; r ≤ ξ 0 ; r > ξ ; ρ -1 0 = 4π 3 ξ 3 (12)
which is also the model considered by Schlembach and Tiemann [START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF]. The associated potential is

φ H n (r) =      Ze 8πε0ξ (3 -r 2 ξ 2 ) ; r ≤ ξ Ze 4πε0r ; r > ξ (13)
More widely employed in relativistic molecular calculations is the Gaussian model [START_REF] Ishikawa | On the use of an extended nucleus in Dirac-Fock Gaussian basis set calculations[END_REF][START_REF] Visser | The use of Gaussian nuclear charge distributions for the calculation of relativistic electronic wavefunctions using basis set expansions[END_REF] in which the charge distribution is given by

ρ G n (r) = ρ 0 exp -r 2 /ξ 2 G ; ρ -3 0 = π 1/2 ξ G (14)
with the associated potential

φ G n (r) = Ze 4πε 0 r erf(r/ξ G ). ( 15
)
The radius parameter ξ G can be connected to the radius ξ of the homogeneous sphere by requiring identical second radial momenta

< r 2 n >= r 2 n ρ n (r n ) dτ n = 3 2 ξ 2 G = 3 5 ξ 2 ; ⇒ ξ G = 2 5 ξ (16)
Schlembach and Tiemann [START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF] consider the modification of the internuclear potential when going from a point nucleus to a finite nucleus for each isotope, that is

δE el A0 = ρ e (r e ; R) φ A (r e ; ξ A ) -φ A (r e ; 0) dτ e . ( 17
)
in order to introduce mass-independent Dunham coefficients U kl for a fictitious molecule of point-like nuclei. However, such an approach is problematic in the relativistic case since the electronic density display a weak singularity at point nuclei. As such, it is perhaps better to modify the Dunham expansion by the introduction of a reference isotopomer, as suggested by Le Roy [START_REF] Le Roy | Improved parametrization for combined isotopomer analysis of diatomic spectra and its application to HF and DF[END_REF]. However, we can formally write Eq. ( 10) as

δE el A A = δE el A 0 -δE el A0 . ( 18
)
where the electron density will be that of a molecule with finite size nuclei. In passing we note that the difference potential in Eq. ( 17), using the homogeneous charged sphere model, looks curiously similar to the potential of the Rutherford atom, consisting of a positive point charge at the origin and the homogeneous electronic charge within a sphere of radius ξ [START_REF] Rutherford | The Scattering of α and β particles by matter and the structure of the atom[END_REF]. We can formally write Eq. ( 17) as

δE el A0 = ρe φ A (r e ; ξ A ) -φ A (r e ; 0) dτ e . ( 19
)
where we have introduced a constant effective electronic charge density ρe . Due to the extreme short-range nature of the above difference potential, as seen for instance from Eq.( 13), the effective density ρe is typically approximated by the contact density, that is the value of the electronic charge density at the nucleus:

ρe ≈ ρ e (0) . (20)
For light atoms this is certainly an excellent approximation; for the more extended nuclei of heavier atoms this leads to an overestimation [START_REF] Fricke | Calculation of isomer shift in Mössbauer spectroscopy[END_REF]. In the case of the Mössbauer isomer shift the deviation is quite systematic in nature, which suggests that it can be handled by a correction factor [START_REF] Knecht | Mössbauer spectroscopy for heavy elements: a relativistic benchmark study of mercury[END_REF], a feature that will also be investigated in the present study.

Using Eq.( 18) we obtain the expression

δE el A A = ρe Ze 6ε 0 δ < r 2 n > A A ; δ < r 2 n > A A =< r 2 n > A -< r 2 n > A (21)
which holds for both the homogeneous charged sphere and the Gaussian model of the nuclear charge distribution.

Following Filatov [START_REF] Filatov | On the calculation of Mössbauer isomer shift[END_REF] an alternative approach to δE el A A is to approximate the potential difference of Eq.( 10) by a first-order Taylor expansion, that is

δE el A A ≈ ρ e (r e ; R) ∂φ A (r e ; ξ) ∂ξ ξ A ξ A -ξ A dτ e = ∂ ∂ξ dE el dR ξ A ξ A -ξ A (22)
thus connecting to the expression for the first-order equilibrium distance, Eq. ( 9).

Using this connection we find that the change in equilibrium bond distance due to change in nuclear size upon isotope substitution A → A can be expressed as

δR e = - Ze 2 6ε 0 k A e δ < r 2 n > A A dρ e dR R A e ( 23
)
where R A e and k A e corresponds to the equilibrium bond distance and the force constant of the reference isotopomer.

The corresponding modification of the equilibrium rotational constant B e is

δB e = B A e -2 R e δR e = B A e V A δ < r 2 n > A A (24)
where appears the factor

V A = Ze 2 3ε 0 k A e R A e dρ e dR R A e (25)
It has exactly the same form as Eq. ( 25) of the 1982 paper by Schlembach and Tiemann [START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF] except that the spectroscopic constants R A e and k A e is that of the reference isotopomer and not of a fictitious molecule with point-like nuclei.

Curiously, in a paper [START_REF] Knöckel | Molecular-beam-laser studies of the states X 1 Σ + and A0 + of PbS[END_REF] from 1985 by the same group the formula for the isotopic field shift factors is given with an additional factor π 2

V A → π 2 V A (26)
which to our opinion is not justified. However, as we shall see in section IV it contributes to a strange twist in our story.

III. COMPUTATIONAL DETAILS

All molecular calculations reported in this paper have been carried out with the DIRAC10 program package [START_REF] Dirac | a relativistic ab initio electronic structure program, Release DIRAC10[END_REF].

A. Electron correlation methods

The absolute magnitude of the electron density in the vicinity of a given nucleus is well described within a meanfield approach to the electron-electron interaction since the dominant part is governed by the influence of the nuclear potential and deformations of core orbitals are expensive in energy. The significant relative changes in the contact density, which yield the most sensitive contributions to the isotopic field shift, are, on the other hand, affected by subtle alterations in the valence electronic structure in a varying chemical environment. An accurate description of these valence contributions thus requires proper inclusion of electron-electron correlation, as has been shown recently in the context of isomer shift predictions for Mössbauer spectroscopy [START_REF] Knecht | Mössbauer spectroscopy for heavy elements: a relativistic benchmark study of mercury[END_REF][START_REF] Kurian | Calibration of 57 Fe isomer shift from ab initio calculations: can theory and experiment reach an agreement ?[END_REF]. For the three closed-shell diatomics TlI, PbTe and PtSi we therefore performed single-reference coupled-cluster (CC) calculations with a full iterative treatment of single and double excitations (CCSD) and including perturbative corrections for triple excitations (CCSD(T)) [START_REF] Raghavachari | A fifth-order perturbation comparison of electron correlation theories[END_REF][START_REF] Visscher | Formulation and implementation of a relativistic unrestricted coupled-cluster method including noniterative connected triples[END_REF][START_REF] Pernpointner | Parallelization of four-component calculations. II. Symmetry-driven parallelization of the 4-Spinor CCSD algorithm[END_REF].

The central object of our study is the first geometrical derivative of the contact density for selected nuclei

X calculated at the equilibrium internuclear distance R e ρ [1]X 0 = dρ X (0) dR Re (27)
which we in the following shall refer to as the contact density gradient. Presently there is no analytic implementation of CC expectation values in the DIRAC10 program package. For the calculation of contact densities at the CC level we therefore pursued a finite-field (ff) approach using a computational protocol which is described at length elsewhere [START_REF] Knecht | Mössbauer spectroscopy for heavy elements: a relativistic benchmark study of mercury[END_REF] and shall only be briefly sketched in the following. Exploiting the additivity of the contact density contributions ρ(0) = ρ HF (0)+ρ corr (0), namely, (i) the Hartree-Fock expectation value of ρ HF (0), which can be evaluated analytically in DIRAC10, and (ii) the electron-electron correlation term ρ corr (0), we are only left with the determination of the latter contribution in a finite-field scheme. In line with our previous ff-coupled cluster property calculations [START_REF] Knecht | Mössbauer spectroscopy for heavy elements: a relativistic benchmark study of mercury[END_REF] we employed an optimal finite-field parameter of 10 -8 for all diatomic systems under consideration and took take advantage of the central-difference method [START_REF] Hildebrand | Introduction to Numerical Analysis[END_REF] using a seven-point stencil for the numerical differentiation. In order to obtain the contact density gradient ρ

[1]X 0 we performed a second numerical differentiation by means of the central-difference method using a step size of 0.0125 Å as "finite-field" parameter. Computing the contact density gradient as derivative of a polynomial data fit function yielded equally identical results.

Besides the isotopic field shift evaluations based on the wave function Hartree-Fock and CC methods we also carried out four-and two-component density functional theory (DFT) calculations. Aiming at an assessment of both the accuracy and internal consistency within the DFT contact densities (and therefore the isotopic field shift) we employed an ample set of exchange-correlation functionals, namely LDA (VWN5) [START_REF] Dirac | Note on exchange phenomena in the Thomas Atom[END_REF][START_REF] Vosko | Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis[END_REF], BLYP [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF][START_REF] Lee | Development of the Colle-Salvetti correlation-energy formula into a functional of the electron-density[END_REF][START_REF] Miehlich | Results obtained with the correlation-energy density functionals of Becke and Lee, Yang and Parr[END_REF], B3LYP [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF][START_REF] Stephens | Ab-initio calculation of vibrational absorption and circular-dichroism spectra using density-functional force-fields[END_REF][START_REF] Hertwig | On the parametrization of the local correlation functional: What is Becke-3-LYP?[END_REF], CAMB3LYP [START_REF] Yanai | A new hybrid exchange-correlation functional using the Coulomb-attenuating method (CAM-B3LYP)[END_REF], PBE [START_REF] Perdew | Generalized gradient approximation made simple[END_REF], PBE0 [START_REF] Perdew | Rationale for mixing exact exchange with density functional approximations[END_REF], and furthermore the SAOP model potential [START_REF] Schipper | Molecular calculations of excitation energies and (hyper)polarizabilities with a statistical average of orbital model exchange-correlation potentials[END_REF]. Of main concern for electron density evaluations in the core region is the use of a sufficiently dense integration grid in the numerical integration of the exchange-correlation evaluation [START_REF] Neese | Prediction and interpretation of the 57 Fe isomer shift in Mössbauer spectra by density functional theory[END_REF]. We met this particular requirement by employing throughout all DFT calculations an ultrafine grid.

B. Hamiltonian

An important aspect of the present study is furthermore to address successive approximations to the four-component DC Hamiltonian and their validity in the evaluation of the field shift effect for a given nucleus. In particular, we therefore compare four-component relativistic results to:

• relativistic and scalar-relativistic data using the eXact 2-Component (X2C) Hamiltonian [START_REF] Iliaš | An infinite-order two-component relativistic hamiltonian by a simple one-step transformation[END_REF],

• four-component spin-orbit free (sf) [START_REF] Dyall | An exact separation of the spin-free and spin-dependent terms of the Dirac-Coulomb-Breit Hamiltonian[END_REF][START_REF] Kutzelnigg | Basis set expansion of the Dirac operator without variational collaps[END_REF] (scalar-relativistic) results,

• non-relativistic (NR) values employing the Lévy-Leblond [START_REF] Lévy-Leblond | Nonrelativistic particles and wave equations[END_REF] Hamiltonian.

The X2C calculations have been carried out either including 2-electron spin-same-orbit corrections provided by the AMFI [START_REF] Amfi | An Atomic Mean-Field Code. B. Schimmelpfennig[END_REF][START_REF] Hess | A mean-field spin-orbit method applicable to correlated wavefunctions[END_REF] code (relativistic) or by retaining only the spin-free terms in the one-electron Hamiltonian prior to the transformation to two-component basis (scalar-relativistic; in the following denoted X2C-sf).

For reasons of computational efficiency, a molecular mean-field approximation 4c DC * * to the four-component DC Hamiltonian 4c DC was applied in the majority of CCSD(T) calculations, where our notation strictly follows the Hamiltonian hierarchy introduced in Ref. [START_REF] Sikkema | The molecular mean-field approach for correlated calculations[END_REF]. The relative deviation in the contact density gradient compared to the exact 4c DC-CCSD(T) value was in all cases tested less than 0.05%.

C. Basis sets

All molecular calculations were performed using atom-centered large-component basis sets consisting of uncontracted scalar Gaussian type orbitals (GTO) where the small-component basis functions were generated, where appropriate, by the restricted kinetic balance condition as implemented in DIRAC10 [START_REF] Visscher | Approximate relativistic electronic structure methods based on the quaternion modified Dirac equation[END_REF]. In the case of the heavy elements Pb, Te, Tl, I and Pt we used the triple-ζ (TZ) and quadruple-ζ (QZ) basis sets of Dyall [START_REF] Dyall | Relativistic quadruple-zeta and revised triple-zeta and double-zeta basis sets for the 4p, 5p, and 6p elements[END_REF][START_REF] Dyall | Relativistic double-zeta, triple-zeta, and quadruple-zeta basis sets for the 5d elements Hf-Hg[END_REF][START_REF] Dyall | Revised relativistic basis sets for the 5d elements Hf-Hg[END_REF]. The basic large-component SCF set of primitives was further augmented with the recommended correlating and polarizing functions in order to properly account for correlation contributions from the nsp valence and outer-core (n-1)d shells. For the Pt TZ basis we also added a set of 2g1h primitives to allow for correlation of the (n-2)f electrons.

Basis set saturation at the heavy nuclei of interest with respect to an accurate computation of the contact density was achieved by further augmenting the set of primitives in an even-tempered fashion with two more tight s and one tight p functions, which is in line with previous studies [START_REF] Knecht | Mössbauer spectroscopy for heavy elements: a relativistic benchmark study of mercury[END_REF][START_REF] Mastalerz | Basis set representation of the electron density at an atomic nucleus[END_REF]

D. Active space considerations

We have shown in a recent benchmark work for mercury-containing compounds [START_REF] Knecht | Mössbauer spectroscopy for heavy elements: a relativistic benchmark study of mercury[END_REF] that an accurate assessment of correlation effects to the contact density ρ(0) by means of wave-function based methods necessitates the inclusion of both core-valence and valence correlation contributions. Our active space in all finite-field coupled-cluster calculations therefore comprised the 5d6s6p shells of Tl (Pb) and 4d5s5p shells of I (Te) of the thallium iodide (lead telluride) compound. For the platinum silicide diatom the active space has been adapted accordingly with an explicit correlation treatment of the Pt 4f 5p5d6s6p and Si 2s2p shells.

The size of the virtual spinor space for all three molecular systems was tailored to contain all recommended core-and valence correlation as well as valence dipole polarization functions. This corresponds for the TlI and PbTe compounds to a threshold of 40 hartree whereas for the PtSi diatom the cutoff is fixed at 62 hartree. The validity of this choice was confirmed by calculating the contact density gradient ρ

[1]T l 0

in thallium iodide at the CCSD(T) level using an enlarged virtual space threshold of 134 hartree. Gradually saturating the unoccupied space by this means led to a (for the present purpose) negligible decrease of 0.1% in the final value of ρ for this molecule.

We conclude this paragraph by noting that all geometrical derivatives, subject to discussion as follows, have been taken at the experimental equilibrium internuclear distances R exp e summarized in Table I. This practice ensures a fair comparison with the computational results by Cooke et al. [START_REF] Cooke | The calculation of field shift effects in the rotational spectra of heavy metal-containing diatomic molecules using density functional theory: comparison with experiment for the Tl-halides and Pb-chalcogenides[END_REF][START_REF] Cooke | The rotational spectrum, nuclear field shift effects, 195 Pt nuclear spin-rotation constant, and electric dipole moment of PtSi[END_REF] that is to be discussed in Section IV E.

B. On the use of the contact density

Before embarking on the analysis and interpretation of our geometrical derivative data of the electronic charge density in the light of experimental and earlier theoretical predictions, we first investigate the validity of our approximation [START_REF] Fujii | Nuclear field vs. nucleosynthetic effects as cause of isotopic anomalies in the early Solar System[END_REF] of the effective density ρe by the contact density ρ e (0). This approach leads to a modified expression of Eq. 25, which describes the factor V A connected to the molecular field shift:

V A ≈ Ze 2 3ε 0 k A e R A e ρ [1]A 0 (28)
This is the formula from which we have derived our V A values that enter the discussion in Section IV E. and p 1/2 shells (as expected from theory) whereas the effective density yields additional contributions considerable in their magnitude first and foremost from p 3/2 orbitals. The latter findings can be understood as a result of sizable p 3/2 orbital values within the nuclear volume. (ii) though contributions from the s 1/2 and p 1/2 orbitals are throughout significantly higher in the contact density than in the effective density approach, they appear to be on a rather systematic basis of ≈ +10%, a fact, which has also been observed for the mercury atom [START_REF] Knecht | Mössbauer spectroscopy for heavy elements: a relativistic benchmark study of mercury[END_REF]. We expect the contact density to be an increasingly better approximation of the effective density with decreasing nuclear charge; for xenon the deviation is found to be around 5% [92].

Significant for the present discussion of the molecular field shift is, however, the first geometrical derivative rather than the absolute magnitude of the electronic charge density at a given nucleus and internuclear distance. Table III summarizes for our three reference molecules the effective as well as contact density gradients ρ [1]X computed at the four-component HF/TZ+2s1p and DFT/TZ+2s1p level, respectively. Yet again we find considerable deviations for the first derivative between both approaches. The use of the contact density approximation for the study of nuclear size effects in the rotational spectra of Pt, Tl, and Pb diatomics, may, nevertheless, be well justified since the deviations exhibit a systematic nature on the order of 10% and, equally important, are independent of the level of electron-electron correlation included.

13

C. Projection analysis of the contact density

At the SCF level the electronic density at nuclear center X may be written as

ρ X = -e nocc i ψ i |δ (r -R X )| ψ i (29)
Inserting a Linear Combination of Atomic Orbitals (LCAO) expansion of the molecular orbitals into the above expression leads to a projection analysis [START_REF] Knecht | Mössbauer spectroscopy for heavy elements: a relativistic benchmark study of mercury[END_REF][START_REF] Bast | Analysis of parity violation in chiral molecules[END_REF] of expectation values particularly useful for local properties such as the contact density. In practice the LCAO-expansion is limited to atomic orbitals p which are occupied in the ground states of atoms P

|ψ i = pP ψ P p c P pi + ψ pol i ( 30
)
and to which is added the orthogonal complement ψ pol i which is denoted the polarization contribution. Expectation values are accordingly decomposed into inter-and intra-atomic as well as polarization contributions. As shown in Ref. [START_REF] Knecht | Mössbauer spectroscopy for heavy elements: a relativistic benchmark study of mercury[END_REF] the contact density ρ X 0 is completely dominated by the intra-atomic contribution of same center X and can accordingly be expressed as

ρ X 0 ≈ -e pq R L p R L q + R S p R S q r=R X D XX qp ; D QP qp = i c Q qi c P * pi (31)
where R L and R S are the large and small component radial functions, respectively. Non-zero contributions to the contact density are provided exclusively by the large components of s 1/2 orbitals and the small components of p 1/2 orbitals, whereas for the effective density ρe other orbitals may come into play, as seen in section IV B. The radial functions R X are evidently independent of internuclear distance, so all geometry-dependence arises from the density matrix D XX qp expressed in terms of atomic orbital expansion coefficients.

From Table IV it can be seen that the intra-atomic contribution from the Thallium atom dominates the contact number density gradient ρ

[1]T l 0 in TlI. Assuming generality of this result, we can accordingly express the contact density gradient as

ρ [1]X 0 ≈ -e pq R L p R L q + R S p R S q r=R X dD XX qp dR Re . ( 32
)
It is furthermore seen that a negative contribution from diagonal elements (p = q), which would contribute to the atomic expectation value, are overwhelmed by positive off-diagonal contributions, which come into play due to the breakdown of atomic symmetry in the molecule. The two opposing contributions reflect re-organization of the electron density as a function of internuclear distance. A more detailed breakdown of the expectation value is seen in Table V. Looking first at just the geometric D [1]T l derivative of the atomic density matrix, it can be seen that the abovementioned density re-organization takes predominantly place, as expected, amongst the valence orbitals. However, due to the very much larger values of the radial functions of core orbitals at the nucleus, the contribution of core orbitals to the overall expectation value is not negligible. Finally, it may be noted that there are non-zero contributions from the p 1/2 orbitals, but significantly smaller than the contributions from s 1/2 orbitals.

D. Effect of Hamiltonian, basis and method

We summarize the results of our calculations on TlI, PbTe and PtSi in Tables VI , VII and VIII, respectively.

Before comparing with available experimental data we consider the effect of Hamiltonian, method and basis sets on the calculated geometrical density derivatives. Starting from the relativistic DFT(SAOP) value ρ

Any 2-component one-electron property operator should be obtained by the same procedure, that is 2c Ω = U † 4c ΩU ++ , rather than simply taking the large-large (LL) block of the parent 4-component operator 2c Ω ≈ [ 4c Ω] LL , an approximation that may lead to significant errors, in particular for properties probing electron density near nuclei (cf. Ref. [START_REF] Pernpointer | Accurate nuclear quadrupole moments of the gallium isotopes 69 Ga and 71 Ga within the PCNQM model[END_REF], as is the case here. The correct expression of the electron charge density in some point P at the 2-component relativistic SCF level is

ρ 2c = -e Nocc i ψ 2c i | U † δ (r -P) U ++ |ψ 2c i = -e Nocc i ψ 2c † i (P)ψ 2c i (P) (34)
Use of the untransformed operator, that is, the ψ † ψ expression known from the 4c DC and NR levels, leads to an overestimation of ρ in PbTe and PtSi, respectively, showing that for this property picture change errors are significant larger than using the correctly transformed operator, [START_REF] Fricke | Calculation of isomer shift in Mössbauer spectroscopy[END_REF], leads to errors less than 1 %, which is quite acceptable. In those calculations the major source of deviation from the 4 DC is probably the incomplete transformation of the 2-electron operator. For ρ

[1]P t 0

in PtSi results obtained with the Douglas-Kroll-Hess Hamiltonian are available [94]. At the DKH(2,0) level, that is, using a 2nd order DKH Hamiltonian for the generation of orbitals and a 0th order (untransformed) property operator, we get picture change errors on the same order as above, whereas at the DKH [START_REF] Heilig | Changes in mean-square nuclear charge radii from optical isotope shifts[END_REF][START_REF] Heilig | Changes in mean-square nuclear charge radii from optical isotope shifts[END_REF] level the deviation with respect to the reference B3LYP/ 4 DC value is within 3 %.

The effect of the basis set can be seen from the HF/ 4 DC results obtained for ρ A first indication of the importance of electron correlation effects is obtained by comparing SAOP and HF results obtained at the 4c DC level. For all three molecules we observe that introduction of electron correlation through the asymptotically corrected SAOP functional reduces the HF value by 30-40 %, clearly indicating that the inclusion of electron correlation is mandatory for this property. However, taking the CCSD(T)/ 4c DC * * value as reference, we see that the performance of the various DFT functionals is variable and hardly satisfying. For TlI we see that LDA and GGA functionals underestimate the value of ρ

[1]T l 0 by some 20 %, whereas the inclusion of exact exchange through hybrid functionals improves the results somewhat. Interestingly, the long-range corrected functional CAMB3LYP

shows quite different performance from that of B3LYP. For ρ Finally, turning to the CCSD(T) results themselves, we see that for ρ

[1]T l 0 in TlI the molecular-mean field scheme 4c DC * * introduces negligible error compared to the full 4c DC calculation. In the former approach [START_REF] Sikkema | The molecular mean-field approach for correlated calculations[END_REF] all two-electron integral classes are employed in the SCF optimization step, but at the CC level only the (LL|LL) integrals, involving the large components only, are retained and employed in conjunction with the Fock matrix (orbital energies) obtained in the SCF step. Comparing TZ+2s1p and QZ+2s1p values we furthermore see that the value of the density derivative

ρ [1]T l 0
is well converged at the TZ+2s1p level. For TlI we also investigated the convergence with respect to the energy cutoff value for the virtual space and found that the selected value introduces very small errors.

E. Comparison with experimentally derived data

We list in Table IX for the three reference molecules our calculated 4c DC -CCSD(T)/TZ+2s1p contact density gradients in comparison with existing experimental and theoretical data. We furthermore provide estimates for the field shift parameter V X by inserting in Eq. 28 the calculated contact density gradient along with the experimental spectroscopic constants reported in Table I. In passing we note that if we combine our calculated data with the changes in mean-square nuclear radii δ < r 2 n > from optical isotope shifts (only available with unknown screening factor β of order unity) reported by Heilig and Steudel [START_REF] Heilig | Changes in mean-square nuclear charge radii from optical isotope shifts[END_REF], we obtain changes in equilibrium bond lengths for isotope pairs δR Our point of departure for the following discussion will be the experimental Tl (Pb) density gradients, ρ [1]T l =120.6 Å-4 (ρ [1]P b =148 Å-4 ), and field shift parameters, V T l =3.20×10 4 Å-2 (V P b =2.12×10 4 Å-2 ), which were reported by Schlembach and Tiemann in the early 1980s [START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF] and are compiled in the rows "experiment[A]" in Table IX. We find for both heavy metals a very good agreement with our reference CCSD(T) values where we notice slightly larger deviations for the Tl nucleus. Scaling down the density gradient by 10 %, as suggested by our study in Section IV B on the validity of approximating the effective density by the contact density, would bring our CCSD(T) values in even closer agreement with "experiment[A]". Based on our re-derivation of theory and our computed results we are thus confident to have proven wrong a suggested re-scaling of the "experiment[A]" density gradients by a factor 10 which Knöckel and co-workers remarked in a follow-up publication in 1985 [START_REF] Knöckel | Molecular-beam-laser studies of the states X 1 Σ + and A0 + of PbS[END_REF] (denoted as "experiment[B]" in Table IX). More curiously, the authors obtained in the latter work a V P b field shift factor for 208 Pb 32 S which disagrees by approximately one order of magnitude with the first estimate by Schlembach and Tiemann! Consider next the molecular field shift data derived from DFT(SAOP)/QZ4P calculations by Cooke et al. [START_REF] Cooke | The calculation of field shift effects in the rotational spectra of heavy metal-containing diatomic molecules using density functional theory: comparison with experiment for the Tl-halides and Pb-chalcogenides[END_REF].

Their values for the contact density gradients in TlI and PbTe agree reasonably well with the revised experimental data from 1985 [START_REF] Knöckel | Molecular-beam-laser studies of the states X 1 Σ + and A0 + of PbS[END_REF], and thus seem to corroborate those. However, a comparison of their predicted contact density gradient ρ [1]T l =21.1 Å-4 at the Tl nucleus with our four-component CCSD(T) value of ρ [1]T l =142.26 Å-4 reveals a rather striking discrepancy by more than a factor 6. We have demonstrated in Section IV D that DFT density gradients cover a range relative to the CCSD(T) reference of up to 25% indicating that this cannot account for the wide variance of the DFT(SAOP)/QZ4P value. A probable explanation may, on the other hand, be deduced from our non-relativistic (NR) contact density gradient and associated field shift parameter data, which is for completeness added to Table IX for all three diatomics. The evident agreement (apart from minor basis set effects) of the two DFT(SAOP) data sets for either molecule strongly implies that Cooke and co-workers by accident carried out non-relativistic calculations though they were aiming for scalar-relativistic studies based on the zeroth-order regular approximation (ZORA); this has been confirmed for PbS [95]. This conclusion is further corroborated by recalling from Table VI our scalar-relativistic DFT(SAOP)/TZ+2s1p contact density gradient at the Tl (Pb) nucleus ρ [1]T l =134.48 Å-4 (ρ [1]T l =174.59 Å-4 ) which, inserted in Eq. 28 yields a field shift parameter V T l (V P b ) on the order of 3.6 × 10 4 Å-2 (2.52 × 10 4 Å-2 ) rather than 0.61 × 10 4 Å-2 (0.33 × 10 4 Å-2 ). The close agreement of the nonrelativistic DFT(SAOP)/QZ4P data with the (presumably) erroneous late experimental results is merely fortuitous and this fact may explain why the computational mistake in Refs. [START_REF] Cooke | The calculation of field shift effects in the rotational spectra of heavy metal-containing diatomic molecules using density functional theory: comparison with experiment for the Tl-halides and Pb-chalcogenides[END_REF][START_REF] Cooke | The rotational spectrum, nuclear field shift effects, 195 Pt nuclear spin-rotation constant, and electric dipole moment of PtSi[END_REF][START_REF] Etchison | Born-Oppenheimer breakdown effects and hyperfine structure in the rotational spectrum of strontium monosulfide[END_REF] was left undiscovered to date.

With that said, we finally turn to the PtSi field shift parameters compiled at the lower end of Table IX. Here, we see again a close match of the contact density gradient ρ [1]P t and field shift factor V P t derived from our non-relativistic DFT(SAOP) calculations and the purportedly scalar-relativistic DFT(SAOP) calculations, which further supports our present conclusions. The open question remaining is then to explain the considerably large disagreement between our four-component 4c DC CCSDT(T)/TZ+2s1p field shift factor V P t = -5.05× 10 4 Å-2 and the experimental value of -0.72 × 10 4 Å-2 reported by Cooke and co-workers [START_REF] Cooke | The rotational spectrum, nuclear field shift effects, 195 Pt nuclear spin-rotation constant, and electric dipole moment of PtSi[END_REF]. Based on the findings so far, we are confident of our theoretical estimate and therefore propose a careful re-examination of the experimental platinum silicide data fit, in particular in the light of the non-conforming Eqs. 25 and 26, the latter probably being employed in the present data fit.

In summary, our extensive four-component DFT and high-level CCSD(T) reference data clearly shows by comparison with existing experimental results for the contact density gradient and field shift factor in rotational spectra of the heavy-element diatomics TlI, PbTe, and PtSi that the original work of Schlembach and Tiemann in 1982 [START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF] provided correct estimates for both properties. We furthermore provide evidence that the revised experimental data of Knöckel and co-workers [START_REF] Knöckel | Molecular-beam-laser studies of the states X 1 Σ + and A0 + of PbS[END_REF] as well as more recent experimental data by Cooke et al. [START_REF] Cooke | The rotational spectrum, nuclear field shift effects, 195 Pt nuclear spin-rotation constant, and electric dipole moment of PtSi[END_REF][START_REF] Etchison | Born-Oppenheimer breakdown effects and hyperfine structure in the rotational spectrum of strontium monosulfide[END_REF] are most likely incorrect and should be re-examined. We also conclude that previous DFT calculations [START_REF] Cooke | The calculation of field shift effects in the rotational spectra of heavy metal-containing diatomic molecules using density functional theory: comparison with experiment for the Tl-halides and Pb-chalcogenides[END_REF][START_REF] Cooke | The rotational spectrum, nuclear field shift effects, 195 Pt nuclear spin-rotation constant, and electric dipole moment of PtSi[END_REF][START_REF] Etchison | Born-Oppenheimer breakdown effects and hyperfine structure in the rotational spectrum of strontium monosulfide[END_REF], which appeared to support the more recent experimental values, are seriously in error since they are accidentally based on a non-relativistic Hamiltonian.

V. CONCLUSIONS AND PERSPECTIVES

The objective of this study has been to provide an independent evaluation of both theoretical and experimental studies of the isotopic field shift in rotational spectra using high-level relativistic four-component electronic structure methods. The first experimental evidence for an isotopic field shift effect in rotational spectra of the heavy-element Pb-chalcogenides and Tl-halides were reported by Knöckel and Tiemann [START_REF] Knöckel | Isotopic field shift in the transition A 0 + -X 1 Σ + of PbS[END_REF] and Schlembach and Tiemann [START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF] who identified it as "hidden" terms in the adiabatic mass-dependent correction to the Born-Oppenheimer ("clamped nuclei") approximation (see for example [START_REF] Watson | The isotope dependence of diatomic Dunham coefficients[END_REF]). The new atom-specific correction terms were shown to be proportional to the gradient of the electron (contact) density at a given nucleus with respect to the internuclear distance [START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF], a quantity that can be extracted from electronic structure calculations. For the present benchmark study we therefore chose three molecules TlI, PbTe, and PtSi, for which experimental and other theoretical data is available.

In the original formulation of the molecular isotopic field shift by Schlembach and Tiemann [START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF] the effective electron density at a given nucleus was identified as the contact density. We have investigated the validity of this approximation.

Using a Gaussian model of the nuclear charge distribution we find that the contact density approach yields a systematic overestimation of about 10% in all three molecules compared to the electron density gradient derived from the effective density, a deviation that can clearly not be ignored in high-precision work. However, the systematic nature of the observed error suggests that it can be handled by a correction factor, although the implementation of the more accurate approach in a computer code is straightforward and recommended.

Since heavy-element compounds are subject to relativistic effects we have provided a thorough assessment of the reliability of different relativistic Hamiltonians for the calculation of the contact density and its gradient. Our extensive calibration studies show that the inclusion of spin-orbit coupling in an exact-two-component or four-component framework is mandatory if reasonably high accuracy is aimed for. Whereas scalar-relativistic Hamiltonians account for a major part of the relativistic effects with errors on the order of 10-12%, a non-relativistic ansatz yields meaningless contact density gradients for all three molecules reduced by up to a factor 6 compared to their four-component reference values. Calculations based on 2-component relativistic Hamiltonians show good performance, provided that the property operator is correctly transformed. Otherwise we find that the picture change errors are larger than relativistic effects. We furthermore point out the importance of considering electron correlation effects in the property evaluation by comparing results both from a choice of density functionals and high-level CCSD(T) calculations with Hartree-Fock data. As summarized recently in a related study on the isomer shift in mercury fluorides [START_REF] Knecht | Mössbauer spectroscopy for heavy elements: a relativistic benchmark study of mercury[END_REF] density functional theory (DFT) with various flavors of functionals exhibits a rather inconclusive overall performance relative to our CCSD(T) reference for the contact density gradient, covering not only a range of ±25% but showing also a strong system-dependency based on the quality of a result for a given functional type. One of the appealing feature of DFT is of course its low computational cost and it may thus, nevertheless, be applicable to further qualitative studies of field shift effects in other heavy-element diatomics. In this context, we have demonstrated how detailed information about changes in electronic structure upon changes in internuclear distances can be extracted from HF/DFT expectation values by projection analysis.

We finally evaluated the known experimental and theoretical electron density gradients and field shift factors for the three nuclei Tl, Pb, and Pt on the basis of our benchmark four-component CCSD(T) data. The comparison holds a surprising twist in the tale and may be best summarized as follows:

1. Our present four-component CCSD(T) electron density gradients and field shift factors are in very good agreement with the experimental predictions given by Schlembach and Tiemann [START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF] in 1982 with slightly larger deviations for Tl in TlI than Pb in PbTe.

2. As a consequence our results do not agree at all with the suggested scaling of the original data from 1982 which was published in 1985 by Knöckel et al. [START_REF] Knöckel | Molecular-beam-laser studies of the states X 1 Σ + and A0 + of PbS[END_REF] along with a revised formula to derive the field shift factor from a fit of rotational spectra data. We have carefully re-derived the appropriate formulas and can not see that such a scaling is justified, nor were any arguments provided for it in Ref. 36.

3. On the other hand, previous DFT(SAOP) predictions [START_REF] Cooke | The calculation of field shift effects in the rotational spectra of heavy metal-containing diatomic molecules using density functional theory: comparison with experiment for the Tl-halides and Pb-chalcogenides[END_REF][START_REF] Cooke | The rotational spectrum, nuclear field shift effects, 195 Pt nuclear spin-rotation constant, and electric dipole moment of PtSi[END_REF] of electron density gradients and field shift parameters of the respective heavy-metal centers in the present three molecules agreed quite well with the revised experimental data by Knöckel and co-workers [START_REF] Knöckel | Molecular-beam-laser studies of the states X 1 Σ + and A0 + of PbS[END_REF]. We conclusively show, though, that these calculations are plagued by a serious computational error. Having aimed at relativistic studies based on the ZORA Hamiltonian, these calculations turned out to be of non-relativistic nature as a comparison with our present non-relativistic DFT(SAOP) data unambiguously reveals. It so happens, and this is the twist in the tale, that the relativistic effect is of about the same order as the scaling factor proposed by Knöckel et al. [START_REF] Knöckel | Molecular-beam-laser studies of the states X 1 Σ + and A0 + of PbS[END_REF] to be applied to the original and correct data by Schlembach and Tiemann [START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF], which is probably why these inconsistencies were not discovered previously.

4. Our predicted field shift parameter for Pt in PtSi differs by a factor 7 from the experimental value reported by Cooke and co-workers [START_REF] Cooke | The rotational spectrum, nuclear field shift effects, 195 Pt nuclear spin-rotation constant, and electric dipole moment of PtSi[END_REF] which we currently attribute to be caused by a mistake in the experimental data fit.

We therefore recommend a careful re-examination of experiment.

We hope that our re-derivation of the key equations accompanied by the computational findings will put the theory of molecular isotopic field shift in rotational spectra back on solid footings and stimulate further experimental and theoretical work in this field. In future work we aim at investigating nuclear size effects of vibrational spectra, in which the second geometrical derivative of the effective/contact density comes into play. We will also continue projection analysis of these quantities to see what information about molecular electronic structure they contain. Ultimately it will be interesting to see if a combination of state of the art correlated electronic structure calculations combined with experiment can be used to extract information about changes in nuclear volume upon isotope substitution. (in Å-4 ) in TlI calculated at the 4-component DFT(PBE0) level using the TZ+2s1p basis. For each quantity contributions from s 1/2 and p 1/2 orbitals are given in the lower and upper triangle of the table, respectively.

D [1]T l 2p 1/2 3p 1/2 4p 1/2 5p 1/2 6p 1/2
1s 1/2 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 2p 1/2 2s 1/2 0.0000 0.0000 0.0000 0.0000 0.0000 0.0007 3p 1/2 3s 1/2 0.0000 0.0000 -0.0001 0.0000 -0.0006 -0.0050 4p

  . The final large-component basis sets (denoted TZ+2s1p (QZ+2s1p) in the following) thus read as [33s28p18d12f 1g] ([37s33p22d18f 4g1h]) for Pb and Tl, [29s22p16d4f 1g] ([34s28p19d6f 4g1h]) for Te and I, and [32s25p15d11f 4g1h] for Pt, respectively. As Si basis set the correlation-consistent Dunning basis set [90] of triple-ζ quality (ATZ) was likewise chosen in fully uncontracted form and augmented with diffuse functions.

  structuresTable I compiles the equilibrium bond distances R comp e and vibrational frequencies ω comp e that have been derived from our CCSD(T)/TZ+2s1p data as a by-product of the contact density calculations. Considering first the equilibrium bond distance R e , we find a very good agreement with the experimental data, with our values consistently being only 0.02 Å longer than their reference. Since our primary concern in the present study was not to reproduce experimental spectroscopic constants to highest precision we did not further pursue a basis set superposition error correction.Turning to the vibrational frequencies ω e , we observe a similar good performance of our CCSD(T)/TZ+2s1p data with the largest discrepancy being less than 2% for thallium iodide. Moreover, the present four-component molecular mean-field CCSD(T) values for platinum silicide evidently improve upon earlier theoretical estimates based on a CASPT2 study byBarysz and Pyykkö [91] who reported spectroscopic constants of R e = 2.1 Å and ω e = 531 cm -1

  4-component relativistic Dirac-Coulomb ( 4 DC) Hamiltonian obtained for TlI at the TZ+2s1p basis level, we see that going to the non-relativistic (NR) limit gives a reduction of -82%, or one order of magnitude, showing the importance of relativistic effects. Treating in turn exclusively scalar-relativistic effects at the spin-orbit free level (sf), results in a slight overshoot by +12% for ρ [1]T l 0 at the DFT(SAOP)/TZ+2s1p level, indicating the necessity to account for spin-orbit coupling contributions for this property irrespective of the closed-shell character of the molecular species. Consider next the ability of the eXact 2-Component (X2C) Hamiltonian to reproduce the 4 DC results. Of particular concern in 2-component relativistic calculations are picture change errors[42-44]. The 2-component relativistic oneelectron Hamiltonian 2c h is obtained by block diagonalization U of the parent 4-component Hamiltonian 4c h 2c h = U † 4c hU ++ .

[1]T l 0 in

 0 TlI by 283 % at the SAOP/TZ+2s1p level. Errors of similar size are observed for ρ

[1]T l 0 in

 0 TlI. Starting from the TZ value of 209.30 Å-4 , we see that adding tight 2s1p functions increases the value by 1.4 %, whereas going to the QZ+2s1p level has only a minor effect. In the following we therefore restrict attention to results obtained at the TZ+2s1p level.

[1]P t 0 in

 0 PtSi the hybrid functionals show the worst performance, and for ρ [1]P b 0 in PbTe the simplest functional LDA agrees best with the reference CCSD(T)/ 4c DC * * value.

  10 -17 m and δR P t(194,196) e =+3.7•10 -17 m, for TlI, PbTe and PtSi, respectively, which are exceedingly small numbers indeed.

Table

  II compiles for the Thallium atom individual orbital contributions to the effective density relative to the contact density, which were computed at the HF/TZ+2s1p level. The evaluation of effective densities has recently been made available in a development version of DIRAC10 and its implementation is described in full detail in Ref.35.

Two trends are clearly discernible from Table II: (i) the contact density is solely composed of contributions from s 1/2

TABLE I :

 I Spectroscopic constants for the 1 Σ + 0 ground state of 205 Tl 127 I, 208 Pb 130 Te, and 195 Pt 28 Si computed at the molecularmean-field 4 DC -CCSD(T)/TZ+2s1p level of theory (the reader may refer to the text for more details on the computational setup). Note that our values are not corrected for basis-set superposition errors.

	molecule R comp e	[ Å]	R exp e	[ Å]	ω comp e	[cm -1 ] ω exp e	[cm -1 ]
	TlI	2.838	2.8136 a	147.5	≈ 150 a
	PbTe	2.612	2.5949 a	211.0	211.9 a
	PtSi	2.083	2.0615 b (2.0629) c	544.9	549 d
	a Ref. [96]

b Ref.

[START_REF] Cooke | The rotational spectrum, nuclear field shift effects, 195 Pt nuclear spin-rotation constant, and electric dipole moment of PtSi[END_REF]

c Ref.

[START_REF] Shao | Resonant two-photon ionization spectroscopy of jet-cooled PtSi[END_REF]

d Ref.

[START_REF] Paul | Cavity ringdown laser absorption spectroscopy and time-of-flight mass spectroscopy of jet cooled platinum silicides[END_REF]

TABLE II :

 II Atomic matrix elements for the Thallium atom (HF/TZ+2s1p), comparing the contact ρe(0) and effective ρe (number) density . All values are in atomic units a -3 0 .

		ρe(0) ρe -ρe(0)
	1s 1/2 2112536.88 -215363.37
	2s 1/2 322388.19 -33035.87
	3s 1/2	74350.56 -7625.34
	4s 1/2	18799.34 -1928.44
	5s 1/2	3695.82	-379.14
	6s 1/2	371.46	-38.11
	2p 1/2	24695.53 -2449.09
	2p 3/2	0	2 x 0.56
	3p 1/2	6395.44	-634.80
	3p 3/2	0	2 x 0.16
	4p 1/2	1598.33	-158.67
	4p 3/2	0	2 x 0.08
	5p 1/2	280.27	-27.82
	5p 3/2	0	2 x 0.01
	6p 1/2	2.46	-0.24
	6p 3/2	0	< 0.01
	Total 2565114.28 -261639.27

TABLE III :

 III Effective and contact electron density gradients, respectively, at the nuclei Tl, Pb, and Pt calculated at the fourcomponent Dirac-Coulomb level using the TZ+2s1p basis. The derivatives (in Å-4 ) are taken at the respective experimental geometries of TlI, PbTe, and PtSi.

	Method	ρ	[1]X 0	ρ[1]X	∆(ρ [1]X -ρ	[1]X 0	)
				TlI
	HF	212.15 190.29	-10.3%
	DFT/PBE	114.53 102.67	-10.4%
	DFT/SAOP 119.96 107.60	-10.3%
				PbTe
	HF	240.18 215.00	-10.5%
	DFT/PBE	155.74 139.29	-10.6%
				PtSi
	HF	-1183.83 -1070.52	-9.6%
	DFT/PBE	-656.39 -593.34	-9.6%

TABLE IV :

 IV Projection analysis of the contact density gradient ρ

	[1]T l 0	(in Å-4 ) in TlI calculated at the 4-component DFT(PBE0)
	level using the TZ+2s1p basis.
	Intra-atomic contribution Tl total		103.47
	diagonal (p = q)	-135.91
	hybridization (p = q) 243.99
	I total		0.00
	Inter-atomic contribution		0.30
	Polarization contribution		32.63

TABLE V :

 V Detailed analysis of individual orbital contributions to geometric derivative of atomic density matrix D

		[1]T l qp	and
	contact density gradient ρ	[1]T l 0

  -0.0001 -0.0004 -0.0013 0.0040 -0.0079 -0.0722 6p 1/2 6s 1/2 0.0001 0.0007 0.0024 -0.0077 0.0332 0.0222 1s 1/2 2s 1/2 3s 1/2 4s 1/2 5s 1/2 6s 1/2 -17.67 43.59 -74.25 118.29 -104.54 -3.90 6p 1/2 6s 1/2 10.54 -26.03 44.61 -72.71 142.75 30.86 1s 1/2 2s 1/2 3s 1/2 4s 1/2 5s 1/2 6s 1/2

						1/2
	4s 1/2 0.0000 0.0000 0.0003 -0.0005 -0.0028 -0.0443 5p 1/2
	5s 1/2 ρ	[1]T l 0
		2p 1/2 3p 1/2 4p 1/2 5p 1/2 6p 1/2
	1s 1/2 -0.75	0.00	0.00	0.00	0.00	0.00 2p 1/2
	2s 1/2 3.69	-4.55 -0.01	0.00	0.00	-0.72 3p 1/2
	3s 1/2 0.00	15.49 -13.19 -0.16	1.35	2.70 4p 1/2
	4s 1/2 0.00	10.00 42.01 -33.45 -2.84 -10.34 5p 1/2
	5s 1/2

TABLE VI :

 VI Contact density gradient ρ[1]T l 0 (in Å-4 ) of the TlI molecule computed at various level of theory. The derivative is calculated at the experimental equilibrium interatomic distance R exp e = 2.8136 Å[START_REF] Huber | Constants of Diatomic Molecules[END_REF].

		TlI
	method	ρ e evaluation Hamiltonian basis set	ρ	[1]T l 0
	HF	ψ † ψ	4c DC	TZ	209.30
	HF	ψ † ψ	4c DC	TZ+2s1p 212.15
	HF	ψ † ψ	4c DC	QZ+2s1p 212.07
	DFT/SAOP	ψ † ψ	NR	TZ+2s1p 21.09
	DFT/SAOP	ψ † ψ	X2C	TZ+2s1p 461.57
	DFT/SAOP	ρ Tl	X2C	TZ+2s1p 120.50
	DFT/SAOP	ψ † ψ	4c DC-sf	TZ+2s1p 134.48
	DFT/SAOP	ψ † ψ	4c DC	TZ+2s1p 119.96
	DFT/SAOP	ψ † ψ	4c DC	QZ+2s1p 119.75
	DFT/LDA	ψ † ψ	4c DC	TZ+2s1p 106.64
	DFT/PBE	ψ † ψ	4c DC	TZ+2s1p 114.53
	DFT/BLYP	ψ † ψ	4c DC	TZ+2s1p 107.09
	DFT/PBE0	ψ † ψ	4c DC	TZ+2s1p 140.63
	DFT/B3LYP	ψ † ψ	4c DC	TZ+2s1p 127.46
	DFT/CAMB3LYP	ψ † ψ	4c DC	TZ+2s1p 154.75
	CCSD(T)	ff	4c DC	TZ+2s1p 142.47
	CCSD(T)	ff	4c DC * *	TZ+2s1p 142.43
	CCSD(T)	ff a	4c DC * *	TZ+2s1p 142.26
	CCSD(T)	ff	4c DC * *	QZ+2s1p 142.19
	Exp.[27]			120.6(38)

a energy threshold for virtual spinor: 134 Hartree

TABLE VII :

 VII Contact density grdient ρ [1]P b 0 (in Å-4 ) of the PbTe molecule computed at various level of theory. The derivative is calculated at the experimental equilibrium interatomic distance R exp e = 2.5949 Å[96].

		PbTe
	method	ρ e evaluation Hamiltonian basis set	ρ	[1]P b 0
	HF	ψ † ψ	4c DC	TZ	236.84
	HF	ψ † ψ	4c DC	TZ+2s1p 240.18
	DFT/SAOP	ψ † ψ	NR	TZ+2s1p 19.64
	DFT/SAOP	ψ † ψ	X2C	TZ+2s1p 585.06
	DFT/SAOP	ρ Pb	X2C	TZ+2s1p 162.19
	DFT/SAOP	ψ † ψ	4c DC-sf	TZ+2s1p 174.59
	DFT/SAOP	ψ † ψ	4c DC	TZ+2s1p 161.41
	DFT/LDA	ψ † ψ	4c DC	TZ+2s1p 148.16
	DFT/PBE	ψ † ψ	4c DC	TZ+2s1p 155.74
	DFT/BLYP	ψ † ψ	4c DC	TZ+2s1p 148.33
	DFT/PBE0	ψ † ψ	4c DC	TZ+2s1p 181.72
	DFT/B3LYP	ψ † ψ	4c DC	TZ+2s1p 169.96
	DFT/CAMB3LYP	ψ † ψ	4c DC	TZ+2s1p 197.55
	CCSD(T)	ff	4c DC * *	TZ+2s1p 163.19
	Exp.[27]				148(11)

TABLE VIII :

 VIII Contact density gradient ρ[1]P t 0 (in Å-4 ) of the PtSi molecule computed at various level of theory. The derivative is calculated at the experimental equilibrium interatomic distance R exp e = 2.0615 Å[START_REF] Cooke | The rotational spectrum, nuclear field shift effects, 195 Pt nuclear spin-rotation constant, and electric dipole moment of PtSi[END_REF].

			PtSi
	method	ρ e evaluation Hamiltonian basis set	ρ	[1]P t 0
	HF	ψ † ψ	4c DC	TZ	-1168.75
	HF	ψ † ψ	4c DC	TZ+2s1p	-1183.83
	DFT/SAOP	ψ † ψ	NR	TZ+2s1p	-135.64
	DFT/SAOP	ψ † ψ	X2C-sf	TZ+2s1p	-2390.42
	DFT/SAOP	ρ Pt	X2C-sf	TZ+2s1p	-732.09
	DFT/SAOP	ψ † ψ	4c DC-sf	TZ+2s1p	-729.65
	DFT/SAOP	ψ † ψ	4c DC	TZ+2s1p	-738.54
	DFT/LDA	ψ † ψ	4c DC	TZ+2s1p	-643.41
	DFT/PBE	ψ † ψ	4c DC	TZ+2s1p	-656.39
	DFT/BLYP	ψ † ψ	4c DC	TZ+2s1p	-625.41
	DFT/PBE0	ψ † ψ	4c DC	TZ+2s1p	-789.17
	DFT/B3LYP	ψ † ψ	4c DC	TZ+2s1p	-728.07
	DFT/B3LYP a	ψ † ψ	DKH(2,0) ANO-RCC -2182.95
	DFT/B3LYP a	ψ † ψ	DKH(8,8) ANO-RCC -749.34
	DFT/CAMB3LYP	ψ † ψ	4c DC	TZ+2s1p	-837.70
	CCSD(T)	ff	4c DC	TZ+2s1p	-599.39
	CCSD(T)	ff	4c DC * *	TZ+2s1p	-599.83

a M. Reiher and R. Mastalerz, private communication.

TABLE IX :

 IX Comparison of both the calculated and, where available, measured electron density gradient ρ[1]X (in Å-4 ) at the Tl, Pb, and Pt nuclei and the field shift parameter V X (in 10 4 Å-2 ). Besides our best theoretical estimates computed at the fourcomponent molecular mean-field CCSD(T)/TZ+2s1p level ( 4c DC ) we also list non-relativistic (NR) DFT/SAOP/TZ+2s1p data.

	method	Hamiltonian ρ [1]X	V X
			TlI
	experiment[A] a		120.6(38)	3.20(10)
	experiment[B] b		12.06(38)
	DFT/SAOP	NR	21.09	0.57
	CCSD(T)	4c DC	142.26	3.81
	DFT/SAOP c	(ZORA)	21.1	0.61 d
			PbTe
	experiment a		148(11)	2.12(16)
	experiment b		14.8(11)
	DFT/SAOP	NR	19.64	0.28
	CCSD(T)	4c DC	163.19	2.36
	DFT/SAOP c	(ZORA)	21.1	0.33 d
			PtSi
	experiment[C] e			-0.72(12)
	DFT/SAOP	NR	-135.64	-1.14
	CCSD(T)	4c DC	-599.83	-5.05
	DFT/SAOP e	(ZORA)	-136.5	-1.10

a original result published in Ref.

[START_REF] Schlembach | Isotopic field shift of the rotational energy of the Pb-chalcogenides and Tl-halides[END_REF]

in 1982 b revised data published in Ref. 36 in 1985 c Ref. 38; a QZ4P basis set was used in the DFT calculations d Ref. 40 e Ref. 39; a QZ4P basis set was used in the DFT calculations

Acknowledgments

We dedicate this paper to Debashis Mukherjee for his profound insight into the electron correlation problem and for his many important contributions on how to treat it. S. K. gratefully acknowledges postdoctoral research grants from l' Université de Strasbourg (UdS) and the Natural Science Foundation (FNU) of the Danish Agency for Science, Technology and Innovation. We would like to thank Markus Reiher and Remigius Mastalerz at ETH Zürich for providing us their unpublished results on the molecular field shift effect in PtSi. We would also like to thank Radovan Bast (Toulouse) for graphical assistance. This work has been supported through ample computing time at the supercomputer centers of UdS and the Danish Center for Scientific Computing at SDU Odense. This work is part of the project NCPCHEM funded by the Agence Nationale de la Recherche (ANR, France).