1 Introduction

It is now well accepted that accurate predictions of electronic properties of heavy atoms cannot be achieved without the introduction of relativistic corrections. To include these corrections the relativistic counterpart of the N-electron Schroedinger equation would be an extension of the two particle Salpeter and Bethe [1] equation generalized to include an external field [2]. Unfortunately such an equation suffers from two problems: firstly, it involves an integral kernel that cannot be written in a closed form; and, secondly it depends on more than one time variable and consequently cannot be reduced to a Hamiltonian form. The first of these problems originates from the fact that any relativistic quantum theory conserves only the total charge of the system but not the number of particles. The second problem arises from the requirement of a fully relativistic covariance. Because of these difficulties, almost all practical calculations are done using an effective electron-electron interaction.

The most commonly used approximation is the so called ”no pair” approximation in which the effective Hamiltonian explicitly excludes electron-positron pairs. Thus writing formally the Hamiltonian of a N electron system as

(1)

only the first term of the expansion is taken into account. Then, as in the non relativistic case, we consider a sum of one-electron and two-electron interactions. Both of these interactions include an expansion in term of the coupling constant (i.e. the fine structure constant ). The choice of the Dirac one-electron operator enables to sum the whole series in powers of for the one-electron interaction with the nuclear field considered as a classical external field (we neglect at this stage field quantization that gives rise to quantum electrodynamics corrections). The electron-electron interaction must obviously be deduced from quantum electrodynamics (QED). This is generally achieved in the bound state picture of the S-matrix with the result that, after time integration, the interaction reduces to a generalized form of the Breit operator [3].

Following the approach outlined above, a computer code has been written to handle relativistic calculations for atoms. The present version is an update of the code first published in 1975 [4]. This version make use of the multiconfiguration approximation for the N electron wavefunction. In this scheme both one-electron orbitals and mixing coefficients between the various configurations are fully optimised during a double step self consistent process. In the first step, the mixing coefficients between the configurations are kept fixed while the radial equations are solved self-consistently. The second step optimises the mixing coefficients between the configurations through a diagonalization of the Hamiltonian matrix.

This chapter is organised as follows: the first section presents a short summary of the multiconfiguration Dirac-Fock method (MCDF). Section 2 describes the reduction of the expectation value of the Hamiltonian to radial integrals. Section 3 is devoted to the numerical solution of the self consistent integro-differential equations. The radiative corrections are surveyed in section 4 and the last section provides a short description of the options available in the package.

2 Summary of the MCDF method

2.1 Wavefunctions

As in many self-consistent field method, the starting point to build up N-electron wavefunctions are central field one-electron orbitals, thus in the relativistic case use is made of the Dirac four component spinors:

(2)

these spinors are simultaneous eigenfunctions of the parity operator, the square of the total angular momentum operator j = l + s and of its z-component. The quantum numbers n,,m are respectively the principal quantum, the relativistic ”kappa” quantum number that includes both the parity and the total angular momentum j and m is the eigenvalue of j_z. The quantum number j,,l are related through:

m

(4)

where the terms under the summation are respectively a Clebsch-Gordan coefficient, a spherical harmonic and one of the two Pauli spinors:

(5)

The Dirac spinors as defined by Eq.(2) are used to build Configuration State Functions (CSF) | JM > that are antisymmetric products of the one-electron wavefunctions and are eigenvalues of the parity , the total angular momentum J and its projection M. The label stands for all other values (angular momentum recoupling scheme, seniority numbers, ...) necessary to define unambiguously the CSF. For a N-electron system, a CSF is thus a linear combination of Slater determinants:

(6)

all of them with the same and M values while the d_i’s are determined by the requirement that the CSF is an eigenstate of J². The method used in the package to build these eigenstates is described in Appendix A. The total MCDF wavefunction is constructed as a superposition of CSF’s, i.e.

(7)

2.2 Hamiltonian

