Table Of ContentNonadiabatic corrections to rovibrational levels of H
2
Krzysztof Pachucki∗
Institute of Theoretical Physics, University of Warsaw, Hoz˙a 69, 00-681 Warsaw, Poland
Jacek Komasa†
Faculty of Chemistry, A. Mickiewicz University, Grunwaldzka 6, 60-780 Poznan´, Poland
(Dated: January 26, 2009)
Theleadingnonadiabaticcorrectionstorovibrationallevelsofadiatomicmoleculeareexpressedin
9 termsofthreefunctionsofinternucleardistance: correctionstotheadiabaticpotential,theeffective
0 nuclear mass, and the effective moment of inertia. The resulting radial Schr¨odinger equation for
0 nuclear motion is solved numerically yielding accurate nonadiabatic energies for all rovibrational
2
levels of H2 molecule in excellent agreement with previouscalculations byWolniewicz.
n
a PACSnumbers: 31.15.-p,31.15.ac,31.50.-x
J
4
2 I. INTRODUCTION der of the perturbative expansion, an effective nuclear
Hamiltonian for the ground electronic state of diatomic
] molecules, in terms of the nonadiabatic potential, vibra-
In the fully nonadiabatic approach the total nonrela-
h
tional and rotational masses. Assuming constant vibra-
p tivistic energy of a molecular state is obtained by solv-
tional and rotational masses and neglecting the nona-
- ing the Schr¨odingerequation with kinetic energy of elec-
m diabatic potential, they obtained values of these effec-
trons and of nuclei on the same footing. This approach
tive masses for H and D by fitting to the experimental
e hasbeenappliedtovibrationalstatesofseveralsmalldi- 2 2
h atomic molecules [1, 2, 3, 4, 5]. Much more commonly data. Later Schwenke in [13] used the Bunker and Moss
c [9, 14] effective Hamiltonian to perform ab initio calcu-
though, the total energy of a molecular state is obtained
.
s inatwo-stepprocedurebasedonthe Born-Oppenheimer lations of nonadiabatic corrections for H2 and H2O. His
c results for purely vibrational spectrum of H differ from
(BO) approximation [6, 7] in which a separation of elec- 2
i
s tronic and nuclear motion is assumed. Namely, in the that of Wolniewicz [15] by about 20%due to the inaccu-
y rate numerical representation of the wave function and,
first step, the electronic Schr¨odinger equation with the
h
what we demonstrate in this work, due to the approxi-
p clamped nuclei Hamiltonian is solved for different nu-
mate secondordernonadiabatic potential ofBunker and
[ clear configurations yielding the electronic energy as a
Moss [9]. The accuracy of Schwenke calculations [13] for
functionofthenuclearcoordinates. Thisfunction,called
2 H O is probably not higher, but clearly demonstrates
the potential energy surface (PES), serves as a potential 2
v
wide applicability of the perturbative approach. Very
for the motionofnucleiin the nuclearSchr¨odingerequa-
5 recently Kutzelnigg et al. [11, 12] performed simplified
5 tion. Solving this equation in the second step yields the
calculations of the nonadiabatic correction and both ef-
3 total rovibronic energy of the system.
4 The accuracy of theoretical predictions for molecular fectivemassesasfunctionsoftheinternucleardistancein
1. states, limited due to the BO approximation, can be in- H+2 and H2. In our recent work [16] we have introduced
nonadiabatic perturbative theory and derived formulae
1 creasedbyincludingtheadiabatic,relativistic,andradia-
8 tive corrections without renouncing the notion of PES. for the leading nonadiabatic corrections to energies and
0 wavefunctions. Theformulaforthenonadiabaticenergy
Further increase in the accuracy requires the nonadia-
: from that work, although apparently different, is in fact
v baticeffectstobetakenintoaccount. Adesirablewayof
equivalenttothatofBunkerandMoss[9]. Ourresultsfor
i inclusionoftheseeffectsisintermsofageometrydepen-
X rotationless vibrationalstates have been obtained as the
dent function, which can be added to PES in the same
r expectation value of nonadiabatic corrections with the
manner as all the other corrections. On one hand, the
a
adiabaticwavefunction. Althoughnumericallyaccurate,
nuclear Schr¨odinger equation, when solved with such a
duetotheneglectedthirdordernonadiabaticcorrections
potential,givesthe molecularenergylevelswithspectro-
[see Eq. (47)], our results differed by about 2% from the
scopic precision. On the other, the notion of the PES is
previous calculations by Wolniewicz [15] and by Stanke
preserved with all its advantages.
et al. [17].
Several, more or less successful attempts to construct
such a nonadiabatic correction function for a diatomic In this paper, the nonadiabatic perturbation theory
has been extended in two directions. Firstly, we general-
molecule can be found in literature [8, 9, 10, 11, 12].
izethepreviousderivationtorotationalstates. Secondly,
Bunker and Moss have derived [9], in the second or-
we include the previously missing third order correction,
which has proved significant. Moreover, we present a
rigorous formulation of the nonadiabatic perturbative
∗Electronicaddress: [email protected] theory and include the numerical example of the H2
†Electronicaddress: [email protected] molecule. Thiscanbeextendedtoanydiatomicmolecule
2
and potentially to an arbitrarily large molecule. We de- Intheadiabaticapproximationthetotalwavefunction
rive formulae valid to all orders, present the leading cor- of an arbitrary molecule
rections of order O(µ−2), and express them in terms of
n
thenonadiabaticcorrectiontothepotentialandtheeffec- φ (~r,R~)=φ (~r)χ(R~) (7)
a el
tive R-dependent nuclear mass and the moment of iner-
tia. These three functions enter the nuclear Schr¨odinger is represented as a product of the electronic wave func-
equation, which can be solved numerically for an arbi- tion φel and the nuclear wave function χ. We note, that
trary energy level. As a test of the presented perturba- φ dependsimplicitlyonthenuclearcoordinatesR~. The
el
tive theory, we perform calculations of all 301 rovibra- electronic wave function obeys the clamped nuclei elec-
tional levels of H molecule. We find an excellent 0.1% tronic Schr¨odinger equation
2
agreementwiththeaccuratenonadiabaticcorrectionsfor
states with the angular momentum J ≤ 10, which were H −E (R~) |φ i=0, (8)
el el el
obtained by Wolniewicz in [15], and present for the first
(cid:2) (cid:3)
time results for states with J >10. while the nuclear wave function is a solution to the
Schr¨odingerequationin the effective potentialgenerated
by electrons
II. THE ADIABATIC APPROXIMATION
H +E (R~)+E (R~)−E |χi=0, (9)
n a el a
The total wave function φ is the solution of the sta-
where (cid:2) (cid:3)
tionary Schr¨odinger equation
E (R~)= φ |H |φ . (10)
[H −E]|φi=0, (1) a el n el el
For the diatomic molecu(cid:10)le the nucl(cid:11)ear radial equation
with the Hamiltonian
reads
H =Hel+Hn, (2) 1 ∂ R2 ∂ J(J +1)
− +
2R2 ∂R µ ∂R 2µ R2
split into the electronic and nuclear parts. In the elec- (cid:20) n n
tronic Hamiltonian H
el +E (R)+E (R)−E χ (R) = 0, (11)
a el a J
∇2 (cid:21)
H =− a +V (3)
el 2m where J is the rotational quantum number.
a e
X
nuclear masses are, by definition, set to infinity, and the
III. PERTURBATIVE FORMALISM
potential V includes all the Coulomb interactions with
fixed positions R~ of nuclei. The nuclear Hamiltonian
A
involves kinetic energies of all nuclei The total wave function
H =− ∇R2A . (4) φ=φa+δφna =φelχ+δφna (12)
n
2M
A A is the sum of the adiabatic solution and a nonadiabatic
X
correction. The nonadiabatic correction δφ is decom-
na
Theseparationofcenterofmassmotionandthechoiceof
posed into two parts
the reference frame depends on the particular molecule.
For example, for a diatomic molecule in the space fixed δφ =φ δχ+δ′φ , (13)
na el na
referenceframeattachedtothegeometricalcenteroftwo
nuclei Hn takes the form which obey the following orthogonality conditions
H =− ∇R2 − ∇2el − 1 − 1 ∇~ ·∇~ , (5) hδ′φna|φeliel = 0, (14)
n R el
2µ 2µ M M
n n (cid:18) B A(cid:19) hδχ|χi = 0. (15)
where
Thelastequationmeansthatnormalizationofφisofthe
1 form
∇~ ≡ ∇~ , (6)
el a
2
a hφ χ|φi=1. (16)
X el
R~ =R~AB =R~A−R~B, and 1/µn =1/MA+1/MB is the The total energy
nuclear reduced mass. The last term in Eq. (5) vanishes
for homonuclear diatomic molecules. E =E +δE (17)
a na
3
is the sum of the adiabatic energy E and the nonadi- It is not always true, especially for rovibrational levels
a
abatic correction δE . Using above definitions we pro- close to the dissociation threshold. In spite of this fact,
na
ceed with the derivation of perturbative formulae. weclaimthateachpowerofE +H −E intheseparticu-
el n a
ThestartingpointistheSchr¨odingerequation(1)with lar matrix elements is at least of the order O( m /µ ),
e n
the HamiltonianH, the wavefunction φ, and the energy which we demonstrate in next sections for the leading
p
E decomposed into adiabatic and nonadiabatic parts terms δ(2)E and δ(3)E of the nonadiabatic perturba-
na na
tive expansion.
(H −E )+(E +H −E −δE )]|φ (χ+δχ)+δ′φ i=0.
el el el n a na el na
(18)
(cid:2)
One rewrites this equation to the form
A. Second-order nonadiabatic corrections
(E −H )|δ′φ i = (E +H −E −δE )
el el na el n a na
Intheleadingorderofperturbativetreatmentofnona-
|φ (χ+δχ)+δ′φ i (19)
el na diabatic effects one has
and,sinceδ′φna isorthogonaltoφel,Eq.(14),theformal 1
solution |δ′φnai = (E −H )′ Hn|φelχi, (26)
el el
1 1
|δ′φnai = (Eel−Hel)′ Hn|φel(χ+δχ)i |δχi = (Ea−Eel−Ea−Hn)′ φel Hn δ′φna el,
+(Eel+Hn−(cid:2) Ea−δEna)|δ′φnai , (20) (cid:10) (cid:12)(cid:12) (cid:12)(cid:12) (cid:11) (27)
1
is obtained,wherethe primein the denominatord(cid:3)enotes δ(2)E = φ χ H H φ χ . (28)
na el n (E −H )′ n el
subtraction of the reference state from the Hamiltonian (cid:28) (cid:12) el el (cid:12) (cid:29)
inversion. When δχ andδ′φna onthe righthandside are (cid:12)(cid:12) (cid:12)(cid:12)
The general formula(cid:12)(28), following [16(cid:12)], can be read-
neglected, Eq. (20) becomes the leading nonadiabatic
ilyrearrangedtoamorepracticalform. Fromnowonwe
correction to the wave function. In the next step one
considerthehomonucleardiatomictwo-electronmolecule
takes Eq. (18) and multiplies it from the left by hφ |
el
andseparateoutelectronicmatrixelementsfromthenu-
hφ |E +H −E −δE |φ (χ+δχ)+δ′φ i =0. (21) clear ones
el el n a na el na el
1
Since χ satisfies Eq. (9) the above can be simplified to δ(2)E = d3R χ⋆χ H φ H φ
na n el (E −H )′ n el
Z (cid:20) (cid:28) (cid:12) el el (cid:12) (cid:29)el
= δ(EEenla+|χE+a+δχHin−−hφEeal|)H|δnχ|δi′φnaiel (22) −χ⋆µ∇nRiχ(cid:28)Hnφel(cid:12)(E(cid:12)(cid:12)(cid:12)el−1Hel)′(cid:12)∇(cid:12)(cid:12)(cid:12)Riφel(cid:29)el (29)
and due to Eq. (15) the solution is −∇Riµχ⋆ χ ∇Ri φel(cid:12)(cid:12)(cid:12)(E −1H )′(cid:12)(cid:12)(cid:12)Hnφel
1 n (cid:28) (cid:12)(cid:12) el el (cid:12)(cid:12) (cid:29)el
|δχi = (Ea−Eel−Ea−Hn)′ + ∇Riχµ⋆2∇Rjχ ∇Ri(cid:12)(cid:12)φel (E −1H(cid:12)(cid:12))′ ∇Rjφel
hφ |H |δ′φ i −δE |χ+δχi . (23) n (cid:28) (cid:12) el el (cid:12) (cid:29)el#
el n na el na (cid:12) (cid:12)
(cid:12) (cid:12)
Inthelaststep,(cid:0)onetakesEq.(22),multipliesit(cid:1)fromthe ≡ d3R χ⋆χU(R)−∇(cid:12)Ri χ⋆χ Vi((cid:12)R)
left by hχ|, and obtains Z h +∇i χ⋆ ∇j χ(cid:2)Wij(R(cid:3)) , (30)
R R
δEna =hφelχ|Hn|δ′φnai. (24) i
where the last equation is the definition of potentials
The set of recursive equations (20), (23), and (24) forms U,Vi and Wij. For the Σ electronic state φ
el
the perturbative expansion of the wave functions δ′φ ,
na
δχandenergyδEna. Forexample,startingfrom(24)one Vi = niV, (31)
gets Wij = ninjW +(δij −ninj)W , (32)
k ⊥
1
δEna = hφelχ|Hn(E −H )′ Hn|φel(χ+δχ)i where ~n=R~/R, hence
el el
h
+(E +H −E −δE )|δ′φ i , (25)
el n a na na δ(2)E = d3R χ⋆χδE (R)
na na
i
whichis the sum ofthe leading,Eq.(28), andthe higher Z (cid:26)
+ninj∇i χ⋆ ∇j χW (R)
order nonadiabatic correction, Eq. (44). This perturba- R R k
tive expansion in general assumes that E + H − E
el n a +(δij −ninj)∇i χ⋆ ∇j χW (R) . (33)
is small with respect to the electronic excitation energy. R R ⊥
(cid:27)
4
The function B. Third-order nonadiabatic correction
2 ∂
δE (R)=U(R)+ + V(R) (34) The third order nonadiabatic correction of Eq. (25) is
na
R ∂R
(cid:18) (cid:19)
1
δ(3)E = φ χ H (H +E −E )
is the nonadiabatic correction to the adiabatic energy na el n (E −H )′ n el a
(cid:28) (cid:12) el el
curve Eel(R)+Ea(R) and pseudopotentials U,V,W are: (cid:12)(cid:12)1
×(E −(cid:12) H )′ Hn φelχ
1 el el (cid:12) (cid:29)
(cid:12)
U(R) = (cid:28)Hnφel(cid:12)(cid:12)(Eel−Hel)′(cid:12)(cid:12)Hnφel(cid:29)el, (35) + φelχ Hn (E −(cid:12)(cid:12)1H )′ Hn φelδχ , (44)
V(R) = µ1n (cid:28)Hn(cid:12)(cid:12)φel(cid:12)(Eel−1H(cid:12)(cid:12)el)′(cid:12)~n·∇~Rφel(cid:29)el, (36) where δχ is given(cid:28)in Eq(cid:12)(cid:12)(cid:12)(cid:12). (27).elLet uesl split(cid:12)(cid:12)(cid:12)(cid:12)this su(cid:29)m into
Wk(R) = µ12 ~n·∇~Rφ(cid:12)(cid:12)(cid:12)el (E −1H(cid:12)(cid:12)(cid:12))′ ~n·∇~Rφel , two parts accordingly
n (cid:28) (cid:12)(cid:12) el el (cid:12)(cid:12) (cid:29)e(l37) δ(3)Ena =δ(3)En′a+δ(3)En′′a. (45)
(cid:12) (cid:12)
1 (δij −ninj(cid:12)) (cid:12) While the first part δ(3)En′a involves many terms which
W (R) = are negligible, since they include the third power of µ
⊥ µ2 2 n
n in the denominator, the dominating O(µ−2) term is
n
1
× ∇iφ ∇jφ . (38)
(cid:28) R el(cid:12)(cid:12)(cid:12)(Eel−Hel)′(cid:12)(cid:12)(cid:12) R el(cid:29)el δ(3)En′a = µ12n (cid:28)∇iRφel∇iRχ(cid:12)(Eel−1Hel)′(Hn+Eel−Ea)
In order to simplify th(cid:12)e nonadiabat(cid:12)ic correction of Eq. 1 (cid:12)
(33), one notes that the nuclear wave function χ has a ×(E −H )′ ∇jR(cid:12)(cid:12)φel∇jRχ +O(µ−n3) (46)
definite angular momentum el el (cid:12)(cid:12) (cid:29)
≈ − d3R~n·∇~(cid:12) (χ∗χ)δV(R)+O(µ−3), (47)
(cid:12)R n
χ=χ (R~)=χ (R)Y (~n), (39) Z
Jm J Jm
where
where YJm are spherical harmonics, thus the nonadia- 1 ∂Eel
batic correctioncan be transformed to the form δV(R)= (48)
2µ2 ∂R
n
1
δ(2)Ena = R2dR χ′J2Wk(R) (40) ×(cid:28)~n·∇~Rφel(cid:12)[(Eel−Hel)′]2(cid:12)~n·∇~Rφel(cid:29)el ,
Z (cid:26) (cid:12) (cid:12)
+χ2 δE (R)+ J(J +1)W (R) . and this correction is in(cid:12)(cid:12)cluded into δ(cid:12)(cid:12)Ena of Eq. (34),
J na R2 ⊥ which becomes now
(cid:20) (cid:21)(cid:27)
2 ∂
Let us note, that Eq. (40) can also be expressed in δE (R)=U(R)+ + [V(R)+δV(R)]. (49)
na
R ∂R
terms of an expectation value of an effective nonadia- (cid:18) (cid:19)
batic Hamiltonian δHna The second term δ(3)E′′ can be obtained from the non-
na
perturbative solution of the nuclear equation with the
δ(2)Ena =hχJ|δHna|χJi, (41) nonadiabatic Hamiltonian δHna from Eq. (42). Namely,
for the states with the rotational quantum number J,
where δ(3)E′′ takes the form
na
1 ∂ ∂ J(J +1) δ(3)E′′ = hχ |δH |δχ i (50)
δH = − R2W (R) + W (R) na J na J
na R2 ∂R k ∂R R2 ⊥
1
+δE (R). (42) = χ δH δH χ ,
na J na ′ na J
* (cid:12) Ea−Eel−Ea−Hn (cid:12) +
(cid:12) (cid:12)
(cid:12) (cid:12)
Now, the nonadiabatic correction δχ of Eq. (27) can be which is interpre(cid:12)ted as(cid:0)a second order corre(cid:1)ction d(cid:12)ue to
conveniently rewritten in terms of δH (cid:12) (cid:12)
na δH .
na
1
|δχ i= δH |χ i, (43)
J E −E −E −H ′ na J C. Nuclear equation with effective potentials
a el a n
(cid:0) (cid:1)
where it is understood, that the derivatives with respect Instead of calculating the second order correction
to electronic variables of the function χ do vanish. Eq. (41) and the third order correction of Eq. (50), it
5
is more appropriate to solve nonperturbatively the cor- E. Evaluation of the wave function derivatives
responding equation
TheelectronicmatrixelementsinEq.(29)involvemul-
H +E (R)+E (R)+δH χ =Eχ , (51)
n el a na J J
tiple differentiation of the electronic wave function with
where(cid:2) the total energy E is (cid:3) respect to the internuclear distance R, which is difficult
to calculate directly. Therefore, following Ref. [16], we
E =E +δ(2)E +δ(3)E +... (52)
a na na rewrite these terms to a more convenient form, where
differentiationistakenoftheCoulombpotential,namely
andthisisdoneinthiswork. Moreover,fromcomparison
of subsequent terms of Eq. (42) with those of Eq. (11),
1
one concludes that the first term can be interpreted as ∇~ φ = ∇~ (V)φ , (59)
R el (E −H )′ R el
an R-dependent correction to the nuclear reduced mass el el
µ , 1
n ∇2φ = ∇2(V)φ +2∇~ (V −E )(60)
R el (E −H )′ R el R el
1 1 el el (cid:26)
≡ +W (R), (53)
2µk(R) 2µn k (E −1H )′ ∇~R(V)φel +φelhφel|∇R2|φeliel,
whereas the second term—as an R-dependent correction el el (cid:27)
to the inverse of the moment of inertia µ⊥(R)R2 The derivatives of the potential V
1 1
≡ +W (R). (54) 1 1 1 1 1 1
2µ (R) 2µ ⊥ V =− − − − + + , (61)
⊥ n r r r r r R
1A 1B 2A 2B 12 AB
With these newly defined functions one can write the
are the following
radial equation as
1 ∂ R2 ∂ J(J +1) 1 ~r ~r ~r ~r R~
(cid:20)−2R2 ∂R µk(R) ∂R + 2µ⊥(R)R2 +Y(R)(cid:21)χJ(R) ∇~R(V) = 2(cid:18)−r113AA + r113BB − r223AA + r223BB(cid:19)− R3 ,
=EχJ(R), (55) (62)
∇2(V) = π δ(~r )+δ(~r )+δ(~r )+δ(~r )
where R 1A 1B 2A 2B
−(cid:2)4πδ(R~), (cid:3)(63)
Y(R)=E (R)+E (R)+δE (R) (56)
el a na
forms the effective nonadiabatic potential. andthe matrixelementswiththese operatorsarereadily
evaluated. The presence of the Dirac delta operators in
Eq. (60) may potentially decrease the accuracy of the
D. Asymptotics of the effective masses evaluationofthosequantitieswhichcontain∇2φ . Ifwe
R el
note that
The adiabatic correction E (R), Eq. (10), and the
a
nonadiabaticcorrectionδEna(R), Eq.(34), do notvanish ∇2el(V)=π[δ(~r1A)+δ(~r1B)+δ(~r2A)+δ(~r2B)], (64)
atlargeinternucleardistances. Forexample,forthelarge
atomic separation in the hydrogen molecule, E (R) and wecangetridofDiracdeltas bya simple rearrangement
a
δE (R) are equal to m /m and −(m /m )2, respec- of the nuclear part of the Hamiltonian Hn to the form
na e p e p
tively,whichcorrespondsto the firstterms inthe expan-
1 1
sion of the atomic reduced mass µ = (1/mp+1/me)−1 Hn =−2µ ∇R2 −∇2el − µ ∇2el. (65)
in the electron-nucleus mass ratio n n
(cid:0) (cid:1)
µ m /m m m 2 m 3 The difference in parenthesis collects the terms of
e p e e e
1− = = − + −...
m 1+m /m m m m Eq. (63) and (64) which cancel out, up to the negligi-
e e p p (cid:18) p(cid:19) (cid:18) p(cid:19) (57) ble −4πδ(R~) term, so that we can write it down as
LargeRasymptoticsofthe pseudopotentialsW (R)and
k
W⊥(R) are equal to −me/m2p, which is related to the ∇R2 −∇2el φel = (66)
changeinEqs.(53)and(54)ofthe reducednuclearmass 2 1
(cid:0) (cid:1)∇~ (V −E ) ∇~ (V)φ
µn to the reducedmassµA =(mp+me)/2 oftwohydro- (Eel−Hel)′ R el (Eel−Hel)′ R el
gen atoms :
2 1
− ∇~ (V) ∇~ (V)φ +λφ
1 1 1 (E −H )′ el E −H el el el
= = el el el el
2µ (∞) 2µ (∞) 2µ
k ⊥ A
with some constant λ. In this way the Dirac delta terms
1 1 m
= = 1− e +··· . (58) are eliminated at the expense of employing additional
mp+me mp(cid:18) mp (cid:19) basis sets for evaluation of the last resolventin Eq.(66).
6
IV. NUMERICAL CALCULATIONS parameters were optimized variationally with respect to
theclampednucleienergywiththetargetaccuracyofthe
Inordertoformtheradialnonadiabaticequation(55), order of a fraction of microhartree. The bases B, C, and
apartfromtheclampednucleienergyE (R)andtheadi- D were intended for evaluation of the resolvents present
el
abaticcorrectionEa(R),itisnecessarytoknowthepseu- inEqs.(59)and(66). ThetwoΣ+g bases(BandD)were
dopotentials U(R),V(R)+ δV(R),W (R), and W (R). optimized in the presence of the basis A: the first 600
k ⊥
The evaluation of these functions is the main numeri- terms were taken from φel wave function and their non-
cal task described here. Calculations were performed at linear parameters were kept fixed during the optimiza-
80 points including long and very short internuclear dis- tion, only the remaining 600 terms were actually opti-
tances R. At each point several electronic wave func- mized. This ensures that the internal wave function φel
tions of different symmetry were generated. All these is well represented at every step of optimization. Then,
functions were represented as linear expansions of prop- the subtraction of the reference state, denoted by the ′
erlysymmetrizedtwo-electronbasisfunctions. Thebasis symbolwithintheresolvent,wasachievedbyorthogonal-
functions were taken in the form of exponentially corre- ization of ∇~R(V)|φeli to the internal |φeli. In the final
lated Gaussians (ECG) calculations the three bases A, B, and D were assembled
together to form a 1800-term Σ+ basis applied not only
g
ψ (~r ,~r ) = (1+Pˆ )(1±ˆı)Ξ (67) to evaluate the pertinent resolvents but also to expand
k 1 2 12 k
theexternalgroundstatefunctionφ . Thetwoungerade
2 el
bases (E and F) were employedto evaluate the resolvent
× exp − A (~r −~s )(~r −~s ) ,
ij,k i i,k j j,k
and to form the components of the scalarproduct in the
i,j=1
X second term of the right hand side of Eq. (66).
where the matrices A and vectors ~s contain nonlin- The adiabatic potential of the nuclear Schr¨odinger
k k
ear parameters, 5 per basis function, to be variationally equation(12)wascomposedoftheclampednucleienergy,
optimized. The antisymmetry projector (1 + Pˆ12) en- Eel(R), and the adiabatic correction Ea(R). For Eel(R)
suressingletsymmetry,the spatialprojector(1±ˆı)—the we used the analytic potential constructed by L ach [18]
gerade (+) or ungerade (−) symmetry, and the Ξ pref- on the basis of the energy points computed by Cencek
k
actorenforcesΣstateswhenequalto1,orΠstateswhen from 1200-term ECG wave functions [19] and Sims and
equal to y —the perpendicular Cartesian component of Hagstrom from Hylleraas wave functions [20]. Their en-
i
the electron coordinate. ergy points where converged up to 13 significant digits.
Six different basis sets have been prepared to repre- The adiabatic correction Ea(R) was evaluated as an ex-
sent the variety of the electronic wave functions. To en- pectation value of the Hamiltonian Hn, Eq. (10),
sure high accuracy of the potentials, the basis sets have
1
been variationally optimized with respect to pertinent E (R) = − hφ |∇2 +∇2|φ i . (68)
a 2µ el R el el el
goal functions according to the specification in Table I. n
Toavoidthecumbersomedifferentiationoftheelectronic
wave function with respect to the internuclear distance
TABLE I: Goal functions used in optimization of the basis
wereplacedtheexpectationvalueinthefirsttermonthe
sets.
right hand side of Eq. (68) by an equivalent expression
Label Size Symmetry Goal function
hφ |∇2|φ i =−h∇~ φ |∇~ φ i , (69)
el R el el R el R el el
A 600 Σg Eel
which, with the help of Eq. (59), can be further trans-
1
B 600† Σg D∇kR(V)(Eel−Hel)′∇kR(V)E formed to
1
C 1200 Πg D∇⊥R(V)Eel−Hel∇⊥R(V)E hφel|∇R2|φeliel =− ∇~R(V) 1 2∇~R(V) .
1 * (E −H )′ +
D 600† Σg D∇2el(Eel−Hel)′∇2elE el el e(l70)
(cid:2) (cid:3)
1 The formula (70), when evaluated with the optimized
E 600 Σu D∇kel(V)Eel−Hel∇kel(V)E bases A, B, and C, yields the adiabatic correction with
F 600 Πu D∇⊥el(V)Eel−1Hel∇⊥el(V)E acunrvaeccwuarasctyheonf aotbtlaeainsetd1bpypmme.anTshoefa1d0i-apboaitnitcppioetceenwtiisael
† Optimized along with the fixedbasis A. polynomial interpolation.
TheelectronicmatrixelementsU,V+δV,W ,W en-
k ⊥
Particular goal functions have been chosen to reflect tering Eq. (33) were evaluated with the ECG basis sets
the contents of the expressionthe basis set is to be used describedabove,yieldingsmoothfunctionsofR. Because
for. The first basis (labeled A), composed of 600 ECG for the highest vibrational levels the nuclear wave func-
functions (67), was employed to expand the X1Σ+ elec- tions are spread out and the contributions from larger
g
tronic ground state wave function φ . Their nonlinear internuclear distances are non-negligible, the functions
el
7
U(R),V(R),andW(R)wererepresentedbytheirasymp- The resulting radial equation
totic forms:
∂ 1 ∂
η (R)
J
∂R µ (R) ∂R
k
U(R) ≈ u +u /R6+u /R8,
0 6 8 W′(R) J(J +1)
V(R) ≈ v9/R9+v11/R11, =−2 E−Y(R)− kR − 2µ (R)R2 ηJ(R), (73)
δV(R) ≈ v /R7+v /R9, (cid:20) ⊥ (cid:21)
7 0
Wk(R) ≈ wk0+wk12/R12+wk14/R14, is solved numerically for the function ηJ(R) = RχJ(R).
We used the code developed by W. Johnson, described
W (R) ≈ w +w /R12+w /R14, (71)
⊥ ⊥0 ⊥12 ⊥14 recently in his book [21], and modified it to account
for the dependence of the mass on the internuclear dis-
subject to u = w = w = −(m /m )2 restriction tance. In the calculations we used the following con-
0 k0 ⊥0 e p
stants [22]: the proton mass m = 1836.15267247 m
(inatomic units). The remaining,free parametersu ,v , p e
i i
and the energy units conversion factor 1 hartree =
and w were determined by fitting the above functions
i
219474.6313705cm−1.
to the calculated points in the range of h6.0,10.0ibohrs.
Because at distances R > 6, the numerical precision of
the potentials U and V was not high enough, we used
V. RESULTS AND DISCUSSION
lowerR-valuesfortheextrapolation. AttheoriginR=0
all the potentials are finite with V ∼ R, W ∼ R2, and
k
W ∼ R2. Numerical results for δE , W , and W are Inamolecule,themovingnucleiare’coated’withelec-
⊥ na k ⊥
shown graphically in Fig. 1. trons and the amount of additional mass carried by the
nuclei changes with R. For a homonuclear molecule in
a purely vibrational state, the effective mass of the nu-
cleus at given R, m (R), is just twice the reduced mass
x 10-7 k
µ (R) defined in Eq. (53). Analogously, for a rigid ro-
k
0 tating molecule, the effective nuclear mass m (R) is re-
⊥
lated to µ (R) of Eq. (54). Thus, their R-dependence
⊥
-1
can be determined explicitly from the potentials W (R)
ε k
-2 δ na(R) W⊥(R) ianndthWe t⊥w(oRe),ffercetsipveectniuvcelleya.rFmiga.ss2esillwuistthratthees itnhteercnhuacnlegaers
-1/m2
p distance in H . The functions m (R) and m (R) join
-3 2 k ⊥
smoothly the proton mass at the united atom limit with
-4 the hydrogen atomic mass (mp +me) at the separated
atoms limit. Interestingly, for R ≥ 2.41 a.u., the effec-
-5 WII(R) tive mass mk(R) is greater than the sum of proton and
electron masses, reaching m +1.6m at the maximum
p e
-6
located near R=3.8 a.u.
0 2 4 6 8 10
The radialequation(73)hasbeensolvedforallbound
R/a.u.
states with three versions of the potential Y(R):
FIG. 1: (Color online) The nonadiabatic potentials (in a.u.):
δEna(R) (bolded, black), Wk(R) (lower dotted, blue), and Y(R)=Eel(R), BO
W⊥(R) (upper dotted, red). The potentials asymptote goes Y(R)=E (R)+E (R), adiabatic
at −1/m2 =−2.966077×10−7. el a
p Y(R)=E (R)+E (R)+δE (R), nonadiabatic
el a na
yielding three sets of dissociation energies. The corre-
The nonadiabatic correction to energy levels can be sponding dissociation thresholds were −1 hartree in the
computedperturbativelyfromEq.(40)ashasbeendone BOapproximation,−1+m /m hartreeinadiabaticap-
e p
in our previous paper [16], or, in a more accurate way, proximation, and −1 + m /m − (m /m )2 hartree in
e p e p
theradialequation(51)canbesolvedforthe totalnona- the nonadiabatic level of theory. The results are listed
diabaticenergy. Thissecondmethodisdescribedshortly in Table II, where for each pair of quantum numbers v
below. The first term including the second order deriva- and J three entries are given(in cm−1): the BO dissoci-
tive is rewritten to the more convenient form ation energy, the adiabatic correction, and the nonadia-
batic correctionto the dissociationenergy. Thus,the to-
talnonrelativisticdissociationenergycanbeobtainedby
1 ∂ 1 ∂
− R2 +W (R) summing up all three entries. The only exception is the
R2 ∂R 2µ k ∂R
(cid:18) n (cid:19) state with v = 14 and J = 4, for which a nonadiabatic
1 ∂ 1 ∂ Wk′(R) level lying just beneath the dissociation threshold has
=− +W (R) R+ . (72)
R ∂R 2µ k ∂R R beenpredictedalthoughneitherBOnoradiabaticbound
(cid:18) n (cid:19)
8
states exist. The entry given for this state is the energy ods [21]. Although representation of the nonadiabatic
separationfrom the nonadiabatic dissociation threshold. correction by pseudopotentials has been studied previ-
ously(see[9,11,12]andreferencestherein),wehaverig-
orously derived new expressions for W , W , and δE
k ⊥ na
1838.0 functions of the internuclear distance, which give nona-
m (R) diabaticcorrectionswith O(m /µ ) accuracy. Moreover,
|| e n
wehaveperformedexplicitnumericalcalculationsforthe
1837.5 simplestexampleofH moleculeinordertoverifytheob-
2
m +m tained perturbative formulae. All the electronic matrix
p e
elements with differentiation of the wave function over
1837.0
the internuclear distance were rewritten in a convenient
form involving differentiation of the Coulomb potential.
m⊥(R) Suchanapproachenablesachievingmuchhighernumer-
1836.5
ical precision even for the well known adiabatic correc-
m tion. The final accuracy of all three nonadiabatic func-
p
tionsislimitedonlyby theneglectedhigherorderterms,
1836.0
0 2 4 6 8 10 namely O(me/µn) relative to the leading order, and re-
R/a.u. sultsforrovibrationallevelsagreewithinthisuncertainty
with former results of Wolniewicz [15]. Having accurate
FIG.2: (Color online) TheR-dependenceoftheeffectivenu-
nonrelativistic energies one can include relativistic and
clear masses m (R) and m (R) (in a.u.).
k ⊥ QEDcorrections,whichbecomesignificantforthe states
close to the dissociation threshold. Relativistic correc-
Our nonadiabatic corrections agree very well with tions for the electronic ground state of H are known
2
those computed by Wolniewicz [15] for rovibrational accuratelyforawiderangeofinternucleardistances[23].
states of J ≤ 10. For all rotational states of the lowest TheirlargeRasymptotics,includingQEDcorrections,is
vibrational level the difference is merely 0.0002 cm−1 or presently investigated by the Jeziorski group [24]. Com-
less. In all the cases the difference is smaller than 0.1%, biningalltheknowledgetogetherwouldenableachieving
which corresponds to the magnitude of the neglected at least an order of magnitude increase in the precision
higher order terms of relative order O(1/µn). This is of theoretical prediction for all molecular states of H2.
the first confirmation of the Wolniewicz’s results [15] for
Analogous calculations can be performed for D and
J >0 rovibrationalstates and, simultaneously,a numer- 2
HDmolecules. Itisworthnoting,thatintheformercase
ical validation of the nonadiabatic perturbation theory
there is no need to recompute the pseudopotentials and
presented in this work. Let us note, that the nonadia-
only a proper rescaling due to different reduced masses
baticcorrectionstothedissociationenergyfromTableII
is required. The calculations on the latter system, due
differinsignfromthecorrectionstotheadiabaticenergy
to difference in mass of the proton and deuteron, would
of [15] and their absolute values differ by the constant
involve additional correction to δE coming from the
1/m2 = 0.06509787 cm−1 corresponding to the asymp- na
p last term of Eq. (5). While this perturbative approach
totic value of the nonadiabatic potential δE (R), dis-
na can be further extended to larger diatomic molecules, it
cussed in Subsec. IIID.
would be more challenging to investigate three-atomic
molecules such as H+—a system of great astrophysical
3
interest. We expect no principal difficulties in such cal-
VI. SUMMARY
culations except for much more increased demands in
computer resources needed to perform optimization of
We havepresentedthe nonadiabatic perturbationthe-
pertinent wave functions.
ory applicable to any molecule in an arbitrary rovibra-
tional state. The leading nonadiabatic corrections for
the diatomic molecule are expressed in terms of three
R-dependent functions depicted in Fig. 1: the nuclear
Acknowledgments
reduced mass in Eq. (53), the moment of inertia in
Eq. (54), and the correction Eq. (49) to the adiabatic
potential, which enter the radial Schr¨odinger equation We are indebted to L. Wolniewicz for valuable com-
for nuclei Eq. (51). This equation can be solved for an ments. Part of the computations has been performed in
arbitrary molecular states by standard numerical meth- Poznan´ Supercomputing and Networking Center.
[1] W. Kol os and L. Wolniewicz, Rev. Mod. Phys. 35, 473 [2] S. Bubin and L. Adamowicz, J. Chem. Phys. 118, 3079
(1963). (2003).
9
[3] M. Cafiero and L. Adamowicz, J. Chem. Phys. 122,
184305 (2005).
[4] S.Bubin,L.Adamowicz,andM.Molski,J.Chem.Phys.
123, 134310 (2005).
[5] M. Stanke, D. K¸edziera, S. Bubin, M. Molski, and
L. Adamowicz, Phys. Rev.A 76, 052506 (2008).
[6] M. Born and J. R. Oppenheimer, Ann. Phys. 84, 457
(1927).
[7] M. Born and K. Huang, Dynamical Theory of Crystal
Lattices, Oxford University Press (1955) Appendix8.
[8] R.M.HermanandA.Asgharian,J.Mol.Spectr.19,305
(1966).
[9] P. R. Bunker and R. E. Moss, Molec. Phys. 33, 417
(1977), ibid 33, 425 (1977).
[10] R.M.HermanandJ.F.Ogilvie.Adv.Chem.Phys.103,
187 (1998).
[11] W. Kutzelnigg, Mol. Phys. 105, 2627 (2007).
[12] R. Jaquet and W. Kutzelnigg, Chem. Phys. 346, 69
(2008).
[13] D.W. Schwenke,J. Phys.Chem. A 105, 2352 (2001).
[14] P.R.BunkerandR.E.Moss, J.Mol.Spectrosc.80,217
(1980).
[15] L. Wolniewicz, J. Chem. Phys. 103, 1792 (1995).
[16] K.PachuckiandJ.Komasa,J.Chem.Phys.129,034102
(2008).
[17] M. Stanke, D. K¸edziera, S. Bubin, M. Molski, and
L. Adamowicz, J. Chem. Phys. 128, 114313 (2008).
[18] G.L ach, PhD thesis, Universityof Warsaw (2008).
[19] W. Cencek, private communication.
[20] J. S. Sims and S. A. Hagstrom, J. Chem. Phys. 124,
094101 (2006).
[21] W. Johnson, Atomic Structure Theory (Springer, 2007).
[22] P. J. Mohr and B. N. Taylor, Rev. Mod. Phys. 77, 1
(2005).
[23] L. Wolniewicz, J. Chem. Phys. 99, 1851 (1993).
[24] K. Piszczatowski, G. L ach, and B. Jeziorski, Phys. Rev.
A 77, 062514 (2008).
10
TABLEII:Thedissociation energy of therovibrational states ofH2 (in cm−1). Foreach pair ofvibrational (v) androtational
(J) quantum numbers, three entries are given: the BO dissociation energy, the adiabatic correction, and the nonadiabatic
correction. The sum of thethreenumbersgives the total nonrelativistic dissociation energy of the(v, J) state.
v\J 0 1 2 3 4 5 6 7
0 36112.5927 35994.0372 35758.0143 35406.6660 34943.1216 34371.4040 33696.3152 32923.3094
0 5.7711 5.8348 5.9611 6.1481 6.3931 6.6924 7.0418 7.4367
0 0.4339 0.4406 0.4539 0.4739 0.5005 0.5338 0.5736 0.6200
1 31949.1892 31836.5534 31612.3337 31278.6001 30838.3755 30295.5437 29654.7385 28921.2187
1 7.1740 7.2311 7.3443 7.5119 7.7310 7.9983 8.3096 8.6605
1 1.2704 1.2761 1.2874 1.3044 1.3271 1.3555 1.3896 1.4294
2 28021.4345 27914.5889 27701.9154 27385.4159 26968.0130 26453.4605 25846.2343 25151.4125
2 8.3336 8.3840 8.4839 8.6314 8.8240 9.0583 9.3304 9.6358
2 2.0271 2.0318 2.0414 2.0557 2.0749 2.0989 2.1279 2.1618
3 24324.4498 24223.3066 24022.0055 23722.4851 23327.5744 22840.9061 22266.8106 21610.1998
3 9.2420 9.2855 9.3714 9.4982 9.6632 9.8633 10.0945 10.3525
3 2.7088 2.7127 2.7206 2.7326 2.7485 2.7686 2.7929 2.8214
4 20855.2072 20759.7245 20569.7150 20287.0583 19914.4976 19455.5548 18914.4286 18295.8815
4 9.8890 9.9251 9.9966 10.1017 10.2378 10.4019 10.5902 10.7985
4 3.3187 3.3219 3.3283 3.3379 3.3508 3.3670 3.3868 3.4100
5 17612.7145 17522.9036 17344.2107 17078.4626 16728.3253 16297.2225 15789.2368 15209.0018
5 10.2629 10.2915 10.3476 10.4298 10.5356 10.6618 10.8049 10.9606
5 3.8576 3.8599 3.8647 3.8719 3.8816 3.8938 3.9087 3.9262
6 14598.2891 14514.2237 14346.9982 14098.3935 13771.0098 13368.1874 12893.9129 12352.7160
6 10.3525 10.3729 10.4130 10.4712 10.5449 10.6313 10.7268 10.8271
6 4.3218 4.3232 4.3262 4.3306 4.3365 4.3440 4.3531 4.3639
7 11815.9479 11737.7779 11582.3242 11351.3310 11047.3466 10673.6480 10234.1517 9733.3170
7 10.1486 10.1606 10.1839 10.2169 10.2573 10.3020 10.3476 10.3901
7 4.7029 4.7031 4.7037 4.7044 4.7055 4.7067 4.7081 4.7095
8 9272.9561 9200.9278 9057.7445 8845.1263 8565.5890 8222.3732 7819.3619 7360.9906
8 9.6483 9.6517 9.6577 9.6648 9.6707 9.6726 9.6670 9.6500
8 4.9838 4.9824 4.9796 4.9752 4.9692 4.9614 4.9516 4.9393
9 6980.5984 6915.0837 6784.9238 6591.8323 6338.3203 6027.6323 5663.6708 5250.9187
9 8.8590 8.8537 8.8422 8.8232 8.7945 8.7533 8.6963 8.6196
9 5.1364 5.1325 5.1246 5.1124 5.0958 5.0743 5.0472 5.0138
10 4955.2699 4896.8118 4780.7739 4608.8912 4383.7137 4108.5514 3787.4126 3424.9438
10 7.8021 7.7882 7.7597 7.7151 7.6525 7.5691 7.4616 7.3260
10 5.1160 5.1082 5.0925 5.0684 5.0354 4.9926 4.9386 4.8720
11 3220.0418 3169.4253 3069.1017 2920.8723 2727.3991 2492.1652 2219.4373 1914.2414
11 6.5140 6.4918 6.4469 6.3777 6.2821 6.1573 5.9996 5.8045
11 4.8566 4.8429 4.8151 4.7725 4.7140 4.6379 4.5418 4.4224
12 1806.9489 1765.3260 1683.0605 1562.1049 1405.3733 1216.7368 1001.0464 764.2101
12 5.0372 5.0068 4.9451 4.8505 4.7203 4.5506 4.3360 4.0685
12 4.2657 4.2429 4.1965 4.1254 4.0271 3.8984 3.7345 3.5281
13 760.3903 729.5279 668.9437 580.9383 469.0254 338.0600 194.5461 47.4825
13 3.3933 3.3526 3.2697 3.1417 2.9631 2.7253 2.4123 1.9886
13 3.2221 3.1850 3.1095 2.9927 2.8294 2.6113 2.3232 1.9317
14 141.7951 124.7523 92.3077 48.0033
14 1.5343 1.4739 1.3479 1.1416
14 1.5847 1.5226 1.3933 1.1825 0.0887†
† This state appears as a resonance in BO and adiabatic approximations. Theentry is a dissociation energy of this
nonadiabatic level.
