# Deviations from Born-Oppenheimer mass scaling

in spectroscopy and ultracold molecular physics

###### Abstract

We investigate Born-Oppenheimer breakdown (BOB) effects (beyond the usual mass scaling) for the electronic ground states of a series of homonuclear and heteronuclear alkali-metal diatoms, together with the Sr and Yb diatomics. Several widely available electronic structure software packages are used to calculate the leading contributions to the total isotope shift for commonly occurring isotopologs of each species. Computed quantities include diagonal Born-Oppenheimer corrections (mass shifts) and isotopic field shifts. Mass shifts dominate for light nuclei up to and including K, but field shifts contribute significantly for Rb and Sr and are dominant for Yb. We compare the ab initio mass-shift functions for Li, LiK and LiRb with spectroscopically derived ground-state BOB functions from the literature. We find good agreement in the values of the functions for LiK and LiRb at their equilibrium geometries, but significant disagreement with the shapes of the functions for all 3 systems. The differences may be due to contributions of nonadiabatic terms to the empirical BOB functions. We present a semiclassical model for the effect of BOB corrections on the binding energies of near-threshold states and the positions of zero-energy Feshbach resonances.

###### keywords:

Born-Oppenheimer approximation, adiabatic correction, ultracold molecules, contact density, isotopic field shift, Born-Oppenheimer breakdown function, fifth force###### Pacs:

31.15.ae, 31.15.vn, 83.10.-y## 1 Introduction

The Born-Oppenheimer approximation (BOA) lies at the heart of chemical and molecular physics. It underpins the concepts of potential energy curves and surfaces that are universally used to understand and interpret molecular structure and dynamics. For many purposes, the BOA is adequate and it is not necessary to go beyond it. However, for light nuclei and high-precision work, deviations from the BOA are important. Quantitative investigations of such deviations go back at least to the theoretical work of Kołos and Wolniewicz on H Kolos:1964 (); Kolos:1965 (), which stimulated reinterpretation of the experimental spectrum by Herzberg Herzberg:1970 (). They are also important for H Polyansky:1999 (), and they have been characterized spectroscopically for hydrides such as HeH Coxon:HeH+:1999 (), BeH Coxon:BeH+:1997 (), HF LeRoy:1999 (); Coxon:HF:2006 (); Coxon:HX:2015 (), HCl Coxon:HCl:2000 (); Coxon:HX:2015 (), HBr and HI Coxon:HX:2015 (), AgH LeRoy:2002 (); LeRoy:2005 (), LiH Coxon:LiH:2004 (), BeH LeRoy:2006 () and MgH Henderson:2013 () and for CO Coxon:CO:2004 (), Li WMLjcp2002 (); Coxon:2006 (); LeRoy:2009 (), LiK Tiemann:2009 (), and LiRb Ivanova:2011 (). For molecules without such light nuclei, the deviations have been hard to detect Seto:2000 (); Docenko:2007 (); Strauss:2010 (), although indications of them have been seen in K Falke:2008 (), Rb Seto:2000 (); VKKpra2009 (), and I Knockel:2004 (); Salumbides:2008 ().

Recent developments in the field of ultracold atoms and molecules offer a new stimulus to understand deviations from the BOA. Key quantities in this field are the binding energies of levels very close to dissociation and the positions of zero-energy Feshbach resonances as a function of magnetic field Chin:RMP:2010 (). The latter are essentially the fields at which the energies of bound molecular states exactly equal those of free atoms. For pairs of heavy atoms, potential curves derived from one isotopolog have been very successfully used to predict resonance positions for another by simply rerunning the scattering calculations with a different reduced mass (and different atomic properties such as nuclear spins and hyperfine splittings) Simoni:2008 (); Cho:RbCs:2013 (); Blackley:85Rb:2013 (). However, Julienne and Hutson Julienne:Li67:2014 () have recently shown that deviations from the BOA are responsible for 4 G of the shift in resonance position between Li+Li and Li+Li, and have obtained potential curves that include the necessary corrections for both the singlet and triplet states.

Breakdown of the BOA is also crucial for attempts to use the spectroscopy of ultracold molecules to explore fundamental physics. For example, Kitagawa et al. Kitagawa:2008 () have measured the binding energies of near-dissociation states of several isotopologs of Yb and similar experiments are underway for Sr Zelevinsky:2006 (); Stellmer:2012 (); McGuyer:2013 (); Borkowski:2014 (); McGuyer:2014 (). For Yb, binding energies are generally in good agreement with the predictions of Born-Oppenheimer mass scaling. However, there are proposals to use the small deviations from such scaling to place limits on the magnitude of a “fifth force” that may exist in addition to the familiar electromagnetic, gravitational and strong and weak nuclear forces (see, e.g., Refs. Salumbides:2013 (); Niu:2014 (); Salumbides:2014 ()). Before any such effects can be ascribed to novel physics, it is crucial first to consider deviations from mass scaling that arise from conventional physics.

The purpose of the present paper is to explore the capabilities of current theoretical methods for calculating deviations from the BOA, and to investigate their magnitude for species of importance in the study of ultracold molecules. These include both homonuclear and heteronuclear alkali-metal diatomics and molecules such as Sr and Yb. For the heavier molecules, almost nothing is known about the corrections needed, and even order-of-magnitude estimates are valuable. We consider effects due to both finite nuclear mass (isotopic mass shifts) and finite nuclear volume (isotopic field shifts).

## 2 Theory

Electronic structure theory provides a framework for computing isotope shifts. Atomic and molecular calculations are usually performed assuming a priori that the nuclei are both infinitely heavy and infinitely tiny, i.e., they are treated as point charges. The former approximation is commonly called the BOA although it resembles the treatment of Born and Huang BHbook1956 () more than it does the original proposition of Born and Oppenheimer BOap1927 ().

For atomic systems, the theory of isotopic shifts has been rigorously developed within a relativistic formalism GNRjpb2011 () and dedicated programs are widely available for their ab initio calculation (see, e.g., Refs. NGGGJcpc2013 (); JGBcpc2013 ()). Meanwhile in molecular systems, bonding-induced isotopic shifts produce relatively small effects in spectroscopic results. To help diagnose whether isotopic mass shifts are important in a given application, Born and Huang derived a first-order perturbative correction to the BOA energy BHbook1956 (), which can take the form,

(1) |

where indicates the isotopolog, is the kinetic energy operator for nucleus and is the normalized electronic wavefunction for state obtained within the BOA, with explicit dependence on electronic positions and parametric dependence on nuclear positions . The quantities are known as adiabatic corrections or diagonal Born-Oppenheimer corrections (DBOCs).

The approach embodied by Eq. 1 appears to be problematic from a formal perspective, because the equations involve manipulating continuum functions as if they were normalizable. Fortunately, as was later shown by Kutzelnigg Kutzelnigg1997 (), this form of the adiabatic correction is correct and, in fact, can be derived rigorously by avoiding the actual specification of relative coordinates in the center-of-mass separation.

The equations of the pragmatic ansatz of Eq. 1 were first solved using Hartree-Fock (HF) wavefunctions by Sellers and Pulay SPcpl1984 (), and later by Handy et al. Handy1986a (). The expressions are evaluated by holding all but one nucleus fixed and calculating analytic derivatives of the wavefunction with respect to the coordinates of the remaining nuclei. The significance of the correction was established by early investigations of its effect on molecular bond lengths and vibrational frequencies Sellers (); Handy1996a (), thermochemical reaction barrier heights Truhlar (); Handy1996b () and the singlet-triplet gap in methylene Handy1986b (). Adiabatic corrections for water were a critical component of models that demonstrated the existence of water on the sun BWjcp1978 (); water1 (); water2 (); water3 (). Computing adiabatic corrections at the HF level has by now become so routine that they are included in standard composite methods for high-accuracy thermochemical calculations HEAT1 (); HEAT2 (); W4 (); W2F12 ().

Valeev and Sherrill showed that the inclusion of electron correlation via configuration interaction in the wavefunction can lead to changes in the absolute DBOC of a few percent for some systems, with the most pronounced effects occurring in hydrides VSjcp2002 (); TVSjpca2004 (). When changes with respect to geometry are considered instead, correlation effects can contribute much more significantly GLRcpl2001 (), in the same way as correlation can contribute more significantly to relative than to total energies. The CFour package CFour () has made available codes for the analytic evaluation of adiabatic corrections using wavefunctions from coupled-cluster (CC) calculations GTKjcp2006 () and Möller-Plesset perturbation theory MPDBOC (). Schwenke also evaluated adiabatic and non-adiabatic corrections Schwenke1 () using internally contracted multireference CI wavefunctions Schwenke2 (); Schwenke3 (), but unfortunately his program was not widely distributed. It is reported IASHjpca2016 () that the next release of the GAMESS software package will have the capability to compute DBOCs based on scalar-relativistic Hamiltonians.

Nonadiabatic corrections, originating from the off-diagonal matrix elements of the nuclear kinetic energy operator, can be essential for understanding molecular dynamics when different electronic states come close together Tbook1976 (). Even for nondegenerate electronic states, they can make significant contributions to spectroscopic line positions, as discussed below. Nevertheless, an adiabatic representation is convenient because it retains the concept of a potential energy surface IUPACCCT (); Bbbo2006 () and can be obtained using standard analytic derivative techniques Handy1986a (); GTKjcp2006 (). For these reasons, we focus here on obtaining adiabatic corrections, though we note that good progress has recently made for obtaining highly accurate nonadiabatic corrections by Pachucki and Komasa PKjcp2008 (); PKjcp2009 ().

There is a fundamental difference between the interpretations of the Born-Oppenheimer approximation in common use among spectroscopists and electronic structure theorists. These are most easily illustrated by considering the two different treatments of a diatomic molecule.