As pointed out in the introduction, the relativistic effective Hamiltonian for an N-electron atomic system is taken as the Dirac Breit Hamiltonian given by:

(8)

with for the Dirac operator

(9)

here and throughout this chapter we make use of atomic units. Thus e = h/2 = m_e = 1 and c = 1/ where e is the absolute value of the electron charge, h the Plank constant, m_e the mass of the electron, c the speed of light and 1/137 the fine structure constant. The 4x4 Dirac matrices _i and are given by

(10)

in terms of the 2x2 Pauli matrices _i and the 2x2 zero and unit matrix I. V _N(r) is the nuclear potential. For the electron-electron interaction the Breit operator

(11)

includes the instantaneous Coulomb repulsion, the magnetic interaction and the retardation in the electron-electron interaction due to the finite value of the speed of light. ₁₂ is the frequency of the exchange photon and, in an independent model description is given by , the ’s being the electron binding energies. At this point it is worth to point out that:
- the electron-electron interaction is expressed in the Coulomb gauge and this the correct one to use as long as a finite summation (in terms of the number of photon lines) is included in the framework of the Dirac-Fock potential. As pointed out first by Love [5] and more recently studied by Sucher [6] and Lindgren [7] the use of another gauge, like the Lorentz one, will introduce at each order of the perturbation series a spurious contribution that will be cancelled at the next order of the expansion. Obviously this gauge dependence arises only from the approximation inherent to any practical calculation and furthermore is specific to the approximation used. For example when a local potential is used to define the single particle states, the result is gauge independent [8]. This gauge problem is indeed the same as the one, well known in non relativistic calculations, that occurs when computing dipole transitions in the so-called length and velocity forms.
- the retardation part beyond the zero frequency limit (i.e. the limit obtained by setting ₁₂ = 0 in Eq.(11)) can only be used at the Dirac-Fock level since Koopman’s theorem can be used to obtain an approximation of the binding energies. On the other hand as soon as more than one configuration is included, the MCDF method does not provided any easy way to get approximated one-particle energies. This is an example of one yet unsolved problem in the framework of the multiconfiguration approach to relativistic correlation energies.
- the Hamiltonian as given in Eq.(8) is the so-called ”no-pair” approximation in the sense that it excludes explicitly any contribution from electron-positron pairs. Consequently care as to be taken to avoid contribution from the negative continuum. As no projection operators are included in the present version of the package, problems of convergence may show up for large values of the nuclear charge and/or correlated orbitals of high principal quantum numbers [9]. An updated version of the program (under development) will try to take care of this problem.

3 Reduction to radial integrals

From what has been discussed in the previous section, the expectation value of the total Hamiltonian will include contributions from:
- the one particle Dirac Hamiltonian
- the Coulomb repulsion
- the Breit interaction
that, using A,B,C,... as label for the one electron orbitals, can be written respectively as

(12)

where the three radial integrals are respectively the contributions of substracting the rest mass energy, the kinematic operator ^. and the nucleus-electron interaction.

For the instantaneous Coulomb repulsion, the usual multipole expansion of 1/r₁₂ leads to:

(13)

where U_k(1,2) = , r_< being the smaller of (r₁,r₂) and r_> the larger. The d^k’s are angular coefficients given, in terms of Wigner 3j symbols, by:

(14)

for which the following triangular conditions must hold:

(15)

while l_A + l_C + k and l_B + l_D + k must both be even.

As given in Eq.(11) the Breit interaction includes the two contributions:

(16)

The reduction to radial integrals gives for the first term V _Br¹():

(17)

with the triangular conditions:

(19)

thus is the orbital quantum number for the small component. At this stage it is also worth to remember that _small = -_large. Furthermore

(20)

The radial integral itself is given by

(21)

j_k and _k are the regular and irregular spherical Bessel functions.

The second term V _Br²() of the Breit operator reduces to:

(22)

