Adel Almoukhalalati

Avijit Shee

Trond Saue

Nuclear size effects in vibrational spectra

We present a theoretical study of nuclear volume in the rovibrational spectra of diatomic molecules which is an extension of a previous study restricted to rotational spectra [Chem. Phys. 401 (2012) 103]. We provide a new derivation for the electron-nucleus electrostatic interaction energy which is basically independent of the choice of model for the nuclear charge distribution.

Starting from this expression we derive expressions for the electronic, rotational and vibrational field shift parameters in terms of the effective electron density and its first-and second derivatives with respect to internuclear distance. The effective density is often approximated by the contact density, but we demonstrate that this leads to errors on the order of 10% and is furthermore not necessary since the contact and effective densities can be obtained at the same computational cost. We calculate the field shift parameters at the 4-component relativistic coupled-cluster singles-and-doubles level and find that our results confirm the experimental findings of Tiemann and co-workers [Chem. Phys, 68 (1982) 21, 1982, Ber. Bunsenges. Phys. Chem., 86 (1982) 821], whereas we find no theoretical justification for a scaling factor introduced in later work [Chem. Phys. 93 (1985) 349]. For lead sulfide we study the effective density as a function of internuclear distance and find a minimum some 0.2 Å inside the equilibrium bond distance. We also discuss Bigeleisen-Goeppert-Mayer theory of isotope fractionation in light of our results.

Introduction

The concept of isotopes grew out of the work by Fajans and Soddy on radioactive decay chains. Soddy coined the name isotope for elements having different atomic mass, but identical chemical properties, so that they should occupy the same place (ίσ ος+τόπος) in the periodic table 1 . There are in fact subtle chemical differences between isotopes of the same element, affecting both chemical rates and equilibria, which has been exploited for their separation and which explains the evolution of relative isotope abundances in geochemical cycles [2][START_REF]Isotope Effects In Chemistry and Biology[END_REF][START_REF] Wolfsberg | Isotope Effects: in the Chemical[END_REF][START_REF] Eiler | [END_REF] . The first successful observation of isotope effects in spectroscopy was a shift of 0.0043 Å of the 4058 Å line of neutral "radio-lead" with respect to ordinary lead reported by Aronberg in 1918. 6 Although small, it was clearly larger than what could be expected from atomic mass differences alone. Bohr in 1922 correctly suggested that the shift was due to "a slight difference in the field of force surrounding the nucleus, arising from the difference in the internal structure of the lead isotopes". [7][START_REF] King | Isotope Shifts in Atomic Spectra[END_REF][START_REF] Kragh | [END_REF] A crude, non-relativistic estimate of the energy shift with a change of radius of the nucleus, described as a charged conducting sphere, was given by Bartlett a Laboratoire de Chimie et Physique Quantiques, UMR 5626 CNRS -Université Toulouse III-Paul Sabatier 118 route de Narbonne, F-31062 Toulouse, France. Fax:+33/(0)561556031 ; Tel: +33/(0)561556065; E-mail: trond.saue@irsamc.upstlse.fr in 1931. 10 Racah derived an expression for the energy shift in terms of the non-relativistic density at the nuclear origin, but with relativistic corrections. 11 The formula is evidently only valid for s 1/2 -orbitals since only non-relativistic s-orbitals have density at the nucleus. Rosenthal and Breit confirmed the formula of Racah and extended the analysis to relativistic orbitals other than s 1/2 , in particular p 1/2 orbitals that do have density at the nucleus (see for instance the analysis in ref. 12). An estimate of the importance of p 1/2 contributions, based on the work of Borch 13 and Bodmer [START_REF] Bodmer | Proceedings of the Physical Society. Section A[END_REF][START_REF] Bodmer | [END_REF] , was provided by Fradkin 16 , who concluded that for mercury the np 1/2 -contribution would be about 10% of the ns 1/2 one. Interestingly, none of the above theoretical studies cites the pioneer experimental paper by Aronberg, rather the subsequent experimental observations by Schüler, Kopfermann and co-workers 17,18 . Historical reviews and further discussion of the isotope field shift in atomic spectra are found in refs. 19

and 8.

A theoretical framework for the study of isotope effects on chemical equilibria was provided by the seminal work of Urey and co-workers [START_REF] Urey | [END_REF][21][22] as well as Bigeleisen and Goeppert-Mayer 23 and later extended to reaction rates by Bigeleisen. 24 However, these formulations only consider the isotopic mass effect and is therefore primarily appropriate for light elements. In 1996 Bigeleisen took into account the isotope field effect on chemical equilibria 25 , based on the observation by Fujii et al. of an anomalous isotope effect of 235 U in the U(IV)-U(VI) exchange reaction. 26 Theoreti-cal studies along these lines have subsequently been reported by Schauble, Abe, Liu, Wiederhold and co-workers. [27][28][29][30][31][32][START_REF] Schauble | Proc. Natl. Acad. Sci[END_REF][START_REF] Nemoto | [END_REF][35] A review on the nuclear field shift effect in chemical exchange reactions has been given by Fujii et al.. 36 In the present contribution we focus on the nuclear volume effect (NVE) in the rovibrational spectra of diatomic molecules. Such effects were proposed by Knöckel and Tiemann 37 in 1982 after anomalously large adiabatic correction factors appeared in the Dunham analysis of the rotational spectrum of diatomic molecules containing lead and thallium. 38 The fact that the nuclear volume effect was first observed in these diatomic molecules is perhaps not surprising since lead and thallium are the heaviest elements having more than one stable isotope. Schlembach and Tiemann subsequently derived expressions for the modification of rotational and force constants due to nuclear volume effects 39 and re-examined the spectra reported in ref. 38. Curiously, the Tiemann group in 1995 revised their results by a factor ten 40 , results that were corroborated by calculations (and experimental work) by Cooke and co-workers. [41][42][43] However, Knecht and Saue 44 in 2012 provided an independent evaluation of both theoretical and experimental studies of the isotopic field shift in rotational spectra using high-level relativistic 4-component electronic structure methods. They concluded that there was no justification for the scaling factor introduced by Tiemann and co-workers and could show that the calculations of Cooke et al. were unwittingly non-relativistic and therefore flawed. On the other hand, the results of Knecht and Saue were generally in good agreement with the experimental predictions given by Schlembach and Tiemann in 1982. 39 In the present work we extend the previous study of isotopic field shift by Knecht and Saue to vibrational spectra. We provide a new derivation of the nuclear volume effect on spectroscopic constants that is independent of the chosen model of the nuclear charge distribution. We furthermore report 4-component Hartree-Fock (HF) and Density Functional Theory (DFT) calculations of Pb-chalcogenides and Tl-halides and benchmark them by our newly implemented module for the analytic calculation of expectation values at the relativistic coupled cluster level. [START_REF] Shee | Analytic Gradient at the 4component Relativistic Coupled Cluster Level with Inclusion of Spin-Orbit Coupling[END_REF] We also investigate the density at the heavy nucleus as a function of equilibrium bond length and discuss the impact of our results on Bigeleisen-Goeppert-Mayer theory of isotopic fractionation.