Electronic structure theorists normally consider nuclei moving on the potential energy curves or surfaces. For a diatomic molecule, the internuclear distance is and the reduced mass for nuclear motion is , where and are nuclear masses. This separation gives the form of the Born-Oppenheimer approximation described above. However, it has the disadvantage that the nonadiabatic corrections are nonzero even as , and this presents problems for scattering theory because the asymptotic wavefunctions are not simple products of the wave functions of the separated atoms.

In spectroscopy and scattering theory, by contrast, it is common to consider atoms moving on effective potential energy curves. Electrons are considered to be parts of the atoms and to move with them. This makes good physical sense, at least for core electrons, which are tightly bound to the nuclei. For a diatomic molecule, the reduced mass for atomic motion is , where and are atomic masses.

The formal justification of this approach is based on work by Bunker and Moss Bunker:effective:1977 () and Watson Watson:1980 (). Bunker and Moss derived an effective Hamiltonian for a single electronic state of diatomic molecule, taking nonadiabatic couplings into account by means of a contact transformation. Their treatment introduces nonadiabatic corrections through both an effective potential term (which scales as ) and separate -dependent reduced masses and for the vibrational and rotational motion, which differ from by terms that scale as . In subsequent work, Bunker, McLarnon and Moss Bunker:H2:1977 () showed that, if the -dependence is neglected, the value of that gives an optimum fit to the full nonadiabatic energies of H is closer to than to . Watson Watson:1980 () showed that the effective Hamiltonian of Ref. Bunker:effective:1977 () may be rearranged to a form containing the “charge-modified reduced mass” , where is the molecular charge; reduces to for a neutral molecule. The replacement of with substantially reduces the remaining nonadiabatic corrections if the electrons move with their respective nuclei; that is, and are much closer to than to . It essentially corresponds to a change of viewpoint: if is used as the reduced mass, then the nonadiabatic terms account for the extent to which the electrons fail to move with their “parent” atom, whereas if is used then the nonadiabatic terms need to account for the extent to which the electrons move with the nuclei at all.

Watson Watson:1980 () showed that it is not possible to determine both the adiabatic correction and the nonadiabatic correction terms in simultaneously from transition frequencies alone, and gave an expression for an effective adiabatic correction that absorbs the nonadiabatic corrections in (but not those in ).

Watson’s nonadiabatic corrections become asymptotically zero when the electrons of each atom move with it and the atomic reduced mass is used. However, his effective adiabatic corrections are not necessarily asymptotically zero. Nevertheless, we are free to choose the zero of energy for each isotopolog, and it is convenient to choose it as the energy of the free atoms at the threshold of greatest interest (which is the ground-state dissociation energy in the present work). With this choice, the value of at infinity for each isotopolog is absorbed into the definition of the origin, and only is explicitly included as an adiabatic correction. Pachucki and Komasa PKjcp2008 (); PKjcp2009 () give a perturbative treatment of nonadiabatic effects that reaches the same conclusion. Non-zero adiabatic corrections at infinity are still required for any electronic states that dissociate to different limits.

In the context of electronic structure theory, Handy and Lee Handy1996a () recommended that atomic masses be used rather than nuclear masses when computing adiabatic corrections. Kutzelnigg Kutzelnigg2007 () also discussed this question in detail, and concluded that at least inner-shell electrons should be considered to move with the nuclei. Nevertheless, the usual convention in electronic structure theory is to use nuclear masses GLRcpl2001 (), and we follow that convention in the remainder of this paper.

For atomic systems, the quantity that corresponds to the molecular adiabatic correction is the nuclear mass shift. Traditionally, the total atomic mass shift was separated into a normal mass shift and a specific mass shift Hughes:1930 (). The normal mass shift is obtained simply by replacing all electron masses with reduced masses in the calculation, and results in a simple scaling of all atomic state energies (and transition frequencies) by a factor . The specific mass shift, however, varies from state to state and its calculation involves mass polarization terms written in terms of products of momentum operators on pairs of electrons. More recently, it has been recognized that this approach can be unreliable, and that a more complete result may be obtained in terms of the relativistic nuclear recoil operator Fricke:1973 (). At least in the nonrelativistic case, the adiabatic correction as evaluated by electronic structure programs such as CFour is asymptotically equivalent to the sum of the total mass shifts of the constituent atoms.

For heavy atoms, isotope shifts are usually dominated by the nuclear field shift, which results from the finite volume of the nucleus, rather than by the mass shift. In quantum mechanics, electrons can penetrate nuclei, and the electric potential they experience inside the nucleus is less negative than . Consequently, the energies of penetrating orbitals are shifted upwards due to the finite size of the nucleus. The magnitude of this effect depends on the structure of the nucleus involved, but can be shown to a rough approximation OttenBook () to scale with the atomic number and mass number as

(2) |

In atomic spectroscopy the interplay between mass and field shifts has been studied extensively and a crossover point is estimated to occur at BDNps2013 ().

The theory of isotopic field shifts was first formulated for atomic systems by Rosenthal and Breit RBpr1932 () and by Racah Rn1932 () in 1932. Assuming a spherical nuclear charge distribution and performing a power expansion of the electron density within the nucleus, the conventional model gives the first-order perturbative correction to the energy shift of level in going from isotope to for an atom as

(3) |

Here is the density at the nucleus, often known as the contact density, and to a first approximation the so-called “nuclear parameter”, , is

(4) |

where is the difference in nuclear rms charge radii between isotopes and , sometimes indicated by the notation . A tabulation of nuclear mean-square charge radii is available Angeli:1999 (), but the results are not to very high precision.

The quantity usually shows an odd-even staggering (a “saw-tooth” pattern), with exceptionally small values occurring for compact nuclei with neutron magic numbers 20, 28, 50, 82, and 126 Aadndt1987 (). It is for this reason that magic-number isotopes such as K, Rb, and Sr are often chosen as reference isotopes when tabulating values of AMadndt2013 (). However, we do not uniformly follow this convention here. In order to maintain consistency with work performed by other authors, we deviate from it for Li and Rb, where Li and Rb are chosen as reference isotopes.

Tiemann et al. Tiemann:1982 () investigated non-Born-Oppenheimer effects in the rotational spectra of group III/VII and IV/VI diatomic molecules. They found good agreement with Watson’s expressions, but with anomalously large corrections for heavy atoms (Tl, Pb). Schlembach and Tiemann Schlembach:1982 () subsequently showed that the anomalous values can be attributed to nuclear field shifts, and that these are the dominant isotope-dependent effect for these species. For the vibronic spectra of PbS, Knöckel and Tiemann Knoeckel:1982 () found that the field shift by itself could explain the isotope dependence and terms due to mass shifts were negligible. In a subsequent reevaluation of the experimental results, Knöckel, Kröckertskothen and Tiemann KKTcp1985 () revised the magnitude of the field shift downwards substantially, but retained the overall conclusion that field shifts dominate adiabatic corrections in the rotational spectra of heavy-atom systems.

In the molecular case, the field shift for the free atoms can again be absorbed into the zero of energy for each isotopolog, but its -dependence contributes to molecular binding energies and level spacings. For rotational spectra, the key quantity is the derivative of the contact density with respect to , evaluated near the equilibrium geometry. Cooke et al. have determined Dunham-type parameters from the rotational spectra of a wide variety of diatomic molecules containing heavy elements, and identified large Born-Oppenheimer breakdown effects that they interpreted as isotopic field shifts CGpccp2004 (); CMGjms2004 (); CGCcp2004 (); CKGpccp2005 (); EDCcp2007 (); SPDcpl2007 (); GBCpccp2008 (); DEGjms2008 (); KCRcjp2009 (); GBSjcp2011 (). For a subset of these studies they used density-functional theory with scalar relativistic corrections to calculate contact density derivatives, and found reasonable agreement with the experiments. However, Knecht and Saue KScp2012 () have carried out 4-component relativistic calculations on the TlI, PbTe, and PbS systems; they obtained substantial disagreement with both the results of Cooke et al. and the experiments, and questioned whether Cooke et al. had actually included relativistic corrections in Refs. CMGjms2004 (); CGCcp2004 (); EDCcp2007 ().

Expectation value and derivative approaches for obtaining approximate contact densities are currently under development within the vibrational and Mössbauer spectroscopy communities FRsb2012 (); Fccr2009 (); KFjctc2008 (); Fjcp2007 (); NHic2006 (); Nica2002 (); LHLjacs2002 (); LLLjacs2001 (); MTams1985 (); NPvDprb1978 (); dVTRjcp1975 (); Srmp1964 (). Within the context of such calculations it has been shown that use of a relativistic Hamiltonian is essential for obtaining accurate contact densities and their geometry dependence for heavy nuclei Fjcp2007 (); Fccr2009 (); KFvMtca2011 (); KScp2012 (); FZCjctc2012 (). Incorporating relativity directly into the electronic Hamiltonian can be done in various ways (see, e.g., Refs. Lmp2010 (); Scpc2011 (); Pcr2012 (); Ajcp2012 (); Lpr2014 () for recent reviews and perspective articles).

When working within a fully relativistic framework, the use of finite-size nuclear models is necessary to avoid singularities in the wavefunction. Finite-size nuclear models may seem at first problematic since the perturbation theory approach for obtaining contact densities outlined above assumes the point-nucleus approximation as a zeroth-order starting point. A more accurate method for obtaining field shifts involves integration of the electron density over the nuclear volume. It has been shown in Refs. KFvMtca2011 () and KScp2012 () that the error introduced by replacing such integration by the finite-nucleus contact density is on the order of 10% for absolute contact densities, with smaller errors for changes with respect to geometry.