with the same triangular and parity conditions as for the B¹ integrals. The radial integrals are here defined as:

(23)

with

(24)

3.1 Zero frequency limit of the Breit operator

When = 0 the two terms of Eq.(16) reduce to

(25)

Using the expansions for the Bessel fonctions

(26)

the radial integral of Eq.(21) reduces to:

(27)

and the definition of the radial operator W simplifies to:

(28)

3.2 Average energy

As in the non-relativistic case, it is useful to introduce the concept of the average energy defined as the centre of gravity of all the levels belonging to a given (jj) configuration and given by

(29)

The usefulness of the average energy is due to the fact that the energy of a given level is just the sum of E_Av plus a small number of extra two electron integrals. As the expression of E_Av is very easy to incorporate in any computer code, this make the input data for a given level rather convenient to specify. The average energy is most easily obtain from the interaction of a single electron with a closed shell. Using the well known closure properties of the Wigner 3j symbols:

(30)

we obtain for the coefficients of the diract Coulomb integrals

(31)

and for the exchange Coulomb integrals:

(32)

For the first contribution to the Breit operator, the triangular relations of Eq.(18) combined with the result of Eq.(31) show that there is no contribution from direct integrals as expected since a closed shell has a zero net magnetic moment. For the exchange integrals we obtain:

(33)

where, using the definition of Eq.(27) for the radial integrals, we have defined:

(34)

and

(35)

With the help of the above relations we can now write the explicit expression of the average energy that, restricted to Coulomb interaction for simplicit, reads:

(36)

the q’s are the occupation numbers of the orbitals and the one-electron integrals I are those defined by Eq.(12). The direct Coulomb integrals F^k are the radial integrals of Eq.(13) in which C = A and D = B while the exchange integrals G^k are obtained from the same equation by setting C = B and D = A.

4 MCDF integro-differential equations

4.1 Radial equations

For simplicity we use only orthonormalized Dirac spinors according to:

(37)

Then the MCDF equations are obtained by writing E = / and using a variational principle on both the c coefficients, and on the large and small radial components P_n(r) and Q_n(r). Variation of the c leads to a Hamiltonian matrix from which the c are deduced by diagonalization. Variation with respect to the radial wavefunctions leads to an integro-differential equation that for a given orbital A reads:

(38)

where V _A is the sum of the nuclear potential and the direct Coulomb potential while X_PA and X_QA include all the two-electron interactions excepted for the direct Coulomb instantaneous repulsion. _A,B are Lagrange parameters used to enforce the orthogonality constraints of Eq.(37) and thus the summation over B run only for orbitals with _B = _A.

4.2 Numerical solution of the radial equations.

Most of the numerical methods used to solve the relativistic Hartree-Fock equations are drawn from what has been known for a long time in the non relativistic case (see [10]). The integro-differential equations are solved by an iterative process that reduces, at each step, the integro-differential equations to simple differential equations by considering the potential terms (direct and exchange) as given source functions computed from the wave function obtained from the previous step. The iterative process is then continued until a given accuracy is reached between two successive iterations. Obviously the need to solve a pair of two first order coupled differential equations instead of a single second order differential equation will introduce some modifications in the recipes available since the pioneer work of Hartree [11].

4.2.1 Predictor Corrector Methods.

These methods are based on the following strategy: to obtain the function at the mesh point n+1 use is made of the known values at the preceding N points to predict a first estimate, this estimate is inserted in the differential equation to obtain the derivative that in turn is used to correct the first estimate, then the final value may be taken as a linear combination of the predicted and corrected values to increase the accuracy. As an example we consider the five points Adams’ method that has been widely selected because of its stability properties [12]. The predicted, corrected and final values are given respectively by:

4.2.2 Finite Differences Methods