Theory

The Dunham analysis of the rovibrational spectra of diatomic molecules is based on the solution of the radial nuclear Schrödinger equation for the rovibrational levels of the electronic state (the rotating vibrator problem) of a (closed-shell) diatomic molecule AB -h2 2µ

d 2 dR 2 + h2 J (J + 1) 2µR 2 +V (R) Ψ ν,J (R) = E ν,J Ψ ν,J (R) (1)
provided by Dunham 46 in 1932 using the semiclassical Jeffreys-Wentzel-Brillouin-Kramers (JWBK) approach. [47][48][49][50] In the above equation R is the internuclear distance, µ the reduced mass and ν and J the vibrational and rotational quantum numbers, respectively. Furthermore V (R) is the electronic potential, corresponding to the total energy obtained upon solving the electronic problem within the Born-Oppenheimer ("clamped nuclei") approximation. The Dunham expansion of the rovibrational energy is given by

E ν,J = h ∑ k,l=0 Y kl ν + 1 2 k [J (J + 1)] l (2)
where the Dunham coefficients Y kl are given in units of frequency.

To lowest order in ν and J we have

Y 10 = ν e 1 + B 2 e 4ν 2 e (. . .) ν e ; 2πν e = k e µ (3)
Y 01 = B e 1 + B 2 e 2ν 2 e (. . .) B e = h 4π µR 2 e (4)
where ν e , k e , B e and R e are the equilibrium harmonic frequency, force constant, rotational constant and bond length, respectively. The above Dunham coefficients have an explicit dependence on reduced mass µ, specific to each isotopomer, which can be generalized to

Y kl = µ -( k /2+l) U kl (5)
where U kl are isotope-independent coefficients. However, further isotope dependence of the Dunham coefficients arises by i) switching from the clamped nucleus to the adiabatic approximation so that the internuclear potential E el (R) becomes specific for each isotopomer, ii) by taking into account the coupling of different electronic states upon breakdown of the Born-Oppenheimer approximation and iii) by correcting for the approximate nature of the JWBK method. Ross et al., in an experimental study of isotopomers of CO, 51 fitted their experimental data to the extended expression

Y kl = µ -( k 2 +l) U kl 1 + m e A kl M A + B kl M B ( 6
)
where M X are atomic masses and X kl are mass-independent correction factors for each atom A and B, and for which Bunker 52 and Watson 53 provided theoretical justification. The resulting energy expression was criticized by Le Roy, in part because of the non-linear appearance of fitting coefficients and that the precise determination of the isotope-independent coefficients U kl requires data from several isotopomers. 54 Le Roy suggested to use the original Dunham expansion, eqn (2), for a reference isotopomer and then obtain corresponding Dunham coefficients for other isotopomers by fitting correction factors. In the present contribution we choose a middle way in that we recognize that the coefficients U kl of eqn (6) have no isotopic mass dependence, but do depend on nuclear size (and possibly shape, which we shall not address here). In the following we shall therefore consider the change of coefficients U 10 ν e and U 10 µB e (7) upon isotopic substitution A → A in a reference isotopomer AB.

In order to avoid a profusion of subscripts we in the following let spectroscopic constants R e , k e , B e and f 3 be those of the reference isotopomer.

The nuclear volume effect arises from the modification of the electrostatic interaction energy between electrons and nuclei upon isotope substitution. We shall therefore start by deriving a convenient expression for this interaction energy. Our approach is inspired by Otten, [START_REF] Otten | The Gregory Breit Centennial Symposium[END_REF] but do not assume the electron density to be constant in the nuclear region. A particular attractive feature of the derivation is that it makes minimal assumptions concerning the nuclear charge distribution. We start by writing the electrostatic interaction energy as

E en = -e ∑ A ρ e (r e )φ A (r e ) d 3 r e , (8)
where we, in line with previous work, let ρ e denote electronic number density, rather than charge density. The scalar potential associated with our target nucleus A is expressed as

φ A (r e ) = 1 4πε 0 r n =r A r n =0 ρ A (r n ) |r e -r n | d 3 r n , (9)
where for convenience the origin of our coordinate system coincide with the center of nucleus A and where r A is an effective nuclear radius outside which the nuclear charge distribution ρ A is zero. From elementary electrodynamics it follows that the scalar potential outside this radius is that of a equivalent point charge.

The interaction energy associated with nucleus A can thereby be rearranged as (r e ) d 3 r e (10) The second term is the interaction energy associated with a point nuclear charge, but it should be kept in mind that the electron density ρ e is the one associated with a molecular system of finite nuclear charge distributions. We now focus on the first term and start by expanding the inverse electron-nucleus distance of eqn (9) in Legendre polynomials. At this point we will assume the nuclear charge distribution to be spherically symmetric. It is then reasonable to assume that the electron density is spherically symmetric as well in the nuclear region, that is, within the radius r A , even for a molecular system. This has two consequences for the further derivation: i) The Legendre expansion will truncate to lowest order (L = 0) and ii) Since the potential difference

E A en = -e
φ A (r e ) -φ point A
(r e ) is negative (and does not change sign) on the integration interval [0, r A ] we can invoke the first mean value theorem for definite integrals and write

Further manipulation then leads to the final expression