Many additional aspects of the computational methods for contact densities have been considered in Refs. MLRcpl2008 (); MWRjcp2010 (); KFvMtca2011 (). The importance of electron correlation has been examined in Ref. KFvMtca2011 (), where it was shown for mercury fluoride systems that correlation effects at the CCSD level contribute as much as % to bonding-induced changes in contact densities. However, the resulting bonding-induced changes were no more accurate than those calculated at the HF level. Including a perturbative correction for triple excitations was found to affect the bonding-induced changes by %. Inclusion of core-valence correlation was also shown to be important, contributing at the same level as the perturbative triples correction. When instead density-functional theory (DFT) methods were considered, multiple studies concluded that DFT contact densities are of comparable accuracy to HF MWRjcp2010 (); KFjctc2008 (); KFvMtca2011 (), though hybrid functionals with HF exchange were shown to give better results than pure functionals KFjctc2008 (). It was also shown in Ref. MWRjcp2010 () that much more accurate results were produced by using basis sets in their fully uncontracted forms.

## 3 Results and discussion

### 3.1 Benchmarking computed isotopic mass shifts

There are not many molecular electronic structure packages that currently have the capability to calculate mass shifts (i.e. adiabatic or diagonal Born-Oppenheimer corrections) at the coupled-cluster level. In the present work we compute isotopic mass shifts for atomic and molecular species using the DBOC facility in the CFour package CFour (). Core orbitals are correlated in all calculations, since core-valence correlation has been shown to contribute significantly to DBOCs GTKjcp2006 (). Parallelized analytic derivative codes are used.

For all atoms except H, DBOC calculations were initially performed using small basis sets, in particular DZP for the alkali metals and Sr DZP1 (); DZP2 (); DZP3 (); DZP4 () and the WTBS basis set for Yb WTBS1 (); WTBS2 (). The larger ANO-RCC basis sets ANORCC1 (); ANORCC2 () were also employed in some cases. The ANO-RCC basis sets have previously been shown to provide excellent dissociation energies, geometrical parameters, and electric properties for systems involving heavy elements, as demonstrated, for example, in a recent study on LiCs SFOjpb2009 (). While these basis sets are generally recommended only for relativistic calculations, the derivative programs required for the evaluation of DBOCs within a relativistic framework are not available in the current public release of CFour. Employing basis sets designed for relativistic calculations in non-relativistic work is somewhat questionable. However, relativistic effects are significant only for heavier elements, where DBOCs are expected to become small compared to the field shift. Our overall conclusions regarding DBOCs are thus unlikely to be affected by the errors due to basis-set incompleteness and neglect of relativity.

CFour allows the use of both spin-restricted HF (RHF) and spin-unrestricted HF (UHF) references. Unrestricted methods seem at first sight an appealing choice, as they offer a better description of the highly stretched molecule near dissociation and are directly applicable to individual doublet atomic species at dissociation. However, they may find wavefunctions that have a lower symmetry than the nuclear framework, and such symmetry-broken solutions are also often spin-contaminated to some degree Stanton1994 (); Krylov2000 (); Schlegel2002 (); Schlegel2009 (). This may cause additional complications in the evaluation of second-order properties such as DBOCs. We explore such phenomena in detail in the following section. Where UHF methods were employed, the lowest-energy UHF eigenstate was located within a reduced computational symmetry (C) by following the appropriate eigenvalue of the orbital rotation Hessian matrix from the totally symmetric RHF solution to the symmetry-broken solution.

Before computing the molecular (bonding-induced) isotopic mass shifts of interest in this work, we first consider whether molecular calculations with CFour can yield accurate absolute values of atomic mass shifts. Here spin-restricted calculations were performed on the corresponding homonuclear diatomic system at large values of the internuclear distance. For the LiLi mass shift, DBOCs were computed at the RHF-CCSD/ANO-RCC and UHF-CCSD/ANO-RCC levels, resulting in values (per atom) of cm and cm, respectively. These compare favorably with the atomic physics literature value of cm, obtained by applying the appropriate Rydberg factors (see, e.g., King’s description in Ref. Kjms1997 ()) to the near-exact total energy of Li calculated by Puchalski and Pachucki PPpra2006 ().

Next we investigate the reliability of restricted and unrestricted references for the computation of molecular DBOCs. We first consider H, for which the nearly exact energies and adiabatic corrections by Kołos and coworkers are available for comparison Kolos:1964 (); Kolos:1965 (); KRjcp1993 (). Various approximate potential energy curves and absolute DBOC functions for HH are shown in Fig. 1. All potential energies are shown with respect to the exact asymptotic value, 0.5 Hartree, and the DBOC values are shown relative to an asymptotic value obtained from a molecular calculation at large . The behavior of the H potential curves as described by the RHF, UHF, RHF-CCSD, and UHF-CCSD methods is discussed in elementary textbooks and we include them in Fig. 1 only to contrast their characteristics with the corresponding DBOC functions.

The computed DBOC functions for HH are shown in Fig. 1. The UHF function is qualitatively wrong, exhibiting an unphysical pole-like feature near 2.2 bohr, which we discuss in detail below. The RHF function is smooth, but in poor quantitative agreement with the results of Kołos and Rychlewski KRjcp1993 (). However, the RHF-CCSD and UHF-CCSD methods are exact for this system, except for basis-set incompleteness, and exhibit errors in DBOCs no larger than 1%.

Figure 2 shows the sensitivity of the computed CCSD DBOC function to the quality of the basis set. Vertical lines mark the radial position of key points on the potential curve. These include the inner turning point at dissociation (so that ), the distance at the potential minimum, and the point where the energy is half way between the minimum and dissociation (). On the scale of the plot, the ANO-RCC and aug-cc-pVQZ basis set curves lie directly on top of the results of Kołos and Rychlewski KRjcp1993 (). Even the much smaller and more affordable DZP basis set gives results that are accurate within 10%. Because of this reasonable accuracy and its availability for most elements across the periodic table, the DZP basis set will be heavily utilized in this work for exploring trends in isotope shifts with atomic number.

To demonstrate the level of accuracy possible for DBOCs obtained with large basis sets, Fig. 2 shows the errors in the computed DBOC functions with respect to the results of Kołos and Rychlewski KRjcp1993 (). The ANO-RCC basis set gives errors that do not exceed 0.13 cm at distances larger than the inner turning point ( bohr). The ANO-RCC basis set is known to perform similarly to the aug-cc-pVQZ basis set, so we also computed DBOC functions with the aug-cc-pV5Z and aug-cc-pV6Z basis sets, where maximum errors were found to drop to 0.042 and 0.025 cm, respectively. The addition of midbond functions (designated in the figure legend as “+mb”) was also investigated and they were found to reduce errors greatly for bohr, while increasing errors somewhat for bohr. The performance of all basis sets tested here degrades rapidly in the region bohr, probably because modern basis sets are tuned for optimum performance near the equilibrium bond length ( bohr).

The unphysical pole-like feature in the UHF function in Fig. 1 is analogous to singularities that have been studied in the context of other properties including quadratic force constants CSASjcp1997 () and indirect nuclear spin-spin coupling constants AGcp2009 (). Such poles arise because the second derivatives of the correlated energies depend upon the orbital rotation parameters, which themselves are not continuously differentiable through the region of the transition from a symmetry-conserved to a symmetry-broken wavefunction. This phenomenon is sometimes referred to as an orbital instability envelope CSASjcp1997 ().

UHF-CCSD may also produce orbital instability envelope artifacts in molecules with more than 2 electrons. This is demonstrated in Fig. 3, which shows Li potential energy curves and DBOCs computed using RHF-CCSD and UHF-CCSD with the DZP and aug-cc-pVDZ basis sets. While the RHF-CCSD and UHF-CCSD potentials are virtually identical for each basis set, the corresponding DBOC functions are not. For both basis sets, the UHF-CCSD results exhibit an unphysical peak in a region where the DBOC radial function should asymptotically approach zero. For this reason we choose to use RHF-CCSD calculations in preference to UHF-CCSD calculations of DBOCs in the following sections.

### 3.2 Isotopic mass shifts for diatomic molecules topical in ultracold physics

Radial isotopic mass-shift functions, obtained from DBOCs computed at the RHF-CCSD/DZP level of theory, are shown in Fig. 4 for the ground states of all the alkali-metal dimers and the molecules Sr and Yb. These are all molecules of interest in the field of ultracold molecules. The mass shifts are defined as differences

(5) |

between the adiabatic corrections for isotopologs and . The radial positions of the corresponding potential minima are indicated by labeled and color-coded vertical lines FDVjcp2014 (); Kjcp2008 (); MPTZRijqc2009 (). The qualitative shape of the curves in Figs. 4 and 4 is inverted compared to those in Figs. 4 and 4, but this is simply because for Li, Sr and Yb the reference isotope is the heaviest, while for K and Rb it is the lightest.

We first consider LiLi mass shifts in alkali-metal diatomics containing Li, which are shown in Figure 4. The functions for Li and LiNa have the same general shape as for H, with a negative segment at long range and a positive segment at short range. The zero-crossing is close to for Li and LiNa, but moves outwards faster than for the heavier alkali metals, and for LiCs the function does not cross zero until about 14 bohr. This may well be due to the increasing contribution of charge transfer to the bonding in LiK, LiRb and LiCs, even at 10 to 15 bohr, as evidenced by their dipole-moment functions Aymar:2005 (). In Section 3.4 we will investigate further whether these mass-shift functions are in agreement with spectroscopically derived isotope shifts. It is worth noting that some mass-shift functions cross zero so close to their respective potential minima that reporting mass shifts at or at any other single point on the curve can be uninformative.