As seen in the previous section, the predictor corrector methods require to solve two systems of equations (the homogeneous and inhomogeneous ones) to obtain a continuous solution while in the non relativistic case the Numerov method associated with tail correction [10] provides directly a continuous approximation (the derivative remains discontinuous until the eigenvalue is found). We consider now alternative methods that easily allow to enforce the continuity of one of the two radial components. Let us define the solution at point n+1 as:

(41)

where _n is a difference correction given, in terms of central differences, by:

(42)

with:

(44)

where r^' stands for dr/dt (to take into account the fact that the tabulation variable t may be a function of r) and X^P(Q) = X_PA(QA) + _BAA,BP_B(Q_B). Furthermore, all the functions of r are calculated, as already said, with the help of the of the wavefunctions obtained at the previous iteration. Then the system of algebraic equations:

(46)

and the two column vectors (PQ) and (uv) defined as:

(47)

as displays in Eq.(46) each row of the matrix M has at most four non zero elements. To solve this system of equation the matrix M is decomposed into the product of two triangular matrices M = L^.T in which L is a lower matrix with only three non zero elements on each row and T an upper matrix with the same property. Introducing an intermediate vector (pq) it is possible to solve L(pq) = (uv) for m, m+1, .., N and then T(PQ) = (pq) for N, N-1, ..,m. The last point of tabulation N is determined by the requirement that P_N should be lower than a specified small value when assuming Q_N = 0. Thus the number of tabulation points of each orbital is determined automatically during the self consistency process. This elimination process produces, as written here, a large component P that is continuous everywhere. The discontinuity of the small component at the matching point r_m can then be used to adjust the eigenvalue _A,A.

4.2.3 Implementation in the package

Because the bound orbitals exhibit a rapid variation near the origin and an exponential decay at large values of r it is more efficient to make the change of variable:

(48)

to tabulate all the radial functions on this grid of points equally spaced with respect to this new variable t. The simple logarithmic distribution of the tabulation points is well adapted for not too high principal quantum numbers n (typically for all occupied orbitals throughout the periodic table) but will not give enough points for high n Rydberg orbitals. For this reason we have choosen the second form that is also convenient for continuum orbitals.

As we have to solve the Dirac-Fock radial equations for different purposes (to get initial starting wavefunctions, during the self-consistent multiconfiguration process or to obtain hydrogenic solutions for QED corrections, see section 5) we have implemented two solvers. The first one is designed for homogeneous equations and used the predictor corrector method described in section 4.2.1. Indeed, for the homogeneous equation, it is always possible to get one the two radial functions continuous while the discontinuity in the other one is used to adjust the eigenvalue. For example if the large component P is made continuous, then the first order correction to the eigenvalue is given by:

(49)

where R is the matching point for the inward and outward integrations. This method has the advantage to be both fast and rather accurate.

As described in section 4.2.2 the finite difference method is well adapted for Dirac-Fock orbitals but less for correlation ones. Indeed it is rather difficult to obtain a good estimate of these correlation orbitals and the rather poor intrinsic accuracy of the finite difference method may result in instabilities during the self consistent field process. We thus proceed in the following way: the outward integration is performed using the predictor corrector method while the outward integration uses the finite difference method to insurance continuity of the large component. But to increase the accuracy an extra loop is introduced in the inward integration scheme in order to compute the central differences of Eq.(43) from the wave function being solved.

5 Radiative corrections

Calculation of many-electron radiative corrections is still one the most difficult numerical problem to deal with for high-precision level predictions. Even if recent progresses have been achieved [13, 14, 15], it is rather difficult to include them as standard procedures in distributed computer codes. Nevertheless as, at least an estimation of the order of magnitude of these corrections, has to be included in present days relativistic calculations, we have provided some estimations for them. These corrections include two contributions:
1) the vacuum polarization, for which effective potentials can be deduced from Quantum Electro-Dynamics (QED)
2) the self-energy, the first-order contribution of which is known exactly (i.e. to all orders in the external field strength Z) for hydrogen-like systems [16, 17, 18]