E A en = -e ρ e (r e )φ point A (r e ) d 3 r e + Z A e 2 6ε 0 ρA e r 2 n A (13)
where appears the mean square radius of the nuclear charge distribution

r 2 n A = ρ A (r n )r 2 n d 3 r n ρ A (r n )d 3 r n ( 14
)
which can for instance be obtained from scattering experiments. [55][56][57] It should be noted that the final expression makes no reference to the effective nuclear radius r A . In the derivation it provides a radius at which the finite and point charge potentials coincide and furthermore furnishes a closed integration interval for the application of the mean value theorem.

Let us now consider the energy shift associated with the change of nuclear volume, as gauged by the mean square radius r 2 n A , upon isotopic substitution A → A . Following Filatov, 58,59 we can express the energy shift by a first-order Taylor expansion, that is

δ E AA en = E A en -E A en = V A 00 δ r 2 n AA ; V A 00 = dE X en d r 2 n r 2 n A = Z A e 2 6ε 0 ρA e , (15)
where we use the notation δ r 2 n AA = r 2 n Ar 2 n A . Notice that this expression is valid not only for isotope field shifts, but also for the isomer shift appearing in Mössbauer spectroscopy, where the change of nuclear radius occurs upon change of the state of the nucleus induced by absorption or emission of gamma radiation. [START_REF]Mössbauer Spectroscopy and Transition Metal Chemistry[END_REF][START_REF] Gütlich | [END_REF] In both cases the effective density ρA e is typically approximated by the contact density ρ A 0 , that is, the electron (number) density at the nuclear origin of atom A. For heavy elements, though, this is a quite severe approximation that may lead to errors, albeit systematic for each element, on the order of 10%. 12,44,62 It is also an unnecessary approximation since the effective density can be obtained at the same computational cost as the contact density. In fact, eqn (15) serves as an operational definition of the effective density with respect to any chosen model for the nuclear charge distribution. In the present work we employ a Gaussian model

ρ G n (r) = ρ G 0 exp -ηr 2 ; ρ G 0 = Ze η π 3/2 (16)
where the exponent η is chosen such that the mean square radius

r 2 = 3 2η (17)
reproduces the empirical formula 63,64

r 2 n 1/2 = 0.836 A 1/3 + 0.570 f m ( 18
)
where A is the atomic mass number 64 . eqn 15 then leads to the following definition of the effective density

ρe = 1 Ze ρ e (r) ρ G n (r) d 3 r = 1 Ze 0 ρ G n 0 . ( 19
)
An alternative to the use of the effective density is to Taylorexpand the electronic density ρ e (r) about the contact density, but this will introduce higher-order momenta of the nuclear charge distribution 62,65 and lead to a more cumbersome formalism. An approximation inherent in eqn (15) is that it is assumed that the electronic density is assumed identical for both isotopes. The error introduced by this approximation was found to be on the order of 0.2 -0.4 % by Fricke and Waber, 62 which is an acceptable error for most purposes.

We are now in a position to consider the nuclear volume effect on spectroscopic properties, in particular the equilibrium force constant k e and the equilibrium rotational constant B e . The NVE arises in part because of a minute change of the internuclear distance upon isotope substitution. In order to derive the desired expressions we shall employ variational perturbation theory, as defined by Helgaker and Jørgensen: [START_REF] Helgaker | Methods in Computational Molecular Physics[END_REF] We consider the electronic energy as a function of both internuclear distance R and nuclear mean square radius r 2 n of atom A

E el ≡ E el R, r 2 n ( 20
)
and assume that the electronic energy has been optimized with respect to bond length for any r 2 n , that is

dE el dR r 2 n = 0. (21)
In this manner the equilibrium distance becomes a function of nuclear mean square radius R e ≡ R e r 2 n . Since the variational condition, eqn (21), is valid for any nuclear mean square radius, its derivative with respect to nuclear mean square radius is zero as well

d d r 2 n dE el dR r 2 n = 0 (22) which implies ∂ ∂ r 2 n dE el dR + d 2 E el dR 2 dR d r 2 n r 2 n = 0. ( 23
)
Upon rearrangement we find that the first-order change of equilibrium distance is

δ R AA e = dR d r 2 n r 2 n A δ r 2 n AA = - Z A e 2 6ε 0 k e ρA[1] e δ r 2 n AA ( 24
)
where appears the equilibrium force constant k e as well as the first geometric derivative of the effective density

ρA[1] e = d ρe dR A . ( 25
)
The first-order change in the rotational coefficient U 01 is

δU AA 01 = dU 01 d r 2 n r 2 n A δ r 2 n AA (26) = U 01 - 2 R e dR d r 2 n r 2 n A δ r 2 n AA = U 01 V A 01 δ r 2 n AA
where appears the field-shift parameter

V A 01 = Z A e 2 3ε 0 1 k e R e
ρA [1] e .

The expression for V A 01 is identical to Eq. ( 25) of both Schlembach and Tiemann 39 and Knecht and Saue, 44 as well as Eq. ( 5) of Tiemann et al.. [START_REF] Tiemann | [END_REF] However, in the 1985 paper by Knöckel et al. 40 Eq.( 7) the parameter is scaled by almost a factor ten

V A 01 → π 2 V A 01 , (28)
for which we can see no justification.

The first-order change in the equilibrium force constant k e with respect to nuclear mean square radius is

δ k AA e = dk e d r 2 n r 2 n A δ r 2 n AA = ∂ k e ∂ r 2 n + ∂ k e ∂ R dR d r 2 n r 2 n A δ r 2 n AA = Ze 2 6ε 0 ρA[2] e - f 3 k e ρA[1] e δ r 2 n AA (29)
where appears the first-and second geometric derivatives of the effective number density ρA e as well as the cubic force constant f 3 . The final expression is equivalent to Eq. ( 18) of Schlembach and Tiemann 39 when one takes into account that the Dunham potential coefficient

a 1 = f 3 R e 3k e (30)
and that they approximate the effective density by the contact density. Schlembach and Tiemann consider the change of the force constant when going from a point to a finite nucleus, whereas we consider the change upon isotopic substitution A → A . A problem with the former approach is that the contact density is singular in the relativistic domain.

The vibrational field-shift parameter V A 10 is found as