Isotopic mass shifts from DBOC calculations at the RHF-CCSD/DZP level are shown for isotopologs of alkali-metal diatomics involving K (with K as reference) in Fig. 4 and for those involving Rb (with Rb as reference) in Fig. 4. As expected, the mass shifts for K are about twice those for K, with the exact ratio determined by the changes in the isotopic mass. KCs curves are omitted from Fig. 4 because many of the corresponding RHF calculations did not converge. The overall magnitude of the mass shift for RbRbRbRb substitution, MHz, is consistent with the order-of-magnitude absolute-value estimate suggested in Ref. VKKpra2009 () for the correction due to Born-Oppenheimer breakdown.

Figure 4 shows mass shifts from DBOC calculations at the RHF-CCSD/DZP and RHF-CCSD/WTBS levels for homonuclear Sr and Yb diatomics, respectively, with Sr and Yb as reference isotopes. Both Sr and Yb have mass shifts that are positive at short range (but well outside the inner turning points, which are both between 7 and 7.5 bohr Stein:2008 (); MPTZRijqc2009 ()). The mass shifts for Sr have a significant negative component at long range, which dominates near . Yb, by contrast, has very weak mass shifts at long range, and the absolute values are at least an order of magnitude smaller at the equilibrium bond length than those for Sr or any of the other systems shown in Figs. 4 to 4. However, it should be noted that the minimal WTBS basis set used for Yb is considerably less flexible than the DZP basis set used for Sr, and is expected to produce too soft a short-range repulsive interaction and provide a poorer description of long-range forces Szabo (). These features are reflected in Fig. 4. Unfortunately, DBOC calculations with the larger basis sets that are currently available for Yb would be computationally very expensive. For now, computations at this level of theory must suffice; they demonstrate that DBOCs of Sr and Yb have similar qualitative features, with those for Yb being much smaller in magnitude.

### 3.3 Isotopic field shifts for diatomic molecules topical in ultracold physics

Field shifts can also contribute significantly to total isotope shifts. The field shift is approximately proportional to the contact density (Eq. 3), so this is the key quantity to compute using electronic structure theory. We have carried out extensive benchmark calculations on LiRb to compare the results of different approaches using the DIRAC, MOLCAS and ADF packages, which are described in the Supplementary Material. The packages all use different treatments of relativity, and offer different options for including electron correlation, either at the coupled-cluster level (DIRAC only) or using density-functional theory (all three programs). Unfortunately the different treatments offer less consistent results than we would have wished, and it is clear that further work on the methods is needed to develop quantitatively accurate procedures.

For the alkali-metal dimers and Sr we have chosen to proceed with a consistent set of full 4-component calculations of contact densities, using DFT calculations with the B3LYP functional Beck93 (); Lee:1988 (), which has a good track record for calculating changes in contact densities due to chemical bonding KFjctc2008 (). These were performed with the DIRAC package DIRAC11 (), since it is the only program used in this work which can offer a potentially exact treatment of relativity. The DIRAC calculations employed the aug-cc-pVTZ-DK basis set for Li Dunning () and the v3z basis sets of Dyall Djpca2009 () for other elements. The basic Gaussian finite nuclear charge distribution model, which has been shown to be sufficiently reliable to yield bond-induced changes in contact densities MLRcpl2008 (), was used with parameters given by Visscher and Dyall VDadndt1997 (). Contact densities were obtained by evaluating the expectation value . An ultrafine grid was employed to ensure converged results in the exchange-correlation evaluation.

For the Yb system DIRAC failed to converge and ADF ADF1 (); ADF2 () was used instead to evaluate . In this case scalar relativistic effects (which are the equivalent of Darwin and mass-velocity terms in the Breit-Pauli Hamiltonian) were included via the zero-order regular approximation (ZORA) CPDps1986 (); vLBSjcp1993 (); vLBSjcp1994 (); FSvLcpl1995 (). The ZORA/QZ4P basis set vLBjcc2003 (), which is an all-electron basis sets of Slater-type orbitals (STO), was used in conjunction with a point-charge nuclear model. ADF parameters were chosen to enable use of the true (exact) electron density in the exchange-correlation potential. As shown in the Supplementary Information, bonding-induced changes in field shifts obtained by this method cannot be regarded as quantitative, and might indeed be only order-of-magnitude estimates, but even this is valuable in understanding which effects are dominant for Yb.

The radial functions for bonding-induced changes in contact densities are shown in Fig. 5 for the alkali-metal dimers. The field shifts for LiLi, KK and RbRb are shown explicitly on the right-hand axes. Note that the changes in the nuclear charge radius are positive for LiLi ( fm) and KK ( fm), but negative for RbRb ( fm).

The overall magnitude of the bonding-induced changes in contact densities increases substantially from Li to K to Rb. This arises mostly because the non-relativistic contact density for an s electron is given approximately by the Goudsmit-Fermi-Segré (GFS) approximation Ko58 (); Ar71 (),

(6) |

where is the outer charge (which is 1 for Group I atoms), is the effective quantum number and is the quantum defect. The effective quantum number for the valence s electron increases only slowly down Group I, taking values of 1.588, 1.626, 1.771 and 1.805 for Li, Na, K and Rb, respectively vdSMbook (), while is small. This gives an overall contact density approximately proportional to . When electrons are treated relativistically, heavy atoms accumulate an additional scaling when the factor becomes non-negligible DFbook (); Moss (); BDNps2013 ().

The overall shapes of the curves in Figs. 5(a) to 5(c) can be explained in terms of simple bonding ideas. For the homonuclear alkali-metal dimers, the contact density changes are negative at larger distances, but positive at short range. This is qualitatively the same as the behavior for H Bader:1968 (). The negative values at longer range occurs partly because covalent bonding results in sp hybridisation, reducing the population in the s orbital and thus reducing the contact density. This effect is partly counteracted by contraction of the core in response to the reduced screening of the nucleus by p electrons. At short range, by contrast, the two atoms interpenetrate one another sufficiently to increase the density at both nuclei. For all the alkali-metal dimers, the overall contact density change is positive at the equilibrium distance, but with a substantial negative gradient.

For the heteronuclear dimers, there are additional effects from charge transfer, which reduce the density on the electropositive atom and increase the density on the electronegative one. These charge transfer effects counteract the covalent reduction slightly for Li in LiNa and overwhelm it for Li in the very polar molecules LiK, LiRb and LiCs. They also reinforce the covalent reduction slightly for Rb in KRb and substantially for K and Rb in LiK, LiRb, NaK and NaRb, and counteract it for K in KRb and KCs and for Rb in RbCs. All these charge transfer effects correlate reasonably well with the corresponding dipole moment functions Aymar:2005 ().

The situation is different for Sr, shown in Fig. 6(a). Here the contact density change due to bonding is positive at long range, but becomes negative around , where there is weak chemical bonding. The 4-component DFT/B3LYP/v3z calculations with DIRAC failed to converge at internuclear distances less than 8.7 bohr, so we supplemented them with ZORA/B3LYP/QZ4P calculations using ADF, shown in blue in Fig. 6(a). These show a weaker negative feature, probably a depth of only cm compared to the experimental value of 1081 cm, though also perhaps because of the incomplete treatment of relativity in ADF. At short range, the ADF contact densities turn upwards, due to interpenetration of the atomic densities. Similar short-range behavior has been seen theoretically in He Bader:1968 (), Ne and Ar Sebastian:1979 ().

For Yb, 4-component DFT/B3LYP/v3z calculations with DIRAC failed to converge entirely, so we used ZORA/B3LYP/QZ4P calculations with ADF instead. The resulting bonding-induced contact density changes are shown in Fig. 6(b). They are positive across the whole range of , without a negative region near . This may be simply because Yb shows weaker covalent bonding than Sr. It may be noted that Yb shows much stronger relativistic effects than Sr, and the contact density change for the lowest transition from multiconfiguration Dirac-Fock calculations Torbohm:1985 () is a factor of 5.3 larger. This is reflected in the overall magnitude of the bonding-induced changes to contact densities, which are much larger for Yb than for Sr.

The increase in the bonding contributions to field shifts in moving from Li to Yb is important, particularly when contrasted with the decrease in the mass shifts through the same series. The effect of bonding on Li field shifts is a few hundred kHz, while the mass shifts for the same systems are on the order of GHz. However, the effects of bonding on field shifts for K and Rb are a few MHz, while the mass shifts are on the order of tens of MHz. Thus, molecules formed from heavier alkali metals may have comparable bonding contributions to both field shifts and mass shifts, and studies of isotope shifts for these systems should consider both effects. The same comments apply to Sr. Continuing the trend, the effects of bonding on field shifts for Yb are computed to be at least two orders of magnitude larger than the corresponding mass shifts. The bonding changes in field shifts for this system are on the order of tens of MHz.

The functions shown in Figs. 4 and 5 indicate that there is indeed a crossover point in atomic mass where the effect of bonding on the field shift becomes larger than on the mass shift. This is analogous to the crossover between field and mass shifts for atoms BDNps2013 (), and should be considered in studies that consider molecular isotope shifts beyond the usual mass scaling.

### 3.4 Comparison with empirical isotope shifts

