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.

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:

| (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^{2}. 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) |

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. _{12} 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 _{12} = 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.

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_{12} leads to:

| (13) |

where U_{k}(1,2) = , r_{<} being the smaller of (r_{1},r_{2}) 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}^{1}():

| (17) |

with the triangular conditions:

while l
| (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}^{2}() of the Breit operator reduces to:

| (22) |

with the same triangular and parity conditions as for the B^{1} integrals. The radial integrals are here defined
as:

| (23) |

with

| (24) |

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) |

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.

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}.

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].

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:

where p

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:

at this stage we may remember that the accuracy of the solution is required only when self consistency is reached. Consequently, the difference correction
| (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)} +
_{BA}_{A,B}P_{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}.

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.

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]

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_{0}(x) is defined as:

| (51) |

Terms of higher order can also be evaluated as the expectation value of potentials. The correction of order
^{2}(Z) is associated to the Kallen and Sabry [21] potential while the term of order (Z)^{3} has been
given by Blomqvist [22]. Various analytical approximations [23, 24] are available to calculate the
various functions like K_{0} 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.

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)^{2} > 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.

Table I: Self-energy screening contribution (in eV) to the ^{3}P_{J} ^{3}S_{1} transitions of helium-like ions.

^{3}P_{0} -^{3}S_{1} | ^{3}P_{2} -^{3}S_{1} | |||||

Z | Z_{eff} | Welton | DeSerio^{(a)} | Z_{eff} | Welton | DeSerio^{(a)} |

12 | 0.0007 | 0.0046 | 0.0048 | 0.0002 | 0.0040 | 0.0034 |

18 | 0.0016 | 0.0105 | 0.0116 | 0.0001 | 0.0089 | 0.0077 |

26 | 0.0038 | 0.0258 | 0.0317 | 0.0001 | 0.0216 | 0.0203 |

^{(a)} Ref. [29]

The comparison, given in the above table, of the self-energy contribution to the energy splittings
between the ^{3}P_{J} ans ^{3}S_{1} 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.

Table II: Screening correction (in eV) to the self-energy in lithium-like Uranium

Orbital | Z_{eff} | Welton | Blundell^{(a)} |

2s_{1/2} | 3.71 | 2.88 | 2.75(2) |

2p_{1/2} | 1.21 | 0.84 | 0.97(2) |

^{(a)} Ref. [15]

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.

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 _{0} 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.

Table III: Finite nuclear size correction (in eV) for hydrogenic levels.

Numbers in parentheses are power of 10

Z | 1s | 2s | ^{2}p_{1/2} | ^{2}p_{3/2} |

25 | 0.0436 | 0.0056 | 4.(-5) | 0. |

50 | 1.9580 | 0.2556 | 0.0075 | 0. |

75 | 42.399 | 5.3110 | 0.3634 | 0. |