5.1 Vacuum polarization

This contribution is related to the creation and annihilation of virtual electron-positron pairs in the field of the nucleus. The first term of order (Z) with respect to the nuclear Coulomb potential can be evaluated, following Wichmann and Kroll [19], as the expectation value of the Uehling potential [20]. In the case of a nuclear spherical charge distribution _N(r) this potential is given by:

(50)

where _e is the Compton wavelength of the electron and the function K₀(x) is defined as:

(51)

Terms of higher order can also be evaluated as the expectation value of potentials. The correction of order ²(Z) is associated to the Kallen and Sabry [21] potential while the term of order (Z)³ has been given by Blomqvist [22]. Various analytical approximations [23, 24] are available to calculate the various functions like K₀ in Eq.(51) involved in the vacuum polarization potentials of various orders.

The contribution of the other bound electrons to the vacuum polarization can be estimated in analogy with that of the nucleus by replacing _N(r) in Eq.(50) by the electronic charge of the electrons.

5.2 Self-energy

This second contribution arises from the interaction of the electron with its own radiation field. For one-electron systems this self self-energy term has been calculated exactly by P.J. Mohr [25, 17, 18] and expressed as:

(52)

where F(Z) is a slowly varying function of Z. Previous approximations of the self-energy in many-electron atoms used an effective Z value to scale the above hydrogenic results, the effective Z being deduced from the comparison of the mean value of the radius of the Dirac-Fock orbital with that of the hydrogenic one. It has been shown [26] that this effective Z approximation may be completly wrong for excited states and must be replaced by a more physical prescription to obtain the effective Z. This prescription is based upon the Welton picture of the Lamb-shift [27] in which the self-energy is due to the perturbations induced in the electron classical trajectory by the fluctutation of the electromagnetic field of the vacuum. Driven by these flucuations, the electron explores the field of the nucleus over a region of extension r, its trajectory being defined by r. The electron then probes a potential V _N(r + r) instead of the unperturbed potential V _N(r). The perturbation potential then reads:

(53)

as the mean value of fields in the vacuum are zero, < r >= 0 and the perturbation reduces to:

(54)

In this picture < (r)² > is infinite and must be renormalized. Taking advantage that exact hydrogenic results are available, the screened self energy correction for s orbitals can be obtained from the following relation:

(55)

where the subscripts DF and Hyd stand for the Dirac-Fock and Hydrogenic Dirac orbitals. If a finite nuclear size is used, this method automatically includes the volume effect in the self energy. For non s orbitals, the main contribution to the self energy is due to the anomalous gyromagnetic ratio of the electron because the electronic density near the nucleus is much smaller than for s orbitals.

We illustrate in the two following tables the respective accuracy of the various approximations to the screening of the self-energy. In the low Z region, the Kabir and Salpeter equation [28] provides a good approximation since relativistic corrections to the radial wavefunctions are small. This equation has been used by De Serio et al. [29] for various nl,n^'l^' levels of helium-like ions.

^(a) Ref. [29]

The comparison, given in the above table, of the self-energy contribution to the energy splittings between the ³P_J ans ³S₁ levels illustrates the poor performance of the effective Z method to correct for the screening. At high Z nuclear charges, accurate results have been recently obtained by Blundell [15] for lithium-like uranium and the next table show the comparison between his accurate values and the previous semi-empirical or phenomenological estimates of the screening corrections.

^(a) Ref. [15]

6 Options of the package

The purpose of this section is not to provide a write up for the package but only to outline some specific options of the code. Besides the fact that, after self consistency has been achieved, radiative corrections are added as a first order correction to the total energy, various options are made available for both the nuclear model and the optimisation of the one-electron wavefunctions.

6.1 Nuclear model