δU AA 10 = U 10 δ k AA e 2k e = U 10 V A 10 δ r 2 n AA ; V A 10 = Ze 2 12ε 0 k e ρA[2] e - f 3 k e ρA[1] e (31
) and is equivalent to Eq. ( 6) of Tiemann et al.. [START_REF] Tiemann | [END_REF] It also agrees with Eq. ( 8) of Knöckel et al. 40 Eq.( 7), but only if our expression, eqn (27) for the rotational field shift parameter V A 01 is taken.

Computational details

All calculations reported in this paper have been performed using the DIRAC package [START_REF] Dirac | a relativistic ab initio electronic structure program, Release DIRAC15[END_REF] and are based on the 4-component relativistic Dirac-Coulomb Hamiltonian as well as uncontracted basis sets of Gaussian-type orbitals (GTO). The default Hamiltonian of the DIRAC package uses a energy correction to avoid the explicit calculation of two-electron integrals containing small component basis functions only. [START_REF] Visscher | [END_REF] However, calibration studies showed that this gave errors on the order of 10% for the gradient and Hessian of the contact and effective densities. We have therefore retained the (SS|SS) class of integrals in all calculations (using the DOSSSS flag). On the other hand, inclusion of the Gaunt two-electron interaction modified results on the order of 1% and was therefore neglected. For the elements O, S, Se, Te, F, Cl, Br and I we used Dyall's triple zeta basis sets including core-valence correlation functions (denoted dyall.acv3z). [70][71][72] For the heavy atoms Pb and Tl, in light of previous studies, 12,73 two more tight s and one tight p function were added to the dyall.acv3z basis. The final largecomponent basis set, denoted TZ+2s1p, are [33s28p18d12f1g] and [37s33p22d18f4g1h] for Pb and Tl, respectively. Small component basis functions were generated from the large component functions using the restricted kinetic balance condition. 74 As explained in the previous section, the nuclear charge distribution was modeled by a Gaussian function. 64 Benchmark values of the effective and contact densities at the Pb-and Tl-centers were obtained at the coupled-cluster singlesand-doubles (CCSD) level using a newly developed extension of the RELCCSD module of the DIRAC package, 75,76 allowing the calculation of CCSD expectation values in an orbital-unrelaxed formalism. [START_REF] Shee | Analytic Gradient at the 4component Relativistic Coupled Cluster Level with Inclusion of Spin-Orbit Coupling[END_REF] For PbO, PbS, TlF and TlCl 20 electrons were correlated, corresponding to the 5d6s6p shells of Pb/Tl and the nsnp shells of the chalcogens/halogens. For PbSe, PbTe, TlBr and TlI the (n-1)d shells of the chalcogens/halogens were included as well, giving a total of 30 correlated electrons. An energy cutoff of 40 E h in the virtual space was introduced, but care was taken to include all core and valence correlation as well as valence dipole polarization functions.

We also carried out 4-component relativistic density functional theory (DFT) calculations, 77 using a wide variety of functionals: PBE, 78 PBE0, [79][80][81] BLYP, 82-84 B3LYP 82,85,86 and CAMB3LYP. 87 All DFT calculations were carried out with an "ultrafine" grid for numerical integration to guarantee well converged results in the exchange-correlation functional evaluation. In practice, for this purpose, we used the .ULTRAFINE option in the DIRAC package which means that radial integration is performed using the scheme proposed by Lindh et al. 88 with a convergence threshold of 2 × 10 -15 , while the angular integration is carried out using the Lebedev scheme with L max = 64 . Geometric first-and second derivatives of the contact and effective densities were obtained by numerical differentiation using a three-point stencil

ρ (x) = ρ(x + h) -ρ (x -h) 2h ; ρ (x) = ρ(x + h) -2ρ(x) + ρ(x -h) h 2 (
32) with a step size h=0.01 Å. 44,[START_REF] Hildebrand | Introduction to numerical analysis[END_REF] For selected molecules calibration studies using five-and seven-point stencils and other step sizes were carried out to assure that the numerical scheme was stable. We also calculated spectroscopic constants (equilibrium bond length r e , harmonic frequencies ω e and anharmonicities ω e x e ) based on a 6th order polynomial fit of 9 energy points and masses of the most abundant isotopes, using the auxiliary program TWOFIT which is included in the DIRAC package. [START_REF] Saue | program TWOFIT[END_REF][START_REF] Fossgaard | [END_REF] 4

Results and discussion

In Tables 1 and2 we present calculated geometric derivatives of contact and effective densities at the position of the heavy nucleus in PbX (X=O, S, Se, Te) and TlY (Y=F, Cl, Br, I), respectively. We have also included calculated and experimental spectroscopic constants (r e , ω e , ω e x e ) to further indicate the performance of each method.

Methodological considerations

We start by considering the performance of SCF methods, that is, Hartree-Fock and Kohn-Sham, using our CCSD data as reference. We have calculated relative errors of calculated geometric derivatives of effective densities with respect to CCSD and present in Table 3 corresponding mean relative errors ∆ and standard deviations ∆ std . For the gradient of the contact density ρ[1]X e the CAMB3LYP functional has the overall best performance in terms of accuracy ( ∆) and precision ∆ std . However, for the Hessian ρ[2]X e of the effective density the same functional is not at all reliable.

In fact, for the lead compounds the DFT functionals show large scattering. One might suspect errors due to numerical evaluation of the functionals, despite the use of ultrafine grids, but HF shows the same degradation of performance when going from the gradient to the Hessian of the effective density. We have also carefully checked our procedure for numerical geometric derivatives, but have not detected any instabilities. What can be noted, though, is that CAMB3LYP overall provides excellent agreement with experiment for equilibrium bond lengths and harmonic constants, but not for anharmonicites, and only the wave function based methods HF and CCSD show reasonable standard deviations for anharmonicites. In conclusion we find that whereas some DFT functionals, in particular CAMB3LYP and PBE0, provide reliable gradients of effective densities ρ [1] e compared to CCSD, this no longer holds true for Hessians of effective densities ρ [2] e . We next consider the effect of approximating the effective density ρX e by the contact density ρ X 0 . For all methods we find that replacing the contact density by the effective density reduces gradients systematically by about 10.5 %; for CCSD we obtain for instance -10.40(11) %. However, for SCF methods the corresponding deviation for Hessians varies substantially, whereas for CCSD we obtain -10.13(54) %, which is more reasonable and seems to indicate that the CCSD values are reliable.

