Table Of Content

Version2.0 Ground state wave function and energy of the lithium atom Mariusz Puchalski∗ and Krzysztof Pachucki† Institute of Theoretical Physics, Warsaw University, Hoz˙a 69, 00-681 Warsaw, Poland 6 0 Abstract 0 2 Highly accurate nonrelativistic ground-state wave function and energy of the lithium atom is obtained n a in the Hylleraas basis set. The leading relativistic corrections, as represented by Breit-Pauli Hamiltonian, J 0 are obtained infair agreement with theformer results. Thecalculational method is based onthe analytical 3 evaluation ofHylleraasintegrals withthehelpofrecursion relations. ] h p - PACSnumbers:31.25.Eb,31.30.Jv31.15.Pf,02.70.-c m o t a . s c i s y h p [ 1 v 3 1 2 1 0 6 0 / s c i s y h p : v i X r a ∗[email protected] †[email protected] 1 I. INTRODUCTION Theoreticalpredictionsfortheenergylevelsoflightfew-electronatomsaremuchlessaccurate thanforthehydrogenicsystems. Itisfortworeasons. Thenonrelativisticwavefunctionhastoin- cludeelectroncorrelationstoahighdegreeofaccuracy. ThiscanbeachievedbyusingaHylleraas basis set, but it is quite difficult to evaluate integrals with Hylleraas functions for three and more electrons. Thesecond reason is thedifficultyin theaccurate treatment ofrelativisticand radiative corrections. The commonly used Dirac-Coulomb Hamiltonian for few-electron atoms does not includerelativisticcorrectionsproperlyasitcannotbederivedfromquantumelectrodynamicthe- ory and its continuous spectrum ranges from to + . One of the possible approaches is the −∞ ∞ derivationofaneffectiveHamiltonian[1]withinthesocalledNRQEDtheory. Matrixelementsof thisHamiltoniangiveexactcorrectiontotheenergyatspecifiedorderinthefinestructureconstant α. However, this Hamiltonian becomes quite complicated at higher orders and for example mα6 corrections has been obtained for few low lying states of helium only [2, 3], not for lithium nor berylliumatoms. Theoretical predictions for light hydrogen-like atoms are at present limited by uncertainty in higher-ordertwo-loopelectronself-energycorrections[4],whichisafewkHzforthe1Sstate. For helium-likeatomspredictionsareapproximately103 timeslessaccurate. Since,thenonrelativistic wave function was computed very accurately using Hylleraas [5] or exponential basis sets [6], the uncertainty in its levels comes mainly from the unknown mα7 terms. These corrections are currently under investigation in the context of helium 23P fine splitting. For lithium atoms, J the Hylleraas functions give very accurate nonrelativistic wave function and energies [7], but the precise calculation of three-electron integrals with Hylleraas functions is very time consuming [8, 9], and so far no result for mα6 corrections have been obtained. For the beryllium atom the most accurate results have been obtained with explicitly correlated Gaussian functions [10]. Althoughitwaspossibletocalculateaccurately theleadingrelativisticandQED corrections[11], the final accuracy is limited by the nonrelativisticenergy. Moreover,this basis cannot be used for higherorder correctionssinceGaussianwavefunctionsdo notfulfillthecusp condition. SofarthemostaccurateresultsforvariousstatesofthelithiumatomwereobtainedbyYanand Drake in Ref. [7]. Here, we present even more accurate results for the lithium ground state, as a demonstrationof an analyticmethod to computethe integralswith Hylleraas functions [12]. This new method is based on recursion relations between integrals with different powers of electron- 2 nucleus and inter-electron distances, which are fast and numerically stable for generating large basissets. Ourresultfor thegroundstateenergy E = 7.4780603239041(+10), (1) − −50 is significantly below the previous one, obtained in [7], which is E = 7.4780603236503(71). − As a further application of the analytic approach, we obtain the leading relativistic corrections to the binding energy by the calculation of the expectation value of Breit-Pauli Hamiltonian in Eq. (13). For this we used recursion relations for extended Hylleraas integrals with 1/r2 and 1/r2 ij i terms. Theyhavebeen derivedin[13]and inthiswork respectively. In the next Section we construct the nonrelativistic wave function, similarly to Ref. [7] and obtain the ground state nonrelativistic energy and the wave function. In Sec. III we compute the leading relativistic correction as given by the Breit-Pauli Hamiltonian. In Sec. IV we derive recursion relations for Hylleraas integrals containing 1/r2 which among others, are necessary i for relativistic matrix elements. In Sec. V we summarize our result and present prospects for calculation of higher order terms as well as the calculation of Hylleraas integrals for 4 and more electrons. II. NONRELATIVISTICWAVEFUNCTIONANDENERGY In the construction of the wave function we closely follow the works of Yan and Drake in [7]. Theground statewavefunction Ψ is expressedas a linearcombinationof ψ, theantisymmetrized productofφ and thespinfunctionχ ψ = [φ(~r ,~r ,~r )χ], (2) 1 2 3 A φ(~r ,~r ,~r ) = e−w1r1−w2r2−w3r3rn1rn2rn3rn4rn5rn6, (3) 1 2 3 23 31 12 1 2 3 χ = α(1)β(2)α(3) β(1)α(2)α(3), (4) − withall n nonnegativeintegersandw R . ThematrixelementoftheHamiltonianH i i + ∈ 3 ~p2 Zα 3 α H = a + , (5) (cid:18) 2 − r (cid:19) r Xa=1 a aX>b=1 ab orofanyspinindependentoperatorcan beexpressedafter eliminatingspinvariables,as ψ H ψ′ = 2φ(1,2,3)+2φ(2,1,3) φ(3,1,2) φ(2,3,1) φ(1,3,2) φ(3,2,1) h | | i h − − − − | H φ′(1,2,3) . (6) | i 3 In this way the calculation of this matrix elements is brought to Hylleraas integrals, namely the integralswithrespect tor oftheform i d3r d3r d3r f(n ,n ,n ,n ,n ,n ) = 1 2 3 e−w1r1−w2r2−w3r3 1 2 3 4 5 6 Z 4π Z 4π Z 4π rn1−1rn2−1rn3−1rn4−1rn5−1rn6−1, (7) 23 31 12 1 2 3 with nonnegative integers n . These are performed analytically for n ,n ,n = 0,1 [14] and by i 1 2 3 recursion relationsforlargern usingformulasderivedin[12]. i Thetotalwavefunctionisgenerated from all φin Eq. (3)withn satisfyingcondition i 6 n Ω, (8) i ≤ Xi=1 forΩbetween3and 12. Foreach Ωweminimizeenergywithrespect tothefreeparameters w in i Eq. (3). Wenoticedthattheuseofonlyonesetofw ’sdoesnotleadtoaccurateresults,therefore, i following Yan and Drake [7], we divide the whole basis set into 5 sectors, each one with its own set ofw ’s. Thisdivisiongoes as follows[7] i sector1: alln , n = 0, n = 0; 3 1 2 sector2: alln , n = 0, n = 0; 3 1 2 6 sector3: alln , n = 0, n = 0; 3 1 2 6 sector4: n = 0, n = 0, n = 0; 3 1 2 6 6 sector5: n = 0, n = 0, n = 0; 3 1 2 6 6 6 To avoid numerical instabilities, within each sector we drop the terms with n > n (or n < n ) 4 5 4 5 and for n = n drop terms with n > n (or n < n ). This division allows for a significant 4 5 1 2 1 2 improvements of nonrelativisticenergies by optimization of all five sets of w ’s. These nonlinear i parameters are obtainedbyNewtonmethodofsearchingzeros usinganalyticderivatives ∂E ∂Ψ ∂Ψ = 2 Ψ H 2E Ψ . (9) ∂w (cid:28) (cid:12) (cid:12)∂w(cid:29)− (cid:28) (cid:12)∂w(cid:29) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) In the numerical calculations, we use(cid:12)sex(cid:12)tuple precision fo(cid:12)r recursion relations and quadruple precision for all other arithmetics to obtain the wave function and the energy up to Ω = 12. The resultsobtained forgroundstateenergies are presentedin TableI. Thepenultimaterow is aresult of extrapolation to infinite length of the basis set, and the last raw are previous results of Yan and Drake[7]. Theresultforthenonrelativisticenergyissignificantlybelowthepreviousestimate[7] andindicatesthatextrapolationtoinfinitebasislengthdoesnotalwaysgivetherightresult. Inthe 4 TABLEI:Groundstatenonrelativistic energies andexpectation valuesofDiracδ-functions obtained using Drachmanformulae[16]forvariousbasislength. Ω No. ofterms E(Ω) δ3(r ) δ3(r ) a a a>b ab P P 3 50 -7.4779815240897 13.84344680398 0.54416435192 4 120 -7.4780523346422 13.84228864167 0.54433156416 5 256 -7.4780594639158 13.84250917463 0.54432787045 6 512 -7.4780602086637 13.84263796667 0.54432526063 7 918 -7.4780603103629 13.84260666238 0.54432478885 8 1589 -7.4780603205076 13.84260824076 0.54432469702 9 2625 -7.4780603234501 13.84261009857 0.54432462945 10 4172 -7.4780603237750 13.84261069867 0.54432462757 11 6412 -7.4780603238610 13.84261077919 0.54432463150 12 9576 -7.4780603238897 13.84261078106 0.54432463205 -7.4780603239041(+10) 13.84261078346(100) 0.54432463396(50) ∞ ∞ −50 Refs. [7,15] -7.4780603236503(71) 13.842609642(55) 0.54432979(31) ∞ sameTablewepresentresultsfortheDiracδ functions,whichalsodiffersfrompreviousresultsin [15]. Weobserve, thethenumberofsignificant digitsforDiracδ is increased by usingDrachman formulae[16], namely 1 1 4π Ψ δ3(r ) Ψ = 2 Ψ (E V) Ψ ~ Ψ ~ Ψ , (10) ab Ψ c c (cid:28) (cid:12)r − (cid:12) (cid:29)− (cid:28)∇ (cid:12)r (cid:12)∇ (cid:29) (cid:10) (cid:12) (cid:12) (cid:11) (cid:12) ab (cid:12) Xc (cid:12) ab(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 1 (cid:12) (cid:12) 1(cid:12) 4π Ψ δ3(r ) Ψ = 4 Ψ (E V) Ψ 2 ~ Ψ ~ Ψ . (11) a Ψ c c (cid:28) (cid:12)r − (cid:12) (cid:29)− (cid:28)∇ (cid:12)r (cid:12)∇ (cid:29) (cid:10) (cid:12) (cid:12) (cid:11) (cid:12) a (cid:12) Xc (cid:12) a(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) whereV isatotalpotentialenergy inEq. (5). 5 III. LEADINGRELATIVISTICCORRECTIONTOBINDINGENERGY The leading relativistic corrections to energy levels are given by the expectation values of the Breit-Pauli HamiltonianH(4). p~4 πZα Zα ~r H(4) = a + δ3(r )+ ~σ a p~ a a a (cid:26)−8m3 2m2 4m2 · r3 × (cid:27) Xa a πα α δij ri rj + δ3(r ) pi + ab ab pj (cid:26)−m2 ab − 2m2 a(cid:18)r r3 (cid:19) b Xa>b Xb ab ab 2πα α σi σj ri rj α ~σ ~σ δ3(r )+ a b δij 3 ab ab + −3m2 a · b ab 4m2 r3 (cid:18) − r2 (cid:19) 4m2r3 ab ab ab 2 ~σ ~r p~ ~σ ~r ~p + ~σ ~r p~ ~σ ~r ~p . (12) a ab b b ab a b ab b a ab a × (cid:20) · × − · × · × − · × (cid:21)(cid:27) (cid:0) (cid:1) (cid:0) (cid:1) ForstateswithvanishingangularmomentumLandspinS = 1/2,theexpectationvalueissimpli- fied to theform p~4 πZα E(4) = Ψ H(4) Ψ = a + δ3(r ) a h | | i (cid:28) (cid:26)−8m3 2m2 (cid:27) Xa πα α δij ri rj + δ3(r ) pi + ab ab pj . (13) (cid:26)m2 ab − 2m2 a(cid:18)r r3 (cid:19) b(cid:27)(cid:29) Xa>b Xb ab ab E(4) has already been obtained in works [15, 17]. Calculations of these matrix elements involves the usual Hylleraas integrals with all n nonnegativeand extended integrals, namely with one pa- i rameter n equal to 1. The direct numerical method to calculate these integrals was presented i − in [8, 9]. Here we apply the analytic approach. Recursion relations for the case of n or n or n 1 2 3 equal to 1 have been obtained in [13]. Hylleraas integrals involving n or n or n equal to 1 4 5 6 − − can in principle be obtained by the integration of the usual Hylleraas integral with respect to the corresponding parameter w [13]. However, some recursion relations may become unstable, for i exampleinthecase ofn = 1 therecursionin n is numericallyunstableforlargew . To avoid 4 1 1 − this problem we derive in the next section stable recursion relations for extended Hylleraas inte- gralswithn = 1 fori = 4,5,6. NumericalresultsformatrixelementsoftheBreitHamiltonian i − using these recursion relations, has been presented in Table I and II. One observes that the lowest convergenceisforthe p4/8term,andinspiteofthedifferencesforseparatematrixelements,the − totalrelativisticcorrection isin goodagreement withtheformerresult in[15]. 6 TABLEII:MatrixelementsoftheBreit-PauliHamiltonianH(4) inatomicunits. Ω 1 4 1 i δij + raibrajb j H(4) a−8∇a a>b 2 ∇a rab ra3b ∇b P P (cid:0) (cid:1) 3 -78.58728669090 -0.43863254584 -12.08067033680 4 -78.55733185961 -0.43609658640 -12.05311194461 5 -78.55635590597 -0.43569734491 -12.05070911655 6 -78.55671450343 -0.43561642650 -12.05038807638 7 -78.55619578085 -0.43560236202 -12.05000429451 8 -78.55616264213 -0.43559952390 -12.04996116216 9 -78.55613747761 -0.43559821744 -12.04992614976 10 -78.55613573401 -0.43559804758 -12.04992141427 11 -78.55613159634 -0.43559796357 -12.04991680081 12 -78.55612863210 -0.43559791050 -12.04991377296 -78.55611288(200) -0.435597765(50) -12.04989786(200) ∞ Ref. [15] -78.55613555(148) -0.435598001(137) -12.04990994(180) IV. RECURSION RELATIONS FOR THREE-ELECTRON EXTENDED HYLLERAAS INTE- GRALWITH1/r2 1 Intheformersectionwecalculatedrelativisticcorrections. Forthisweneededvariousextended Hylleraas integrals, among them, integrals with 1/r2, which are being derived here. To obtain i recursion relations for three-electron Hylleraas integral in Eq. (7), one first considers the integral G 1 G(m ,m ,m ;m ,m ,m ) = d3k d3k d3k (k2 +u2)−m1 (k2 +u2)−m2 1 2 3 4 5 6 8π6 Z 1Z 2Z 3 1 1 2 2 (k2 +u2)−m3(k2 +w2)−m4(k2 +w2)−m5(k2 +w2)−m6, (14) 3 3 32 1 13 2 21 3 which is related to f by: f(0,0,0,0,0,0) = G(1,1,1,1,1,1) . The following 9 inte- |u1=u2=u3=0 gration by part identitiesare valid because the integral of thederivativeof a function vanishingat infinityvanishes, ∂ 0 id(i,j) = d3k d3k d3k ~k (k2 +u2)−1 ≡ Z 1Z 2Z 3 ∂~kih j 1 1 (k2 +u2)−1(k2 +u2)−1(k2 +w2)−1(k2 +w2)−1(k2 +w2)−1 , (15) 2 2 3 3 32 1 13 2 21 3 i 7 wherei,j = 1,2,3. Thereductionofthescalarproductsfromthenumeratorleadstotheidentities forthelinearcombinationoftheGfunctions. Ifanyoftheargumentsisequalto0,thenGbecomes aknowntwo-electronHylleraastypeintegral. Theseidentitiesareusedtoderivevariousrecursion relations. Here, we derive a set of recursions for the case when n , n or n is equal to 1. Let 4 5 6 − us assume that n = 1. The analytic expression for f(0,0,0, 1,n ,n ) involves powers of 4 5 6 − − w w inthedenominatorwhichisnotveryconvenientinhighprecisionnumericalcalculations. 2 3 − Instead, we use recursions for f(0,0,0,0,n ,n ) and numerically integrate with respect to w , 5 6 1 namely ∞ f(0,0,0, 1,n ,n ) = dw f(0,0,0,0,n ,n ). (16) 5 6 1 5 6 − Z w1 These recursions are derived as follows. We take id(i,i) with i = 1,2,3 and put u = 0. Result- i ing three equations are solved against three unknowns: G(1,1,1,2,1,1), G(1,1,1,1,2,1), and G(1,1,1,1,1,2). Thesolutionforthelast twoG functionsisthefollowing 1 G(1,1,1,1,2,1) = G(0,1,1,1,1,2) G(1,0,1,1,1,2) G(1,0,1,2,1,1) w2 − − 2 (cid:2) +G(1,1,0,2,1,1)+G(1,1,1,1,1,1)/2 , (17) 1 (cid:3) G(1,1,1,1,1,2) = G(0,1,1,1,2,1)+G(1,0,1,2,1,1) G(1,1,0,1,2,1) w2 − 3 (cid:2) G(1,1,0,2,1,1)+G(1,1,1,1,1,1)/2 . (18) − (cid:3) By differentiationwithrespect tow and w oneobtainsthefollowingrecursion relations 2 3 1 f(0,0,0,0,n +1,n ) = (n +1)f(0,0,0,0,n ,n )w w 5 6 5 5 6 1 3 w w w 1 2 3 (cid:2) (n +1)n f(0,0,0,0,n ,n 1)w 5 6 5 6 1 − − +n f(0,0,0,0,n +1,n 1)w w 6 5 6 1 2 − n Γ(n ,n 1, 1,w +w ,w ,0) 6 5 6 1 2 3 − − − +n Γ(n 1,n , 1,w +w ,w ,0) 6 6 5 1 3 2 − − Γ(n ,n , 1,w +w ,w ,0)w 6 5 1 3 2 1 − − +Γ(n +n ,0, 1,w +w ,w ,0)w 5 6 2 3 1 1 − +Γ(n ,n , 1,w +w ,w ,0)w 5 6 1 2 3 3 − Γ(n ,n , 1,w +w ,w ,0)w , (19) 6 5 1 3 2 3 − − (cid:3) 8 1 f(0,0,0,0,n ,n +1) = (n +1)f(0,0,0,0,n ,n )w w 5 6 6 5 6 1 2 w w w 1 2 3(cid:2) n (n +1)f(0,0,0,0,n 1,n )w 5 6 5 6 1 − − +n f(0,0,0,0,n 1,n +1)w w 5 5 6 1 3 − +n Γ(n 1,n , 1,w +w ,w ,0) 5 5 6 1 2 3 − − n Γ(n ,n 1, 1,w +w ,w ,0) 5 6 5 1 3 2 − − − Γ(n ,n , 1,w +w ,w ,0)w 5 6 1 2 3 1 − − +Γ(n +n ,0, 1,w +w ,w ,0)w 5 6 2 3 1 1 − Γ(n ,n , 1,w +w ,w ,0)w 5 6 1 2 3 2 − − +Γ(n ,n , 1,w +w ,w ,0)w . (20) 6 5 1 3 2 2 − (cid:3) whereΓ is aknown[18, 19, 20]two-electronintegral d3r d3r Γ(n ,n ,n ,α ,α ,α ) = 1 2 e−α1r1−α2r2−α3r12rn1−1rn2−1rn3−1. (21) 1 2 3 1 2 3 Z 4π Z 4π 1 2 12 The integration in Eq. (16) is performed numerically using adapted points and weights to the functionwhichhas logarithmicend-pointsingularity,namely 1 dx W (x)+W (x) ln(x) , (22) 1 2 Z 0 (cid:2) (cid:3) where W are functions without any singularities. The method to obtain n adapted points and i weights is presented in Appendix A, and this integral is exact for W being polynomialsup to the i ordern 1. Intheactualcalculationsweachievedatleast48digitsprecisionwithonly100points. − Having obtained f(0,0,0, 1,n ,n ) we construct recursion relations in n , n , and n . This is 5 6 1 2 3 − achievedintwosteps. InthefirststepweuseintegrationbypartsinmomentumrepresentationEq. (15), toform thefollowinglinearcombination id(2,2)+id(3,3) id(1,1) = 2 G(0,1,1,1,1,2)+G(0,1,1,1,2,1) G(1,0,1,1,1,2) − − (cid:2) G(1,1,0,1,2,1) G(1,1,1,1,1,1)/2 G(2,1,1,1,1,1)u2 − − − 1 G(1,1,1,1,1,2)(u2 u2)+G(1,2,1,1,1,1)u2 − 1− 2 2 G(1,1,1,1,2,1)(u2 u2)+G(1,1,2,1,1,1)u2 − 1− 3 3 +G(1,1,1,2,1,1)w2 = 0. (23) 1 (cid:3) 9 We integrate with respect to w and differentiate over u , u , u , w , and w to obtain the main 1 1 2 3 2 3 formula 1 f(n ,n ,n , 1,n ,n ) = 1 2 3 5 6 − (n +n n )w w 2 3 − 1 2 3 (cid:2) (n 1)n n f(n 2,n ,n , 1,n 1,n +1) 1 1 5 1 2 3 5 6 − − − − +(n 1)n n f(n 2,n ,n , 1,n +1,n 1) 1 1 6 1 2 3 5 6 − − − − (n 1)n n f(n ,n 2,n , 1,n 1,n +1) 2 2 5 1 2 3 5 6 − − − − − (n 1)n n f(n ,n ,n 2, 1,n +1,n 1) 3 3 6 1 2 3 5 6 − − − − − +(n n n )n n f(n ,n ,n , 1,n 1,n 1) 1 2 3 5 6 1 2 3 5 6 − − − − − +n n f(n ,n ,n ,0,n 1,n 1)w 5 6 1 2 3 5 6 1 − − (n 1)n f(n 2,n ,n , 1,n ,n +1)w 1 1 1 2 3 5 6 2 − − − − +(n 1)n f(n ,n 2,n , 1,n ,n +1)w 2 2 1 2 3 5 6 2 − − − (n n n )n f(n ,n ,n , 1,n ,n 1)w 1 2 3 6 1 2 3 5 6 2 − − − − − n f(n ,n ,n ,0,n ,n 1)w w 6 1 2 3 5 6 1 2 − − (n 1)n f(n 2,n ,n , 1,n +1,n )w 1 1 1 2 3 5 6 3 − − − − +(n 1)n f(n ,n ,n 2, 1,n +1,n )w 3 3 1 2 3 5 6 3 − − − (n n n )n f(n ,n ,n , 1,n 1,n )w 1 2 3 5 1 2 3 5 6 3 − − − − − n f(n ,n ,n ,0,n 1,n )w w 5 1 2 3 5 6 1 3 − − +f(n ,n ,n ,0,n ,n )w w w 1 2 3 5 6 1 2 3 +n δ(n )Γ(n 1,n 1,n +n 1,w +w ,w ,0) 6 3 5 6 1 2 1 2 3 − − − +n δ(n )Γ(n 1,n 1,n +n 1,w +w ,w ,0) 5 2 6 5 1 3 1 3 2 − − − n δ(n )Γ(n +n 1, 1,n +n 1,w +w ,w ,0) 5 1 5 6 2 3 2 3 1 − − − − n δ(n )Γ(n +n 1, 1,n +n 1,w +w ,w ,0) 6 1 5 6 2 3 2 3 1 − − − − δ(n )Γ(n ,n 1,n +n 1,w +w ,w ,0)w 2 6 5 1 3 1 3 2 2 − − − +δ(n )Γ(n +n , 1,n +n 1,w +w ,w ,0)w 1 5 6 2 3 2 3 1 2 − − δ(n )Γ(n ,n 1,n +n 1,w +w ,w ,0)w 3 5 6 1 2 1 2 3 3 − − − +δ(n )Γ(n +n , 1,n +n 1,w +w ,w ,0)w . (24) 1 5 6 2 3 2 3 1 3 − − (cid:3) Thisgeneralformuladoesnotworkinthecasen = n +n . Inthesecondstepweuseintegration 1 2 3 by partidentitiesinthecoordinatespace tofill thishole. We limitourselvesonlytoaspecial case 10