For the inner shell of heavy atoms and more generally for all their s orbitals, it is mandatory to go beyond the point nucleus approximation. For that purpose it is possible to use two nuclear models. The first one assumes that the nuclear charge distribution can be approximated by a uniform charge sphere of radius R_nuc, thus for an atom of atomic number Z the nuclear potential is given by:

(56)

A more realistic description of the proton charge distribution inside the nucleus is provided by the Fermi distribution:

(57)

where t is the thickness parameter of the Fermi distribution and ₀ is a normalisation constant such that the integral of the proton charge equals Z. As the detailed proton charge distribution becomes a critical parameter only for very heavy atoms, the default option of the program is to use an uniformed charged sphere for atomic numbers lower then 45 and a Fermi distribution with a thickness parameter t = 2.3 otherwise. Furthermore the code include a table of nuclear radii as given by W.R. Johnson and G. Soff [30]. This standard option can be overwritten by providing specific nuclear parameters for the atom under consideration. As an illustration of the influence of the finite nuclear size, we give in Table III the values of this correction for the low lying hydrogenic energy levels.

6.2 Self consistent process

Two options are available to achieve self consistency. For the first one, both the mixing coefficients between the (jj) configurations and the radial orbitals are fully optimised and such a calculation is what is called MCDF. In the second option only the radial orbitals are optimised while the weights of the configurations are kept fixed. When the weights are frozen, an useful concept to introduce (and included in the code) is the definition of an extended average energy (EAL). This concept introduce almost simultaneously by various authors [31, 32] has the property than the EAL exactly reduces to the non relativistic average energy. The EAL is defined as the sum of the (jj) average energy defined in Eq. (29), sum that includes all the (jj) subconfigurations arising from a single (LS) configuration. In the sum, each of the (jj) subconfiguration is weighted by its degeneracy. As an example, consider the p² configuration that, in the relativistic case, splits into the three subconfigurations : p_1/2², p_1/2p_3/2, p_3/2². The EAL is then defined as:

(58)

Whatever is the choice selected (MCDF or EAL) for the self consistent process, it is possible to include either only the instantaneous Coulomb repulsion or both this repulsion and the zero frequency limit of the first term of the Breit interaction as given by Eq. (25) also known as the Gaunt interaction. In fact there is no fundamental reason to include only the first term of the Breit interaction since both term have the same Z dependence. On the other hand the numerical value of this first term is typically one order of magnitude larger than the second one given in Eq. (25). In so doing one save computer time and still get the main contribution of the Breit interaction to the wavefunctions. Table IV illustrates the changes in the radial wavefunctions of uranium Beryllium-like when the Gaunt interaction is introduced in the self consistent field process.

(a) Coulomb only (b) Coulomb + Gaunt

To conclude, let us just point out that, contrary to the non relativistic approximation, the multiconfiguration approach is required, in the relativistic case, not only to include correlation effects but also to handle intermediate coupling scheme as soon as open shells are involved.

Acknowledgements

The author is indebted to Dr P. Indelicato for his many contributions to the improvements in the code and especially in the calculation of radiative corrections. The Welton approximation of the screening correction to the self energy belongs to him.

The method we used [33] is unsophisticated in the sense that it did not use any fractional parentage coefficients but it is easy to implement in a computer code and is quite general. Indeed the only limitation is the size of the matrix to be diagonalized. We may also note that the same approach has been used in the non relativistic case [34]. The wavefunction being written as a sum of Slater determinants all of them having the same and M values it is easy to calculate the matrix of J² and then to diagonalize it to obtain the eigenstates of the total angular momentum (i.e. the coeffiecients d_i of Eq.(6)). To compute the matrix of J² first remember that:

(59)

Then, as the one-electron orbitals are assumed to be orthonormalized, it is enough to compute the expectation value of J² between two determinants in the three following cases:
1) the two determinants are the same

(60)

with the short hand notation:

(61)

2) the two determinants are different by one orbital (A A^')

(62)

3) the two determinants are different by two orbitals (A A^') and (B B^')

(63)

References