When characterizing a particular electronic state spectroscopically, it is sometimes possible to isolate small mass-dependent effects not accounted for by rudimentary mass-scaling. Such Born-Oppenheimer breakdown functions have been derived empirically for a number of relatively light diatomic molecules by least-squares fitting of measured line positions for a pair or series of isotopologs Coxon:HeH+:1999 (); Coxon:BeH+:1997 (); LeRoy:1999 (); Coxon:HF:2006 (); Coxon:HX:2015 (); Coxon:HCl:2000 (); LeRoy:2002 (); LeRoy:2005 (); Coxon:LiH:2004 (); LeRoy:2006 (); Henderson:2013 (); Coxon:CO:2004 (); WMLjcp2002 (); Coxon:2006 (); LeRoy:2009 (); Tiemann:2009 (); Ivanova:2011 (). The purpose of this section is to compare examples of Born-Oppenheimer breakdown functions from the literature with the results obtained using ab initio techniques. In particular, we will describe results for the Li, LiK, and LiRb systems, since these are the only alkali-metal diatomics for which we could find Born-Oppenheimer breakdown functions in the spectroscopic literature. Since in Section 3.3 the effects of chemical bonding on isotopic field shifts were found to be on the order of kHz for the substitution Li Li, we include only mass shifts in this section.

For the and states of Li, extensive sets of high-quality line positions have been measured for the LiLi, LiLi, and LiLi isotopologs, and have been used in several studies to derive correction functions for Born-Oppenheimer breakdown. Since Li is the most abundant isotope, it is used in the spectroscopic literature as the reference isotope. The LiLi isotopolog is the one for which the most line positions have been measured, so has the best-determined potential curve. In the following, the correction functions are given for the substitution LiLiLiLi.

Figure 7 compares our mass-shift function , computed at the RHF-CCSD/ANO-RCC level of theory, with the empirical Born-Oppenheimer breakdown corrections developed by Coxon and Melville Coxon:2006 () and by Le Roy et al. LeRoy:2009 (). All three functions are negative in the long-range region, with the function asymptotically approaching zero from below. However, the ab initio function changes sign near the potential minimum and is positive at short range, rising steeply between the potential minimum and the inner turning point at the dissociation energy. The empirical functions, by contrast, remain negative at short range.

In comparing the ab initio and empirical correction functions, it is important to appreciate that the ab initio function represents a true adiabatic correction, whereas the empirical functions in Fig. 7 are effective adiabatic corrections that contain contributions from both adiabatic and nonadiabatic terms. Watson Watson:1980 () showed that it is not possible to separate the adiabatic and nonadiabatic contributions on the basis of line positions alone. However, an alternative formulation of the Hamiltonian by Herman and Ogilvie Herman:1998 () does allow the contributions to be separated, using constraints from the molecular dipole moment function or rotational -factor. This approach has not yet been applied to alkali-metal dimers, but Coxon and Hajigeorgiou have shown that, for HCl Coxon:HCl:2000 () and CO Coxon:CO:2004 (), the empirical true and effective adiabatic corrections have similar values near but very different gradients (actually of opposite sign for HCl). This may be the origin of the qualitative differences in shape in Fig. 7.

We next consider LiK and LiRb, for which extensive spectra were measured in Refs. Tiemann:2009 () and Ivanova:2011 () using high-resolution fluorescence spectroscopy. The spectra were used to obtain ground-state potential energy curves and Born-Oppenheimer breakdown functions by least-squares techniques. Figure 8 compares the empirically determined Born-Oppenheimer breakdown functions with our ab initio mass shift functions, computed at the RHF-CCSD/ANO-RCC level of theory, with the Li isotopolog taken as the reference species in each case. The ab initio functions for LiK and LiRb have qualitatively similar features to that for Li, which in turn is qualitatively similar to that for H. For both systems, the ab initio and empirical functions have similar values at the equilibrium distance , but different gradients and different overall shapes. Tiemann et al. Tiemann:2009 () comment that the adiabatic correction near for LiK is about 5 times larger than the typical uncertainty in line positions at low , so is statistically well determined. As for Li, the difference in gradient between the ab initio and empirical functions around may plausibly be attributed to the fact that the empirical functions represent effective adiabatic corrections that include contributions from nonadiabatic effects.

At very long range, outside the Le Roy radius LeRoy:1973 (), the potential between atoms in S states dies off as . The dispersion coefficient is slightly different for different atomic isotopes, so should also be proportional to PJcp2012 (). For alkali-metal atoms the difference between the coefficients is determined mostly by the valence s electron, whose wavefunction dies off at long range as . Here is the effective nuclear charge and is an effective Bohr radius for an electron with reduced mass , where is the nuclear mass. An atom with a finite-mass nucleus has a larger , larger polarizability and larger (more attractive) coefficients that an atom with an infinite-mass nucleus. This corresponds to being negative at long range. The nonadiabatic contributions to the effective adiabatic correction are proportional to , so die off as and should not affect the long-range sign.

All the ab initio mass-shift functions are negative at very long range, but the empirical functions for LiK Tiemann:2009 () and LiRb Ivanova:2011 () are positive. In addition, should approach a constant outside . However, Fig. 4 of ref. Tiemann:2009 () shows that the empirical for LiK increases nearly linearly between and 23 bohr, with little sign of levelling off. It has reached only about 20% of its asymptotic value at 23 bohr, which is well outside the Le Roy radius. We therefore conclude that the long-range behaviour of the functional form used for the ground- and excited-state Born-Oppenheimer breakdown (BOB) functions should be revisited in future interpretations of the spectra.

## 4 Isotope shifts in ultracold molecular physics

The quantities that are usually measured in ultracold molecular physics are the binding energies of levels very close to dissociation, often as a function of magnetic field, and the positions of zero-energy Feshbach resonances as a function of magnetic field. The levels of interest are often bound by only a few MHz, which is less than 1 part in of the well depth for the alkali-metal dimers. Because of this, they are very strongly dominated by long-range effects, and are insensitive to the shape of the short-range potential. Nevertheless, the binding energies of these levels depend sensitively on the fractional part of the non-integer quantum number at dissociation LeRoy:1970 ().

Considerable insight into the effects of small potential shifts may be gained by writing the quantum number at dissociation semiclassically in terms of a WKB phase integral,

(7) |

where

(8) |

is the internuclear potential for a reference isotopolog and is the inner turning point at the dissociation energy. The usual WKB fraction is replaced by at dissociation if is asymptotically Gribakin:1993 (). The scattering length for a single potential curve is related to the phase integral by

(9) |

where is the mean scattering length of Gribakin and Flambaum Gribakin:1993 (), which depends only on and the reduced mass . The presence of in the integrand of Eq. 8 produces the normal Born-Oppenheimer mass scaling. If a mass-dependent perturbation is now introduced, changes by an additional amount

(10) |

This integral is formally problematic because so the condition is not satisfied very close to the turning point for the reference isotope, where is comparable to . In reality this should be handled by a shift of , but in practice this region makes little contribution to the integral and can be neglected.

Figure 9 compares the integrand of Eq. 10 for our mass-shift function for LiLiLiLi, together with its cumulative integral, with the corresponding quantities for the empirical functions of Coxon and Melville Coxon:2006 () and Le Roy et al. LeRoy:2009 (). It may be seen that all three give contributions to that are positive for this substitution, but that the two empirical functions give considerably larger values than the ab initio function. For comparison, the analysis of Julienne and Hutson Julienne:Li67:2014 () gave for LiLi and would thus give for LiLiLiLi. The ab initio function underestimates this by a factor of about 5, but this might be because it neglects nonadiabatic terms that contribute to the full effective adiabatic correction. In addition, it should be noted that the value of the integral is a delicate balance between short-range and long-range contributions.

Figure 10 shows similar plots for the semiclassical integrands and cumulative integrals in LiK and LiRb, comparing the ab initio mass-shift functions with the empirical ones of refs. Tiemann:2009 () and Ivanova:2011 (). Outside the Le Roy radius, should approach a constant, as described above, so that the semiclassical integrand should die off as . The integrands for the ab initio functions do show this behaviour, but those for the empirical functions remain nearly constant, so that the corresponding phase integrals show no sign of converging by bohr. This is because for the empirical functions increases nearly linearly in this region, and does not approach its asymptotic value until much larger distances.

It would be very valuable to obtain spectra of near-threshold levels or Feshbach resonance positions for LiK and LiRb for both Li and Li, in order to find values of that can be used, with the constraints on the long-range functions established here, to determine improved adiabatic correction functions.

## 5 Conclusions

We have investigated electronic structure calculations of bonding contributions to breakdown of the Born-Oppenheimer approximation for a range of molecules important in ultracold physics. These include the homonuclear and heteronuclear alkali-metal dimers and the Sr and Yb molecules. We have considered both isotopic mass shifts (also known as diagonal Born-Oppenheimer corrections, DBOCs, or adiabatic corrections) and isotopic field shifts (nuclear volume effects).

In a first step, we explored the performance of different electronic structure methods for the elementary systems H and Li. For these systems, it is well known that potential curves from restricted Hartree-Fock (RHF) calculations have incorrect dissociation behavior, whereas those based on unrestricted Hartree-Fock (UHF) reference calculations can have the correct behavior. However, for mass shifts we demonstrated inherent problems with methods based on UHF references. We suggest that, when computing mass shifts for systems that dissociate to open-shell monomers, unrestricted methods should be avoided whenever possible. When there is no alternative to unrestricted methods, care must be taken to avoid the occurrence of orbital instability envelopes.

For all our target molecules, we computed isotopic mass shifts at the CCSD level and field shifts from electron densities at the nucleus (contact densities) obtained with relativistic DFT methods. For many of the alkali-metal dimers, the mass-shift functions change sign in the vicinity of the equilibrium distance, as a result of competing physical effects dominating in the short- and long-range regions. It is thus difficult to identify periodic trends in the magnitude of the values near equilibrium. For the contact densities, the bonding changes are consistently positive at short range, but at longer range result from a combination of covalent effects (which are always negative) and charge-transfer effects, which are positive for the more electronegative atom in a molecule but negative for its partner.