Effective density as a function of internuclear distance

Following Tiemann et al. [START_REF] Tiemann | [END_REF] we plot in Figure 1 our CCSD values for the gradient of the heavy-atom effective density ρ [1] e as a function of reported experimental equilibrium bond lengths. A rather striking linear relationship becomes apparent, probably indicating similar electronic structure of the title compounds. Tiemann and co-workers further speculated that they could use this relationship to estimate a value for the Hessian of the contact density, or, more correctly, the effective density. Linear regression of the data in Figure 1 gives a slope of -237 Å-5 . This is indeed very close to our CCSD value of ρTl [2] = -238 Å-5 obtained for TlCl. However, further inspection of Tables 1 and2 shows that the Hes-sian of the effective density span the range -414 Å-5 to 210 Å-5 . It is furthermore positive for the lead chalcogenides, but negative for thallium halides. An example of the effective density as a function of internuclear distance R is given in Figure 2. We plot the effective density ρPb e of the PbS molecule relative to its value at R = 4.0 Å. The values have been obtained at the Kramers-restricted Kohn-Sham level, using the PBE0 functional, which means that the values in the dissociation limit calls for some caution. We are, however, more interested in the behavior of the effective density around the equilibrium distance, which is 2.29 Å from experiment. From Table 1 we see that both the gradient and the Hessian of the effective density ρPb e is positive in this region, but we also see from Figure 2 that the effective density goes through a minimum around R= 1.96 Å. We have tried to rationalize this behavior using projection analysis. 92,93 We start by expressing the effective density as an SCF expectation value, which in the present case gives

0 ρ G n 0 = ∑ i ψ MO i ρ G n ψ MO i . (33)
We then expand the molecular orbitals ψ MO i in pre-calculated orbitals ψ A k of the constituent atoms

ψ MO i = ∑ kA ψ A k c A ki + ψ pol i . ( 34
)
where A refers to the atomic center and k is an orbital index for that center. The expansion is normally restricted to the orbitals that are occupied in the ground state configuration of each atom and is thereby not guaranteed to fully span each occupied molecular orbital. The orthogonal complement ψ pol i is denoted the polarization contribution and should normally be small for a meaningful analysis. Insertion of the expansion, eqn (34), into the expression for the expectation value, eqn (33), allows a decomposition of the latter into intra-and interatomic contributions. In the present case we find, as expected, that the effective density ρPb e is completely dominated by the intraatomic contribution associated with the atomic orbitals of lead. From previous analysis 94 we know that only the large components of s 1/2 -orbitals and the small components of p 1/2 -orbitals contribute to the contact density; for the effective density other atomic orbitals may contribute as well, but their contributions can be ignored. In the present case we see that the effective density ρPb e is dominated by the contribution from Pb s 1/2 -orbitals, but it is interesting to see that there is a non-negligible p 1/2 -contribution at intermediate internuclear distances. By looking at contributions from individual pairs of lead orbitals, we find that there are a number of large contributions of opposite sign which leads to the minimum of Figure 2, but we have so far not succeeded in formulating a simple mechanism explaining this observation.

Comparison with experimental data

Let us now compare our CCSD results to available experimental data. In Table 4 we report calculated field shift parameters for lead and thallium in the title compounds. The reported values have been obtained by using our CCSD values of