100 | 451.23 | 92.965 | 13.922 | 6.(-5) |

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^{2} configuration that, in the relativistic
case, splits into the three subconfigurations : p_{1/2}^{2}, p_{1/2}p_{3/2}, p_{3/2}^{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.

Table IV: Change in the radial expectation values < r^{n} > due to the Gaunt interaction (uranium
Beryllium-like). The calculation includes the two 1s^{2}2s^{2} and 1s^{2}2p^{2} configurations.

Numbers in parentheses are power of 10.

1s | 2s | ^{2}p_{1/2} | ^{2}p_{3/2} |
||

< r^{-3} > | (a) | n.a. | n.a. | n.a. | 3.404(4) |

(b) | n.a. | n.a. | n.a. | 3.582(4) | |

< r^{-1} > | (a) | 1.230(2) | 3.205(1) | 3.187(1) | 2.349(1) |

(b) | 1.227(2) | 3.202(1) | 3.175(1) | 2.374(1) | |

< r > | (a) | 1.356(-2) | 5.455(-2) | 4.406(-2) | 5.226(-2) |

(b) | 1.359(-2) | 5.462(-2) | 4.419(-2) | 5.185(-2) | |

< r^{2} > | (a) | 2.576(-4) | 3.580(-3) | 2.441(-3) | 3.220(-3) |

(b) | 2.587(-4) | 3.589(-3) | 2.455(-3) | 3.176(-3) | |

(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.

Appendix

Eigenstates of the total angular momentum

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^{2} 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^{2}
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^{2} 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) |

[1] E.E. Salpeter and H.A. Bethe. Phys. Rev., 84, 1232, 1951.

[2] J. Sucher. Phys. Rev., 103, 468, 1956.

[3] A.I. Akhiezer and V.B. Berestetskii. Quantum Electrodynamics. John Whiley & Sons Inc., New York, 1965.

[4] J.P. Desclaux. Comp. Phys. Commun., 9, 31, 1975.

[5] S. Love. Ann. Phys., NY, 113, 153, 1978.

[6] J. Sucher. J. Phys. B: At. Mol. Phys., 21, L585, 1988.

[7] I. Lindgren. J. Phys. B: At. Mol. Phys., 23, 1085, 1990.

[8] J. Hata and I.P. Grant. J. Phys. B: At. Mol. Phys., 17, L107, 1984.

[9] P. Indelicato and J.P. Desclaux. Phys. Scripta, to be published, 1992.

[10] C. Froese-Fischer. The Hartree-Fock Method for Atoms. John Wiley & Sons Inc., New York, 1977.

[11] D.R. Hartree. The Calculations of Atomic Structures. John Wiley, New York, 1957.

[12] F.B. Hildebrand. Introduction to Numerical Analysis. McGraw Hill, New York, 1956.

[13] P.J. Mohr. Phys. Rev. A, 46, 4421, 1992.

[14] N.J. Snyderman. Ann. Phys., NY, 211, 43, 1991.

[15] S.A. Blundell. Phys. Rev. A, 46, 3762, 1992.

[16] P.J. Mohr. Ann. Phys., NY, 88, 26, 1974.

[17] P.J. Mohr. Phys. Rev. A, 26, 2338, 1982.

[18] P.J. Mohr and Y.K. Kim. Phys. Rev. A, 45, 2727, 1992.

[19] E.H. Wichmann and N.M. Kroll. Phys. Rev., 101, 843, 1956.

[20] E.A. Uehling. Phys. Rev., 48, 55, 1935.

[21] G. Kallen and A. Sabry. K. Danske Vidensk. Selsk. Mat.-Fis. Medd., 29, 17, 1955.

[22] J. Blomqvist. Nuc. Phys. B, 48, 95, 1972.

[23] L.W. Fullerton, and G.A. Rinker,Jr. Phys. Rev. A, 13, 1283, 1976.

[24] K.N. Huang. Phys. Rev. A, 14, 1311, 1976.

[25] P.J. Mohr. Phys. Rev. Lett., 34, 1050, 1975.

[26] P. Indelicato, O. Gorceix and J.P. Desclaux. J. Phys. B: At. Mol. Phys., 20, 651, 1988.

[27] T.A. Welton. Phys. Rev., 74, 1157, 1948.

[28] P.K. Kabir and E.E. Salpeter. Phys. Rev., 108, 1256, 1957.

[29] R. De Serio, H.G. Berry, R.L. Brooks, J. Hardis, A.E. Livingston and S.J. Hinterlong. Phys. Rev. A, 24, 1872, 1981.

[30] W.R. Johnson and G. Soff. At. Dat. Nuc. Dat. Tables, 33, 405, 1985.

[31] D.F. Mayers. J. de Physique, 11-12, C4-221, 1970.

[32] J.P. Desclaux, C.M. Moser and G. Verhaegen. J. Phys. B: At. Mol. Phys., 4, 296, 1971.

[33] G. Bessis N. Bessis and J.P. Desclaux. J. de Physique, 11-12, C4-231, 1970.

[34] W. Eissner and H. Nussbaumer. J. Phys. B: At. Mol. Phys., 2, 1028, 1969.