The magnitudes of the mass shifts decrease with increasing atomic number, while the opposite is true for the field shifts. For Li, Na and K, mass shifts are strongly dominant. For Rb and Sr, mass shifts are still generally larger than field shifts, but the latter are not insignificant. For Yb, however, the field shifts are dominant.

For the light molecules Li, LiK, and LiRb, we have compared our ab initio mass-shift functions with Born-Oppenheimer breakdown functions fitted to electronic spectra. For LiK and LiRb the ab initio functions have similar values to the empirical functions at the equilibrium distance, where the empirical functions are most reliably determined. However, in all three cases the ab initio functions have qualitatively different shapes from the empirical functions away from equilibrium. This may be because the empirical functions are effective adiabatic corrections that include contributions from nonadiabatic terms.

The ab initio functions are slightly more attractive at long range for Li than for Li, as expected from the larger polarizability of Li. The empirical functions for LiK and LiRb have the opposite sign to the ab initio functions at long range. The results presented here should help inform the qualitative shape of the functional form used in future analyses of electronic spectra to model Born-Oppenheimer breakdown functions.

We also considered the effect of Born-Oppenheimer corrections on quantities of interest in ultracold physics. Scattering lengths and the positions of near-threshold levels may be related to the non-integer quantum number at dissociation. We developed a theory based on semiclassical phase integrals to give insight into how small perturbations affect this quantity. For Li, LiK and LiRb, the overall effect arises from a subtle balance of short-range and long-range effects. Neither the ab initio function nor the empirical functions for Li are in quantitative agreement with the overall mass shift obtained from studies of Feshbach resonances and near-threshold bound states. A simultaneous treatment of both types of experiment, incorporating insights from the ab initio studies, is needed to resolve the remaining discrepancies.

For molecules such as Sr and Yb, there are proposals to use deviations from Born-Oppenheimer mass scaling to probe a possible “fifth force” that may exist in addition to the familiar electromagnetic, gravitational and strong and weak nuclear forces. One possible force is a “short-range gravity”, proportional to the product of the nuclear masses in the molecule. However, before attributing any deviations from Born-Oppenheimer mass scaling to such forces, it is crucial to consider effects due to conventional isotopic mass shifts and field shifts. These have mass dependences different from short-range gravity but are likely to be difficult to distinguish in experiments. Nevertheless, if their effects can be calculated reliably, they can be taken into account, providing greater sensitivity to a fifth force (or allowing a tighter bound to be placed upon it). This work has provided an initial attempt to investigate the magnitude of mass and field shifts in Sr and Yb.

Data underlying this article are available at http://dx.doi.org/10.15128/r2vd66vz89.

## Acknowledgments

This paper is dedicated to Robert J. Le Roy, who was JMH’s postdoctoral supervisor in 1981-83 and had a formative influence on his early scientific career. Anyone who has had a draft paper “Le Royed” will appreciate its importance.

The authors are grateful for the use of the EPSRC UK National Service for Computational Chemistry Software (NSCCS) at Imperial College London, awarded as NSCCS Project CHEM700. We thank Paul Julienne for initiating this work and John Coxon for bringing the importance of nonadiabatic effects to our attention. We are also grateful to John Stanton, Péter Szalay, and Markus Reiher for helpful discussions regarding specific functionalities, keywords, and limitations associated with their software. This work was supported by the Engineering and Physical Sciences Research Council [grant number EP/I012044/1].

## References