ρ[1]X e , ρ[2]X e
as well as the cubic force constant f 3 and combining them with equilibrium bond lengths r e and harmonic vibrational frequencies ω e taken from experiment. [START_REF] Huber | Constants of Diatomic Molecules (Data Prepared[END_REF] The latter values have been converted to equilibrium force constants k e using the atomic masses of the most abundant isotopes, corresponding to the reference isotopomer.

We first compare our rotational field shift parameters V 01 with the values reported by Tiemann and co-workers in 1982. 39,[START_REF] Tiemann | [END_REF] Although our calculated values are systematically outside experimental error bounds, with relative errors on the order of -20%, they are of the right order of magnitude, correctly reproduce trends observed in the experimental data and therefore seem to corroborate the observation of nuclear size effects in rotational spectra by Tiemann and co-workers. What is curious is that the field shift parameters V Pb 01 for PbS, although 40 On the other hand, their reported vibrational field shift parameter V Pb 10 for PbS is 6.84(41)•10 4 Å-2 , which is in line with our theoretical value 7.87•10 4 Å-2 . In the 1985 paper Knöckel et al. 40 Eq.( 7) report ρ [2]Pb 0 = 47(61)• Å-5 for PbS , which deviates significantly from our CCSD value of 131 Å-5 . However, if we use the rotational field shift parameter V Pb 01 for PbS reported by the Tiemann group in 1982, 39,[START_REF] Tiemann | [END_REF] the vibrational field shift parameter V Pb 10 from 1985 40 combined with our CCSD cubic force constant we obtain the value ρ [2]Pb 0 = 137 Å-5 which is in very good agreement with our result. Serafin et al. completed the series of Pb-chalcogenides by studying the nuclear volume effect in the rotational spectrum of PbO. 43 They report a rotational field shift parameter V Pb 01 which is in reasonable agreement with our CCSD value. On the other hand, they give the 1985 formula for V 01 , eqn (28), and so their value for ρPb [1] is about an order of magnitude smaller than our CCSD value. Interestingly, they cite the 1985 value of Knöckel et al. 40 Eq.( 7) V Pb 01 for PbS, but a factor of ten smaller, so in line with our CCSD results.

In Table 5 we have compiled the change of ground electronic state potential energy, equilibrium distance and force constant, that is, δ E AA , δ R AA e and δ k AA e , upon isotope substitution of the heavy atom in the title compounds. As illustrated for lead in Figure 3, the nuclear radii of the isotopes of a given element are monotonically increasing as a function of atomic number, despite the well-known odd-even staggering. As seen from eqn (15), this implies that the ground electronic state potential energy is monotonically increasing with atomic mass. Interestingly, the changes of bond lengths δ R AA e due to isotope substitution shown in Table 5 are all positive, although nuclear radii decrease. From eqn (24) Table 2 Spectroscopic constants (bond lengths r e , harmonic vibrational frequencies ω e , and anharmonic constants ω e x e ) of diatomic molecules TlY (Y=F, Cl, Br, I) obtained with various method and using TZ+2s1p basis set, as well as the gradient and Hessian of the contact and effective (number) density. All experimental spectroscopic data (r e , ω e , ω e x e ) were taken from ref. 95.

Method

r e Å ω e cm -1 ω e x e cm -1 ρ one sees that the sign of δ R AA e is determined by the gradient of the effective density, which is positive for all title compounds in our study. What is quite amazing is that the NVE modification of spectroscopic constants such as the equilibrium rotational and force constant is due to bond length changes that are on the order of one percent of the nuclear radius of the heavy atom, yet apparently seen in high-resolution rovibrational spectroscopy. The change in the equilibrium force constant upon reduction of the heavy atom nuclear radius is negative for all title compounds. As shown in eqn (29), the change in the equilibrium force constant has contributions from both the gradient and the Hessian of the effective density. Since cubic force constants are typically negative, as for instance inferred from a Morse potential and as shown in Table 4, the two contributions have the same sign if the gradient and the Hessian of the effective density have the same sign. From our CCSD results reported in Tables 1 and2 this is seen to be the case of the lead-chalcogenides, but not for the thallium halides; for the first class of compounds the contribution from ρA [2] e is on the order of 10 %, whereas it is on the order of 30 -40 % for the latter.

[1]T l 0 -4 ρ[1]Tl e -4 ρ [2]T l 0 -5 ρ[2]Tl e -

Isotope fractionation

We close this section by a short discussion of the nuclear volume effect on chemical equilibria. In his seminal 1996 paper 25,[START_REF] Bigeleisen | Isotope Effects In Chemistry and Biology[END_REF] Bigeleisen considered the nuclear volume effect on the isotopic exchange reaction A Y + AX = A X + AY (35) and in particular the isotope fractionation factor

α = (A/A ) AY (A/A ) AX , (36)
Table 3 Statistical analysis of relative errors (the mean relative error and the standard deviation) in bond lengths r e , harmonic vibrational frequencies ω e , and anharmonic constants ω e x e obtained for the eight diatomic molecules PbX (X=O, S, Se, Te) and TlY (Y=F, Cl, Br, I) with different methods and compared to experimental data. [START_REF] Huber | Constants of Diatomic Molecules (Data Prepared[END_REF] We also give relative errors of calculated first-and second derivatives of effective densities using our CCSD results as reference which compares distributions of isotopes A and A in compounds AY and AX. He concluded that the field shift effect would modify the isotope fractionation factor α 0 , as predicted by Bigeleisen-Goeppert-Mayer theory, 23 by

Method

ln α = ln α 0 + ln K f s , (37)
where appears the correction factor

ln K f s = - 1 kT δ E AA (AX) -δ E AA (AY ) , (38)
due to the change of the electronic ground state energies induced by the nuclear volume effect. Using eqn 15 we can re-write the above expression as

ln K f s = - 1 kT Z A e 2 6ε 0 δ r 2 n AA ρA e (AX) -ρA e (AY ) , (39)
which shows that the nuclear volume effect favors the compound having the smallest effective density at center A. Normally this would imply the compound with the smallest s 1/2 -contribution to the contact density, but for heavy elements one can not neglect the possible importance of p 1/2 -contributions, as seen from Figure 2.

A key result of Bigeleisen-Goeppert-Mayer theory is that the isotope fractionation factor α 0 can be expressed in terms of reduced partition-function ratios

α 0 = s s f (AY /A Y ) s s f (AX/A X) , (40)
where the reduced partition-function ratio in turn is defined as

s s f = s s q q N ∏ j m j m j , -3/2 (41)
where appears partition functions q and q , symmetry numbers s and s from the rotational part of the partition functions as well as the atomic mass ratios. Bigeleisen and Goeppert-Mayer in particular noted that the latter factor would cancel out in the calculation of the equilibrium constant for any isotopic exchange reaction, eqn (35). The natural logarithm of the reduced partition-function ratio is conveniently expressed as

ln s s f = 3N-6(5) ∑ i G (u i ) u i -u i ; u i = hω i kT ; G(u) = 1 2 - 1 u + 1 e u -1 ;
(42) this expression is valid when the difference in harmonic frequencies ω i between the two isotopomers is small. In a 1998 paper Bigeleisen considered the effect of the nuclear volume effect on this frequency difference for a diatomic molecule. 97 Bigeleisen based his estimation on a plot by Tiemann et al. [START_REF] Tiemann | [END_REF] similar to that of Figure 1, but based on experimental data. Tiemann and coworkers proposed based from the linearity of the plot that i) the Hessian of the contact (effective) density would be almost constant for this series of molecules and that ii) the slope of the linear fit therefore suggested that the relative field shift effect would be on the order of 10 -6 for vibrational energies. Our results indicate that the first suggestion is wrong, yet the estimates of Bigeleisen are to good measure correct, as we show in the following. We may write the harmonic frequency of the substituted isotopomer as

ω e = k e µ = k e + δ k AA e µ ≈ k e µ + 1 2 δ k AA e µ k e = ω [0] e + 1 2 δ k AA e k e ω [0]
e .

(43) The first term is the harmonic frequency ω

[0]
e in the absence of the nuclear volume effect. From the second term it follows that the nuclear volume effect may be expressed relative to this zerothorder harmonic frequency and that the weight is determined by the ratio between the change of force constant δ k AA e due to NVE and the force constant k e of the reference isotopomer. This ratio is given for all title compounds in Table 5 and is indeed of order 10 -6 . Using the experimental spectroscopic data from Table 1 we find that upon isotope substitution 208 Pb → 207 Pb in 208 Pb 32 S, the change of harmonic frequency due to the mass effect alone is

ω e -ω [0]
e is 0.1396 cm -1 , whereas the nuclear volume effect gives an additional shift of 2.514•10 -4 cm -1 . These numbers agree quite well with those of Bigeleisen and corroborate his conclusion that the nuclear volume effect on the force constant has very little influence on the equilibrium constant in isotopic exchange reactions. h Table 4 Rotational V 01 and vibrational V 10 field shift parameters of the heavy atom in the selected diatomic molecules. Theoretical values have been obtained using CCSD geometric derivatives of effective densities ( ρ [1] e and ρ [2] e ) and cubic force constants f 3 . For r e and ω e we have used experimental data, [START_REF] Huber | Constants of Diatomic Molecules (Data Prepared[END_REF] except for TlI, where we used the CCSD(T) values reported by Knecht and Saue. 44 Experimental values of V 01 are those given by Tiemann and co-workers in 1982, 39,[START_REF] Tiemann | [END_REF]

Conclusion

In the present contribution we have investigated the nuclear volume effect in the rovibrational spectra of diatomic molecules, motivated by the pioneering work of Tiemann and coworkers. 37,39,40,[START_REF] Tiemann | [END_REF] We provide a new derivation of the electrostatic interaction between electrons and nuclei which is responsible for the NVE, as well as the isomer shift of Mössbauer spectroscopy. This derivation has the advantage of not assuming any specific model for the nuclear charge distribution and naturally brings out the nuclear mean square radius r 2 n as a measure of nuclear size. It also leads to the introduction of the effective density ρA e which corresponds to the value of the electronic number density at some radius inside the nuclear volume. The effective density is often approximated by the contact density ρ A 0 , but the present and previous studies 44,62,94 show that for heavy atoms this leads to errors, albeit systematic, on the order of 10%. We furthermore argue that this is an unnecessary approximation, since the effective and contact densities can be calculated at identical computational cost.

Based on the energy expression, eqn (13), we derive, using variational perturbation theory, expressions for the change of energy, rotational constant B e and force constant k e due to the NVE in terms of the field shift parameters V A 00 , V A 01 and V A 10 , respectively. These expressions agree with the 1982 derivations of Tiemann and co-workers, 39,[START_REF] Tiemann | [END_REF] but we find no justification of the later introduction of a scaling factor π 2 in the rotational field shift parameter V A 01 . 40 The field shift parameters are expressed in terms of the effective density ρA e and its first and second derivative with respect to internuclear distance. We have calculated these quantities for lead chalcogenides and thallium halides at the level of 4component relativistic CCSD theory and conclude that our results corroborate the experimental findings of the Tiemann group, although we observe relative errors on the order of 20 %. This is all the more remarkable in that the modifications of spectroscopic constants primarily arise from bond lengths changes on the order of one percent of nuclear radii. We have also carried out 4-component relativistic DFT calculations, but they do not appear to be reliable for the Hessian of the effective density, at least for the lead chalcogenides.

For lead sulfide we investigate the effective density at the lead nucleus as a function of internuclear distance and find that it has a minimum around 0.2 Å inside the equilibrium bond length. Projection analysis shows that the effective density is dominated by the intra-atomic contribution of the atomic orbitals of lead, with contributions from both s 1/2 -and p 1/2 -orbitals, the latter significant at intermediate internuclear distances. We have on the other hand not been able to explain the appearance of the minimum of the effective density at short internuclear distance.

We also discuss modifications of Bigeleisen-Goeppert-Mayer theory of isotope fractionation in light of our results. We conclude, as Bigeleisen, 25,97 that the modification of electronic ground state energies due to NVE has an important bearing on the isotope chemistry of the heavy elements, whereas the modifications of force constants due to NVE can be ignored at this level of approximation.

In future work we would like to further check the convergence of our results in terms of basis sets and correlation treatment. We would in particular like to investigate the effect of going beyond the 4-component relativistic CCSD level of theory in the study of nuclear volume effects, notably the effect of the inclusion of triples. We would also like to investigate possible effects of nuclear spin on heavy element isotope chemistry.

r e =r A r e

 e =0 ρ e (r e ) φ A (r e )φ point A (r e ) d 3 r e e ρ e (r e )φ point A

  E A en;1 = -e ρA e r e =r A r e =0 φ A (r e )φ point A (r e ) d 3 r e ; ρA e = ρ e (r); r ∈ 0, r A (11) where the effective density ρA e corresponds to the value of the electron (number) density at some radius inside the nuclear volume of atom A. Insertion of the Legendre expansion now leads to ) d 3 r e .

Fig. 1

 1 Fig. 1 Gradient of effective densityρ[1] e (in Å-4 ) plotted as function of equilibrium bond length R e (in Å). Linear regression gives y = 802.12 -236.91x

Fig. 2

 2 Fig. 2 Projection analysis of the effective density ρPb e (in Å-3 and relative to values at R = 4.0 Å) of the PbS molecule, using the PBE0 functional, as a function of internuclear distance R (in Å)

Fig. 3

 3 Fig.3Root mean square (rms) nuclear charge radii (in fm) for lead as a function of atomic mass number A. Data is taken from the tables of Angeli and Marinova56