- (1) W. Kołos and L. Wolniewicz, J. Chem. Phys. 41, 3663 (1964).
- (2) W. Kołos and L. Wolniewicz, J. Chem. Phys. 43, 2429 (1965).
- (3) G. Herzberg, J. Mol. Spect. 33, 147 (1970).
- (4) O.L. Polyansky and J. Tennyson, J. Chem. Phys. 110, 5056 (1999).
- (5) J.A. Coxon and P.G. Hajigeorgiou, J. Mol. Spectrosc. 193, 306 (1999).
- (6) J.A. Coxon and R. Colin, J. Mol. Spectrosc. 181, 215 (1997).
- (7) R.J. Le Roy, J. Mol. Spectrosc. 194, 198 (1999).
- (8) J.A. Coxon and P.G. Hajigeorgiou, J. Phys. Chem. A 110, 6261 (2006).
- (9) J.A. Coxon and P.G. Hajigeorgiou, J. Quant. Spectrosc. Rad. Transf. 151, 133 (2015).
- (10) J.A. Coxon and P.G. Hajigeorgiou, J. Mol. Spectrosc. 203, 49 (2000).
- (11) R.J. Le Roy and Y.Y. Huang, J. Mol. Struct. (Theochem) 591, 175 (2002).
- (12) R.J. Le Roy, D.R.T. Appadoo, K. Anderson, A. Shayesteh, I.E. Gordon and P.F. Bernath, J. Chem. Phys. 123, 204304 (2005).
- (13) J.A. Coxon and C.S. Dickinson, J. Chem. Phys. 121, 9378 (2004).
- (14) R.J. Le Roy, D.R.T. Appadoo, R. Colin and P.F. Bernath, J. Mol. Spec. 236, 178 (2006).
- (15) R.D.E. Henderson, A. Shayesteh, J. Tao, C.C. Haugen, P.F. Bernath and R.J. Le Roy, J. Phys. Chem. A 117, 13373 (2013).
- (16) J.A. Coxon and P.G. Hajigeorgiou, J. Chem. Phys. 121, 2992 (2004).
- (17) X.J. Wang, J. Magnes, A.M. Lyyra, A.J. Ross, F. Martin, P.M. Dove and R.J. Le Roy, J. Chem. Phys. 117, 9339 (2002).
- (18) J.A. Coxon and T.C. Melville, J. Mol. Spectrosc. 235 (2), 235 (2006).
- (19) R.J. Le Roy, N.S. Dattani, J.A. Coxon, A.J. Ross, P. Crozet and C. Linton, J. Chem. Phys. 131, 204309 (2009).
- (20) E. Tiemann, H. Knöckel, P. Kowalczyk, W. Jastrzebski, A. Pashov, H. Salami and A.J. Ross, Phys. Rev. A 79 (4), 042716 (2009).
- (21) M. Ivanova, A. Stein, A. Pashov, H. Knöckel and E. Tiemann, J. Chem. Phys. 134, 024321 (2011).
- (22) J.Y. Seto, R.J. Le Roy, J. Vergès and C. Amiot, J. Chem. Phys 113 (8), 3067 (2000).
- (23) O. Docenko, M. Tamanis, R. Ferber, E.A. Pazyuk, A. Zaitsevskii, A.V. Stolyarov, A. Pashov, H. Knöckel and E. Tiemann, Phys. Rev. A 75 (4), 042503 (2007).
- (24) C. Strauss, T. Takekoshi, F. Lang, K. Winkler, R. Grimm, J. Hecker Denschlag and E. Tiemann, Phys. Rev. A 82, 052514 (2010).
- (25) S. Falke, H. Knöckel, J. Friebe, M. Riedmann, E. Tiemann and C. Lisdat, Phys. Rev. A 78, 012503 (2008).
- (26) B.J. Verhaar, E.G.M. van Kempen and S.J.J.M.F. Kokkelmans, Phys. Rev. A 79, 032711 (2009).
- (27) H. Knöckel, B. Bodermann and E. Tiemann, Eur. Phys. J. D 28, 199 (2004).
- (28) E.J. Salumbides, K.S.E. Eikema, W. Ubachs, U. Hollenstein and E. Tiemann, Eur. Phys. J. D 47, 171 (2008).
- (29) C. Chin, R. Grimm, E. Tiesinga and P.S. Julienne, Rev. Mod. Phys. 82, 1225 (2010).
- (30) A. Simoni, M. Zaccanti, C. D’Errico, M. Fattori, G. Roati, M. Inguscio and G. Modugno, Phys. Rev. A 77, 052705 (2008).
- (31) H.W. Cho, D.J. McCarron, M.P. Köppinger, D.L. Jenkin, K.L. Butler, P.S. Julienne, C.L. Blackley, C.R. Le Sueur, J.M. Hutson and S.L. Cornish, Phys. Rev. A 87, 010703(R) (2013).
- (32) C.L. Blackley, C.R. Le Sueur, J.M. Hutson, D.J. McCarron, M.P. Köppinger, H.W. Cho, D.L. Jenkin and S.L. Cornish, Phys. Rev. A 87, 033611 (2013).
- (33) P.S. Julienne and J.M. Hutson, Phys. Rev. A 89, 052715 (2014).
- (34) M. Kitagawa, K. Enomoto, K. Kasa, Y. Takahashi, R. Ciuryło, P. Naidon and P.S. Julienne, Phys. Rev. A 77, 012719 (2008).
- (35) T. Zelevinsky, M.M. Boyd, A.D. Ludlow, T. Ido, J. Ye, R. Ciuryło, P. Naidon and P.S. Julienne, Phys. Rev. Lett. 96, 203201 (2006).
- (36) S. Stellmer, B. Pasquiou, R. Grimm and F. Schreck, Phys. Rev. Lett. 109, 115302 (2012).
- (37) B.H. McGuyer, C.B. Osborn, M. McDonald, G. Reinaudi, W. Skomorowski, R. Moszynski and T. Zelevinsky, Phys. Rev. Lett. 111, 243003 (2013).
- (38) M. Borkowski, P. Morzyński, R. Ciuryło, P.S. Julienne, M. Yan, B.J. DeSalvo and T.C. Killian, Phys. Rev. A 90, 032713 (2014).
- (39) B.H. McGuyer, M. McDonald, G.Z. Iwata, M.G. Tarallo, W. Skomorowski, R. Moszynski and T. Zelevinsky, Nat. Phys. 11, 32 (2014).
- (40) E.J. Salumbides, J.C.J. Koelemeij, J. Komasa, K. Pachucki, K.S.E. Eikema and W. Ubachs, Phys. Rev. D 87, 112008 (2013).
- (41) M.L. Niu, E.J. Salumbides, G.D. Dickenson, K.S.E. Eikema and W. Ubachs, J. Mol. Spectrosc. 300, 44 (2014).
- (42) E.J. Salumbides, W. Ubachs and V.I. Korobov, J. Mol. Spectrosc. 300, 65 (2014).
- (43) M. Born and K. Huang, Dynamical Theory of Crystal Lattices (Oxford University, New York, 1956).
- (44) M. Born and J.R. Oppenheimer, Ann. Phys. 84, 457 (1927).
- (45) E. Gaidamauskas, C. Nazé, P. Rynkun, G. Gaigalas, P. Jönsson and M. Godefroid, J. Phys. B: At. Mol. Opt. Phys. 44, 175003 (2011).
- (46) C. Nazé, E. Gaidamauskas, G. Gaigalas, M. Godefroid and P. Jönsson, Comp. Phys. Comm. 184, 2187 (2013).
- (47) P. Jönsson, G. Gaigalas, J. Bieroń, C.F. Fischer and I.P. Grant, Comp. Phys. Comm. 184, 2197 (2013).
- (48) W. Kutzelnigg, Mol. Phys. 90, 909 (1997).
- (49) H. Sellers and P. Pulay, Chem. Phys. Lett. 103, 463 (1984).
- (50) N.C. Handy, Y. Yamaguchi and H.F. Schaefer, III, J. Chem. Phys. 84, 4481 (1986).
- (51) H. Sellers, Chem. Phys. Lett. 108, 339 (1984).
- (52) N.C. Handy and A.M. Lee, Chem. Phys. Lett. 252, 425 (1996).
- (53) B.C. Garrett and D.G. Truhlar, J. Chem. Phys. 82, 4543 (1985).
- (54) A.G. Ioannou, R.D. Amos and N.C. Handy, Chem. Phys. Lett. 251, 52 (1996).
- (55) N.C. Handy, T.J. Lee and W.H. Miller, Chem. Phys. Lett. 125, 12 (1986).
- (56) R.D. Bardo and M. Wolfsberg, J. Chem. Phys. 68, 2686 (1978).
- (57) N.F. Zobov, O.L. Polyansky, C.R. Le Sueur and J. Tennyson, Chem. Phys. Lett. 260, 381 (1996).
- (58) O.L. Polyansky, N.F. Zobov, S. Viti, J. Tennyson, P.F. Bernath and L. Wallace, Science 277, 346 (1997).
- (59) O.L. Polyansky, A.G. Császár, S.V. Shirin, N.F. Zobov and Barletta, Science 299, 539 (2003).
- (60) A. Tajti, P.G. Szalay, A.G. Császár, M. Kállay, J. Gauss, E.F. Valeev, B.A. Flowers, J. Vázquez and J.F. Stanton, J. Chem. Phys. 121, 11599 (2004).
- (61) Y.J. Bomble, J. Vázquez, M. Kállay, C. Michauk, P.G. Szalay, A.G. Császár, J. Gauss and J.F. Stanton, J. Chem. Phys. 125, 064108 (2006).
- (62) A. Karton, E. Rabinovich, J.M.L. Martin and B. Ruscic, J. Chem. Phys. 125, 144109 (2006).
- (63) A. Karton and J.M.L. Martin, J. Chem. Phys. 136, 124114 (2012).
- (64) E.F. Valeev and C.D. Sherrill, J. Chem. Phys. 118, 3921 (2003).
- (65) B. Temelso, E.F. Valeev and C.D. Sherrill, J. Phys. Chem. A 108, 3068 (2004).
- (66) S. Garashchuk, J.C. Light and V.A. Rassolov, Chem. Phys. Lett. 333, 459 (2001).
- (67) CFOUR, a quantum chemical program package, J. F. Stanton, J. Gauss, M. E. Harding, et al., see http://www.cfour.de.
- (68) J. Gauss, A. Tajti, M. Kállay, J.F. Stanton and P.G. Szalay, J. Chem. Phys. 125, 144111 (2006).
- (69) A. Tajti, P.G. Szalay and J. Gauss, J. Chem. Phys 127, 014102 (2007).
- (70) D.W. Schwenke, J. Phys. Chem. A 105, 2352 (2001).
- (71) D.W. Schwenke, J. Chem. Phys. 118, 6898 (2003).
- (72) S.L. Mielke, D.W. Schwenke and K.A. Peterson, J. Chem. Phys. 122, 224313 (2005).
- (73) Y. Imafuku, M. Abe, M.W. Schmidt and M. Hada, J. Phys. Chem. A 120, 2150 (2016).
- (74) J.C. Tully, Nonadiabatic Processes in Molecular Collisions, in Dynamics of Molecular Collisions, edited by W. H. Miller (Plenum, New York, 1976), pp. 217–267.
- (75) IUPAC. Compendium of Chemical Terminology, 2nd ed. (the “Gold Book”). Compiled by A. D. McNaught and A. Wilkinson. Blackwell Scientific Publications, Oxford (1997).
- (76) M. Baer, Beyond Born-Oppenheimer: Electronic Nonadiabatic Coupling Terms and Conical Intersections (Wiley-Blackwell, New Jersey, 2006).
- (77) K. Pachucki and J. Komasa, J. Chem. Phys. 129, 034102 (2008).
- (78) K. Pachucki and J. Komasa, J. Chem. Phys. 130, 164113 (2009).
- (79) P.R. Bunker and R.E. Moss, Mol. Phys. 33, 417 (1977).
- (80) J.K.G. Watson, J. Mol. Spectrosc. 80, 411 (1980).
- (81) P.R. Bunker, C.J. McLarnon and R.E. Moss, Mol. Phys. 33, 425 (1977).
- (82) W. Kutzelnigg, Mol. Phys. 105, 2627 (2007).
- (83) D.S. Hughes and C. Eckart, Phys. Rev. 36, 694 (1930).
- (84) B. Fricke, Phys. Rev. Lett. 30, 119 (1973).
- (85) E.W. Otten, Nuclear Radii and Moments of Unstable Isotopes, in Treatise on Heavy-Ion Science, edited by D. A. Bromley, Vol. 8 (Plenum Press, New York, 1989), pp. 517–638.
- (86) K. Blaum, J. Dilling and W. Nörtershäuser, Phys. Scripta T153, 014017 (2013).
- (87) J.E. Rosenthal and G. Breit, Phys. Rev. 41, 459 (1932).
- (88) G. Racah, Nature (Lond.) 129, 723 (1932).
- (89) I. Angeli, Table of Nuclear Root Mean Square Radii, INDC(HUN)-033, International Atomic Energy Agency, Vienna 1999 .
- (90) P. Aufmuth, K. Heilig and A. Steudel, At. Data Nucl. Data Tables 37, 455â490 (1987).
- (91) I. Angeli and K.P. Marinova, At. Data Nucl. Data Tables 99, 65 (2013).
- (92) E. Tiemann, H. Arnst, W.U. Stieda, T. Törring and J. Hoeft, Chem. Phys. 67, 133 (1982).
- (93) J. Schlembach and E. Tiemann, Chem. Phys. 68, 21 (1982).
- (94) J. Knöckel and E. Tiemann, Chem. Phys. 68, 13 (1982).
- (95) H. Knöckel, T. Kröckertskothen and E. Tiemann, Chem. Phys. 93, 349 (1985).
- (96) S.A. Cooke and M.C.L. Gerry, Phys. Chem. Chem. Phys. 6, 4579 (2004).
- (97) S.A. Cooke, J.M. Michaud and M.C.L. Gerry, J. Mol. Spec. 695, 13 (2004).
- (98) S.A. Cooke, M.C.L. Gerry and D.P. Chong, Chem. Phys. 298, 205 (2004).
- (99) S.A. Cooke, C. Krumrey and M.C.L. Gerry, Phys. Chem. Chem. Phys. 7, 2570 (2005).
- (100) K.C. Etchison, C.T. Dewberry and S.A. Cooke, Chem. Phys. 342, 71 (2007).
- (101) M.M. Serafin, S.A. Peebles, C.T. Dewberry, K.C. Etchison, G.S. Grubbs, II, R.A. Powoski and S.A. Cooke, Chem. Phys. Lett. 449, 33 (2007).
- (102) B.M. Giuliano, L. Bizzocchi, S. Cooke, D. Banser, M. Hess, J. Fritzsche and J.U. Grabow, Phys. Chem. Chem. Phys. 10, 2078 (2008).
- (103) C.T. Dewberry, K.C. Etchison, G.S. Grubbs, II, R.A. Powoski, M.M. Serafin, S.A. Peebles and S.A. Cooke, J. Mol. Spec. 248, 20 (2008).
- (104) C. Krumrey, S.A. Cooke, D.K. Russell and M.C.L. Gerry, Can. J. Phys. 87, 567 (2009).
- (105) B.M. Giuliano, L. Bizzocchi, R. Sanchez, P. Villanueva, V. Cortijo, M.E. Sanz and J.U. Grabow, J. Phys. Chem. 135, 084303 (2011).
- (106) S. Knecht and T. Saue, Chem. Phys. 401, 103 (2012).
- (107) S. Fux and M. Reiher, Electron Density in Quantum Theory, in Electron Density and Chemical Bonding II: Theoretical Charge Density Studies, edited by D. Stalke, Vol. 147 (Springer, Berlin, 2012), pp. 99–142.
- (108) M. Filatov, Coord. Chem. Rev. 253, 594 (2009).
- (109) R. Kurian and M. Filatov 4, 278 (2008).
- (110) M. Filatov, J. Chem. Phys. 127, 084101 (2007).
- (111) V.N. Nemykin and R.C. Hadt, Inorg. Chem. 45, 8297 (2006).
- (112) F. Neese, Inorg. Chim. Acta 337, 181 (2002).
- (113) T. Lovell, W.G. Han, T. Liu, and L. Noodelman, J. Am. Chem. Soc. 124, 5890 (2002).
- (114) T. Lovell, J. Li, T. Liu, D.A. Case and L. Noodelman, J. Am. Chem. Soc. 123, 12392 (2001).
- (115) V.R. Marathe and A. Trautwein, in Advances in Mössbauer Spectroscopy, edited by B. V. Thosar, J. K. Srivasta, P. K. Iyengar and S. C. Bhargava (Elsevier, Amsterdam, 1985), p. 398.
- (116) W.C. Nieuwpoort, D. Post and P.T. van Duijnen, Phys. Rev. B 17, 91 (1978).
- (117) J.L.K.F. de Vries, J.M. Trooster and P. Ros, J. Chem. Phys. 63, 5256 (1975).
- (118) D.A. Shirley, Rev. Mod. Phys. 36, 339 (1964).
- (119) S. Knecht, S. Fux, R. van Meer, L. Visscher, M. Reiher and T. Saue, Theor. Chem. Acc. 129, 631 (2011).
- (120) M. Filatov, W.L. Zou and D. Cremer, J. Chem. Theory Comput. 8, 875 (2012).
- (121) W. Liu, Mol. Phys. 108, 1679 (2010).
- (122) T. Saue, ChemPhysChem 12, 3077 (2011).
- (123) P. Pyykkö, Chem. Rev. 112, 371 (2012).
- (124) J. Autschbach, J. Chem. Phys. 136, 150902 (2012).
- (125) W. Liu, Phys. Rep. 537, 59 (2014).
- (126) R. Mastalerz, R. Lindh and M. Reiher, Chem. Phys. Lett. 465, 157 (2008).
- (127) R. Mastalerz, P.O. Widmark, B.O. Roos, R. Lindh and M. Reiher, J. Chem. Phys. 133, 144111 (2010).
- (128) A.C. Neto, E.P. Muniz, R. Centoducatte and F.E. Jorge, J. Mol. Struct. (Theochem) 718, 219 (2005).
- (129) G.G. Camiletti, S.F. Machado and F.E. Jorge, J. Comp. Chem. 29, 2434 (2008).
- (130) C.L. Barros, P.J.P. de Oliveira, F.E. Jorge, A.C. Neto and M. Campos, Mol. Phys. 108, 1965 (2010).
- (131) A.C. Neto and F.E. Jorge, Chem. Phys. Lett. 582, 158 (2013).
- (132) S. Huzinaga and B. Miguel, Chem. Phys. Lett. 175, 289 (1990).
- (133) S. Huzinaga and M. Klobukowski, Chem. Phys. Lett. 212, 260 (1993).
- (134) O. Roos, V. Veryazov and P.O. Widmark, Theor. Chem. Acc. 111, 345 (2003).
- (135) P.O. Widmark, P.A. Malmqvist and B.O. Roos, Theor. Chim. Acta 77, 291 (1990).
- (136) L.K. Sørensen, T. Fleig and J. Olsen, J. Phys. B 42, 165102 (2009).
- (137) J.F. Stanton, J. Chem. Phys. 101, 371 (1994).
- (138) A.I. Krylov, J. Chem. Phys. 113, 6052 (2000).
- (139) H.B. Schlegel, Spin Contamination, in Encyclopedia of Computational Chemistry, edited by P. v.R. Schleyer, N. L. Allinger, T. Clark, J. Gasteiger, P. A. Kollman, H. F. Schaefer and P. R. Schreiner, Vol. 4 (John Wiley & Sons, New York, 2002).
- (140) J.L. Sonnenberg, H.B. Schlegel and H.P. Hratchian, Spin Contamination in Inorganic Chemistry Calculations, in Computational Inorganic and Bioinorganic Chemistry, edited by E. I. Solomon, R. A. Scott and R. B. King (John Wiley & Sons, New York, 2009), p. 173.
- (141) F. King, J. Mol. Struct. THEOCHEM 400, 7 (1997).
- (142) M. Puchalski and K. Pachucki, Phys. Rev. A 73, 022503 (2006).
- (143) W. Kołos and J. Rychlewski, J. Chem. Phys. 98, 3960 (1993).
- (144) T.D. Crawford, J.F. Stanton, W.D. Allen and H.F. Schaefer, J. Chem. Phys. 107, 10626 (1997).
- (145) A.A. Auer and J. Gauss, Chem. Phys. 356, 7 (2009).
- (146) D.A. Fedorov, A. Derevianko and S.A. Varganov, J. Chem. Phys. 140, 184315 (2014).
- (147) S. Kotochigova, J. Chem. Phys. 128, 024303 (2008).
- (148) N.S. Mosyagin, A.N. Petrov, A.V. Titov, A.V. Zaitsevskii and E.A. Rykova, Int. J. Quant. Chem. 111, 3793 (2009).
- (149) M. Aymar and O. Dulieu, J. Chem. Phys. 122 (20), 204302 (2005).
- (150) A. Stein, H. Knöckel and E. Tiemann, Phys. Rev. A 78, 042508 (2008).
- (151) A. Szabo and N.S. Ostlund, Modern Quantum Chemistry: An Introduction to Electronic Structure Theory (Dover, New York, 1996), p. 205.
- (152) A.D. Becke, J. Chem. Phys. 98, 5648 (1993).
- (153) C. Lee, W. Yang and P.R. G., Phys. Rev. B 37, 785 (1988).
- (154) DIRAC, a relativistic ab initio electronic structure program, Release DIRAC11 (2011), written by R. Bast, H. J. Aa. Jensen, T. Saue, and L. Visscher, with contributions from V. Bakken, K. G. Dyall, S. Dubillard, U. Ekström, E. Eliav, T. Enevoldsen, T. Fleig, O. Fossgaard, A. S. P. Gomes, T. Helgaker, J. K. Lærdahl, J. Henriksson, M. Iliaš, Ch. R. Jacob, S. Knecht, C. V. Larsen, H. S. Nataraj, P. Norman, G. Olejniczak, J. Olsen, J. K. Pedersen, M. Pernpointner, K. Ruud, P. Sałek, B. Schimmelpfennig, J. Sikkema, A. J. Thorvaldsen, J. Thyssen, J. van Stralen, S. Villaume, O. Visser, T. Winther, and S. Yamamoto (see http://dirac.chem.vu.nl).
- (155) B.P. Prascher, D.E. Woon, K.A. Peterson, T.H. Dunning, Jr. and A.K. Wilson, Theor. Chem. Acc. 128, 69 (2011).
- (156) K.G. Dyall, J. Phys. Chem. A 113, 12638 (2009).
- (157) L. Visscher and K.G. Dyall, At. Data Nucl. Data Tables 67, 207 (1997).
- (158) G. te Velde, F.M. Bickelhaupt, S.J.A. van Gisbergen, C. Fonseca Guerra, E.J. Baerends, J.G. Snijders and T. Ziegler, J. Comput. Chem. 22, 931 (2001).
- (159) C. Fonseca Guerra, J.G. Snijders, G. te Velde and E.J. Baerends, Theor. Chem. Acc. 99, 391 (1998).
- (160) C. Chang, M. Pelissier and M. Durand, Phys. Scripta 34, 394 (1986).
- (161) E. van Lenthe, E.J. Baerends and J.G. Snijders, J. Chem. Phys. 99, 4597 (1993).
- (162) E. van Lenthe, E.J. Baerends and J.G. Snijders, J. Chem. Phys. 101, 9783 (1994).
- (163) S. Faas, J.G. Snijders, J.H. van Lenthe, E. van Lenthe and E.J. Baerends, Chem. Phys. Lett. 246, 632 (1995).
- (164) E. van Lenthe and E.J. Baerends, J. Comp. Chem. 24, 1142 (2003).
- (165) H. Kopfermann, Nuclear Moments (Academic Press, New York, 1958), pp. 161–188.
- (166) L. Armstrong, Theory of the Hyperfine Structure of Free Atoms (Wiley (Interscience), New York, 1971).
- (167) P. van der Straten and H. Metcalf, Atoms and Molecules Interacting with Light: Atomic Physics for the Laser Era (Cambridge U.P., Cambridge, 2016), p. 173.
- (168) K.G. Dyall and J. K. Fægri, Introduction to Relativistic Quantum Chemistry (Oxford University, New York, 2007).
- (169) R. Moss, Advanced Molecular Quantum Mechanics (Springer, Netherlands, 1973), pp. 194–224.
- (170) R.F.W. Bader and A.K. Chandra, Can. J. Chem 46, 953 (1968).
- (171) K.L. Sebastian and A.K. Chandra, Pramana 12, 9 (1979).
- (172) G. Torbohm, B. Fricke and A. Rosén, Phys. Rev. A 31, 2038 (1985).
- (173) R.M. Herman and J.F. Ogilvie, Adv. Chem. Phys. 103, 187 (1998).
- (174) R.J. Le Roy, in Molecular Spectroscopy, Vol. 1, Specialist Periodical Reports (The Chemical Society, London, 1973), pp. 113–176.
- (175) M. Przybytek and B. Jeziorski, Chem. Phys. 401, 170 (2012).
- (176) R.J. Le Roy and R.B. Bernstein, J. Chem. Phys. 52, 3869 (1970).
- (177) G.F. Gribakin and V.V. Flambaum, Phys. Rev. A 48, 546 (1993).