Table 1

 1 Spectroscopic constants (bond lengths r e , harmonic vibrational frequencies ω e , and anharmonic constants ω e x e ) of diatomic molecules PbX (X=O, S, Se, Te) obtained with various method and using TZ+2s1p basis set, as well as the gradient and Hessian of the contact and effective (number) density. All experimental spectroscopic data (r e , ω e , ω e x e ) were taken from ref.95 Ref. 43 (believed to be too small by a factor ten, see text) b Ref. 39 c Ref. 40 (believed to be in error, see text) a fitting parameter, was changed by about an order of magnitude from 2.45(19)•10 4 Å-2 to 26.38(51)•10 4 Å-2 in the 1985 experimental study reported by the same group.

	Method	r e Å	ω e cm -1	ω e x e cm -1	ρ	[1]Pb 0	Å-4	ρ[1]Pb e	Å-4	ρ 0 [2]Pb	Å-5	ρ[2]Pb e	Å-5
	PbO
	HF	1.8717	825.1	5.29		473	423	817	733
	CCSD	1.9060	754.1	4.74		394	353	233	210
	PBE0	1.9017	774.9	5.19		384	343	21		20
	PBE	1.9346	716.7	4.92		362	324	-112	-100
	BLYP	1.9524	691.1	3.69		372	333	-148	-131
	B3LYP	1.9211	744.8	5.38		390	349	-40	-35
	CAMB3LYP	1.8982	788.2	17.93		397	355	314	221
	Exp.	1.9218	712	3.54		-		40.4(31) a	-		-
	PbS
	HF	2.2672	464.6	1.31		310	278	308	276
	CCSD	2.2844	440.3	1.25		268	239	145	131
	PBE0	2.2804	440.4	0.53		262	235	63		57
	PBE	2.3121	412.7	0.39		237	212	4		4
	BLYP	2.3398	396.3	0.99		237	212	-15	-13
	B3LYP	2.3075	419.1	0.34		257	230	31		29
	CAMB3LYP	2.2754	444.2	0.28		280	250	345	305
	Exp.	2.2868	429.4	1.3		-		213(17) b	-		47(61) c
	PbSe
	HF	2.3906	294.2	0.66		281	251	204	183
	CCSD	2.4062	282.6	0.58		236	211	73		66
	PBE0	2.4007	287.1	0.80		227	203	4		6
	PBE	2.4329	272.7	1.47		201	180	-39	-32
	BLYP	2.4616	258.8	1.50		197	176	-54	-46
	B3LYP	2.4290	274.6	1.25		218	195	-18	-15
	CAMB3LYP	2.3952	287.1	-1.76		244	218	-134	-67
	Exp.	2.4022	277.6	0.51		-		177(15) b	-		-
	PbTe
	HF	2.5995	215.3	0.54		240	215	134	121
	CCSD	2.5947	216.0	0.40		191	171	56		51
	PBE0	2.5971	215.1	0.27		181	162	4		5
	PBE	2.6313	201.9	-0.04		155	139	-19	-16
	BLYP	2.6656	191.8	0.01		148	132	-28	-24
	B3LYP	2.6300	202.8	0.07		170	152	-8		-7
	CAMB3LYP	2.5927	214.8	1.63		197	176	-6		-3
	Exp.	2.5949	211.9	0.43		-		148(11) b	-		-

a

  5

	TlF
	HF	2.0801	478.1	2.38	438	393	-411	-368
	CCSD	2.0887	473.3	2.18	368	330	-461	-414
	PBE	2.1190	440.3	1.44	307	276	-521	-468
	PBE0	2.0905	471.8	5.23	343	308	-513	-461
	BLYP	2.1415	421.5	-2.94	308	276	-530	-474
	B3LYP	2.1124	451.3	2.53	335	300	-523	-469
	CAMB3LYP	2.0912	473.2	-9.09	351	315	-375	-366
	Exp.	2.0844	477.3	2.3	-	254(30) a	-	-
	TlCl
	HF	2.5264	268.9	0.88	310	278	-256	-228
	CCSD	2.5064	277.2	0.84	247	222	-266	-238
	PBE	2.5198	264.2	-1.68	190	171	-249	-224
	PBE0	2.4955	273.1	-1.64	222	199	-264	-237
	BLYP	2.5577	259.4	1.06	186	166	-259	-231
	B3LYP	2.5306	265.7	-0.29	210	189	-268	-241
	CAMB3LYP	2.5005	273.1	-11.10	238	213	-306	-307
	Exp.	2.4848	283.7	0.81	-	185(25) a	-	-
	TlBr
	HF	2.6720	180.7	0.40	264	237	-203	-182
	CCSD	2.6393	189.0	0.40	206	184	-207	-186
	PBE	2.6546	190.4	2.42	155	139	-186	-166
	PBE0	2.6325	195.4	1.73	184	165	-200	-180
	BLYP	2.6931	176.4	1.90	149	133	-195	-173
	B3LYP	2.6667	186.8	2.13	171	154	-207	-184
	CAMB3LYP	2.6353	200.2	17.63	200	179	-333	-267
	Exp.	2.6181	192.1	0.39	-	142(42) a	-	-
	TlI
	HF	2.8867	138.8	0.25	212	190	-130	-117
	CCSD	2.8389	147.0	0.25	160	143	-127	-113
	PBE	2.8565	137.7	-0.73	114	102	-104	-92
	PBE0	2.8315	144.2	-0.31	140	126	-119	-107
	BLYP	2.9041	133.1	-0.18	107	95	-109	-97
	B3LYP	2.8731	137.1	-0.52	127	114	-122	-108
	CAMB3LYP	2.8335	145.4	4.70	155	138	-23	-58
	Exp.	2.8136	150	-	-	120.6(38) a	-	-

a Ref.

39

  except for PbO, taken from ref.43.

	Molecules	V 01 10 4 Å-2	f 3 E h /a -3 0	V 10 10 4 Å-2
		theory	exp.	theory	theory	exp.
	PbO	3.28	2.64 (36)	-0.957	10.94
	PbS	2.75	2.45 (19)	-0.456	7.87 6.84(41)
	PbSe	2.66	2.21 (19)	-0.384	7.37
	PbTe	2.47	2.12 (16)	-0.299	7.15
	TlF	5.31	4.09 (19)	-0.385	9.95
	TlC	4.94	4.09 (55)	-0.179	8.13
	TlBr	4.44	3.37 (10)	-0.151	7.47
	TlI	3.96	3.20 (10)	-0.117	7.33

| 1-14

[1][2][START_REF]Isotope Effects In Chemistry and Biology[END_REF][START_REF] Wolfsberg | Isotope Effects: in the Chemical[END_REF][START_REF] Eiler | [END_REF][6][7][START_REF] King | Isotope Shifts in Atomic Spectra[END_REF][START_REF] Kragh | [END_REF][10][11][12][13][START_REF] Bodmer | Proceedings of the Physical Society. Section A[END_REF]

1-14 J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Acknowledgments

Computing time from CALMIP (Calcul en Midi-Pyrenées) is gratefully acknowledged. AA was financially supported by a PhD-grant from the Syrian Ministry of Education as well as kind support

h Table 5 Change of ground electronic state potential energy, equilibrium distance and force constant upon isotope substitution of the heavy atom in the title compounds. Nuclear radii changes δ AA r 2 n are taken from the tables of Angeli and Marinova 56 and are (in fm