Trond Saue
email: trond.saue@irsamc.ups-tlse.fr

Relativistic Hamiltonians for chemistry: a primer

A major breakthrough in relativistic quantum chemistry in recent years is the development of the exact 2-component (X2C) Hamiltonian which reproduces exactly the positive-energy spectrum of the parent 4-component Hamiltonian, yet is constructed in a simple manner using matrix algebra. This mini-review provides an overview of 4-and 2-component relativistic Hamiltonians employed in molecular electronic structure calculations and the underlying physics they carry, thus setting the recent developments in perspective.

INTRODUCTION

Relativity has added a new dimension to chemistry. Pople introduced the concept of theoretical model chemistries [START_REF] Hehre | Ab Initio Molecular Orbital Theory[END_REF][START_REF] Pople | [END_REF] as specic combinations of molecular electronic structure methods and basis sets. Their formal independence means that the space of theoretical model chemistries may be spanned by two orthogonal axis representing the choice of method and basis. In practice balanced combinations of these two are sought, preferably with the possibility of a systematic approach towards exact solutions of the underlying wave equation. Increased accuracy of a model chemistry generally comes with increased computational cost, which can be expressed as N y , where N is the basis set size and y a method-dependent power. A major undertaking in the past few decades has been to bring the computational cost down towards linear scaling with respect to system size [START_REF]Linear-Scaling Techniques in Computational Chemistry and Physics[END_REF]. However, it has also become increasingly clear that the applicability of theoretical model chemistries for the study of chemistry involving elements across the entire periodic table requires that relativistic eects are taken into account [START_REF] Pyykkö | [END_REF]. Such eects arise from the high speed of electrons in the vicinity of heavy nuclei. One distinguishes between scalar relativistic eects, associated with the relativistic mass increase of electrons, and the spin-orbit interaction, generated by magnetic induction, a mechanism that can not be described within a non-relativistic framework [5]. the market, but they can all be dened with respect to the generic form of the electronic Hamiltonian within the Born-Oppenheimer approximation

H = i ĥ (i) + 1 2 i =j ĝ (i, j) + V N N , (1)
where V N N is the classical repulsion of clamped nuclei. The basic formulas of present molecular electronic structure methods can be developed without specication of the one-electron operator ĥ and the two-electron operator ĝ, which implies that the various non-relativistic methods available today can be carried over into the relativistic domain, albeit with possible adaptions. The choice of Hamiltonian introduces accordingly a third axis for the specication of theoretical model chemistries [6], as depicted in gure 1. The independence of the choice of method and of Hamiltonians is often muddled in relativistic terminology where for instance Dirac-Hartree-Fock is presented as a method, whereas it refers to a specic choice of both method and Hamiltonian. The shortened form Dirac-Fock is in addition unfair because Bertha Swirles made the rst formulation of 4-component relativistic Hartree-Fock theory following a suggestion by Douglas Hartree [7].

A convenient specication of Hamiltonians along this axis is with respect to the dimension of the one-electron operator. The non-relativistic one-electron operator is a scalar operator, whereas the fully relativistic Dirac Hamiltonian is a 4 × 4 matrix operator. The increased dimension of the corresponding one-electron wave functions (orbitals) is due to the fact that the Dirac equation explicitly includes spin and describes two kinds of particles: the electron and its anti-particle, the positron, having an entangled existence.

In between the scalar non-relativistic and 4-component relativistic one-electron Hamiltonians there are a number of 2-component relativistic Hamiltonians where the positronic degrees of freedom have been frozen. Whereas scalar relativistic eects may be introduced into a non-relativistic computer code with essentially no extra computational cost, the description of the spin-orbit interaction requires at least a 2-component formalism and will increase computational cost, typically by one order of magnitude, due to the transition from real to complex algebra and the general reduction of symmetry. Full relaxation of the electronic wave function is only possible at the 4-component level, as will be discussed in section II C, but formally increases the computational cost by another order of magnitude. It is, however, important to stress that the additional computational eort associated with the introduction of relativistic

Hamiltonians is a prefactor and thus independent of system size.

In the past few years there has been signicant excitement in the domain of (relativistic) quantum chemistry due to the development of the so-called exact 2-component relativistic (X2C) Hamiltonian which reproduces exactly the positive-energy spectrum of the parent 4-component Hamiltonian. Recently Wenjian Liu, one of the key players in the eld, published a comprehensive review on 2-component Hamiltonians [8]. The review is rather complete and full of insight, but also quite technical and easily gives the impression that the construction of the X2C Hamiltonian is complicated. This has motivated me to write the present mini-review which is a more introductory text. Since I strongly believe that it is important to add some esh in terms of physics to the skeleton of theoretical expressions I start with a discussion of 4-component Hamiltonians. I also hope to clear up some misunderstandings in the literature and share some insight gained over more than twenty years of working with relativistic molecular quantum mechanics. I will furthermore discuss some 2-component approximate relativistic Hamiltonians: i) the Pauli Hamiltonian, because its underlying physics is particularly transparent, and ii) the Zeroth-Order Regular Approximation (ZORA) and Douglas-Kroll-Hess (DKH) Hamiltonians, because they are widely used. I then turn to a discussion of the X2C

Hamiltonian, where I will stress the simplicity of its formulation. I hope to convince the reader of this by describing in some detail the one-step construction of the X2C Hamiltonian as implemented by Miroslav Ilia² in the dirac code [START_REF] Dirac | a relativistic ab initio electronic structure program, Release DIRAC10[END_REF] and recently re-written in more modular form by Stefan Knecht. It should not be forgotten that the 2-component Hamiltonians are typically obtained by some decoupling transformation of a parent one-electron Hamiltonian, and so the proper handling of the two-electron interaction as well as property operators merits separate discussions. Next, The present contribution is a primer and accordingly willingly brief and limited in scope. I hope to leave sucient references for the interested reader to continue exploring this fascinating domain of theoretical chemistry. There are also some excellent textbooks available. In the molecular domain I would recommend the book by Kenneth Dyall and Knut Faegri [START_REF] Dyall | Introduction to Relativistic Quantum Chemistry[END_REF] as well as the book by Markus Reiher and Alexander Wolf [START_REF] Reiher | Relativistic Quantum Chemistry: The Fundamental Theory of Molecular Science[END_REF], which are to some extent complementary. Peter Schwerdtfeger has recently edited a collection of contributions on the fundamentals and applications of relativistic electronic structure theory [START_REF] Schwerdtfeger | Relativistic Electronic Structure Theory. Part 1. Fundamentals[END_REF][START_REF] Schwerdtfeger | Relativistic Electronic Structure Theory. Part 2. Applications[END_REF]. I can also recommend similar compilations by M.

I
Barysz and Y. Ishikawa [START_REF]Relativistic Methods for Chemists[END_REF] as well as U. Kaldor and S. Wilson [START_REF]Theoretical Chemistry and Physics of Heavy and Superheavy Elements[END_REF]. A clear omission in the present mini-review are relativistic pseudopotentials, for which excellent reviews are available [1619].

Unless otherwise stated I employ SI-based atomic units [START_REF] Whien | [END_REF] in the following.

II. 4-COMPONENT HAMILTONIANS

A.

One-electron part

We start from the classical expression for the energy of a relativistic free particle

E 2 = m 2 c 4 + c 2 p 2 (2)
where m is mass, c the speed of light and p linear momentum. Taking the square root on each side of the above expression we are faced with a choice of sign

E = ± m 2 c 4 + c 2 p 2 ; ⇒ E ∈ -∞, -mc 2 ∪ +mc 2 , ∞ (3)
showing that possible values of the energy E lie in two bands of opposite sign separated by a huge energy gap of 2mc 2 . In classical mechanics we can safely discard the negative energy branch as unphysical since the energy can only change in a continuous manner. This is, however, not possible in quantum mechanics, where quantum leaps in energy are, as the name implies, possible.

The relativistic energy expression does not resemble its non-relativistic counterpart

E N R = p 2 2m (4)
at all, but the two can be connected by a Taylor expansion of the positive-energy branch of the former

E = mc 2 1 + p 2 m 2 c 2 = mc 2 rest mass + p 2 2m - p 4 8m 3 c 2 + . . .

kinetic energy (5)

The rst term is the rest mass, the second term is the energy of a non-relativistic free particle, whereas further terms are relativistic corrections that go to zero as c → ∞. Since there is no absolute origin of the energy scale, we can align the relativistic scale with the non-relativistic one by subtracting the rest mass E → E + = E -mc 2 . For negative-energy solutions, on the other hand, the proper change of energy scale in order to obtain nite energies in the non-relativistic limit is

E → E -= E + mc 2 .
The quantization of the non-relativistic energy expression (4) is straightforward: We replace the coordinate and momentum variables by the corresponding operators, perform the heuristic substitution E → i ∂ ∂t and insert a wave function on both sides to obtain the Schrödinger equation for a free particle. The relativistic case is quite a bit more complicated and was not handled successfully until Dirac obtained a linearisation of the energy expression (2) [21,22].

Here show an alternative derivation due to van Waerden [START_REF] Van Der Waerden | Die Grundlehren der mathematischen Wissenschaften Band 37, Die Gruppentheoretische Methode in der Quantenmechanik[END_REF]. It is a useful exercise since it provides an introduction to the use of the Dirac identity

(σ • A) (σ • B) = A • B + iσ • (A × B) . ( 6
)
involving the three Pauli spin matrices collected in the vector σ. A particular instance of this relation is

(σ • p) (σ • p) = p 2 (7)
suggesting that spin may be hidden in the non-relativistic wave equation and only comes into play when an external magnetic eld is introduced. In the following we shall drop the hats indicating operators.

A rst step towards the Dirac equation is to re-arrange the relativistic energy expression (2) and to perform quantization as discussed above

- 1 c 2 ∂ 2 ∂t 2 -p 2 φ 1 = i c ∂ ∂t + (σ • p) i c ∂ ∂t -(σ • p) φ 1 mcφ2 = (mc) 2 φ 1 . (8)
Factorization is obtained through the use of the Dirac relation ( 6), more precisely (7), which implies that the wave function φ 1 is a two-component vector function. We next introduce a second wave function φ 2 to obtain two coupled

equations i c ∂ ∂t -(σ • p) φ 1 = mcφ 2 (a) i c ∂ ∂t + (σ • p) φ 2 = mcφ 1 (b) . (9)
As a third step we take linear combinations

i c ∂ ∂t [φ 1 + φ 2 ] -(σ • p) [φ 1 -φ 2 ] = mc [φ 1 + φ 2 ] (a + b) -i c ∂ ∂t [φ 1 -φ 2 ] + (σ • p) [φ 1 + φ 2 ] = mc [φ 1 -φ 2 ] (b -a) (10)
and introduce the large ψ L and small ψ S components

ψ L = [φ 1 + φ 2 ] ; ψ S = [φ 1 -φ 2 ] . (11)
The resulting equations can be recast on matrix form

i c ∂ ∂t -(σ • p) (σ • p) -i c ∂ ∂t ψ L ψ S = mc ψ L ψ S , (12)
and the conventional form of the Dirac equation

βmc 2 + c (α • p) ψ = i ∂ ∂t ψ; ψ = ψ L ψ S (13)
is obtained by multiplication by βc from the left, where we have introduced the Dirac matrices

α = 0 2 σ σ 0 2 , β = I 2 0 2 0 2 -I 2 . ( 14
)
As in the non-relativistic case we may separate o the time-dependent part to obtain the Dirac equation on timeindependent form, albeit specic for a given inertial frame.

The free-particle Dirac equation is admittedly not very useful for chemistry and we therefore proceed to the introduction of external elds in accordance with the principle of minimal electromagnetic coupling [START_REF] Gell-Mann | [END_REF] p → p -qA; E → E -qφ

where appears particle charge q. The Dirac equation describes both electrons and positrons, having opposite charge. As a rst step towards chemistry we focus on electrons and set q = -e, where e is the fundamental charge. Note, however, that this implies that all solutions, of both positive and negative energy, of the resulting Dirac equation for the given potentials are electronic. The positronic solutions of the same external elds are related to the electronic ones through charge conjugation. For use in the calculation of the electronic structure of molecules we shall invoke the Born-Oppenheimer (clamped nuclei) approximation. We furthermore limit attention to nuclear charge, introducing the notation V = -eφ nuc , where φ nuc is the electrostatic potential of clamped nuclei. The eect of nuclear spin is much weaker and can be treated by perturbation theory, as in the calculation of NMR parameters. We accordingly arrive at the Dirac equation for an electron in a molecular eld

h D ψ = E + ψ; h D = β mc 2 + c (α • p) + V = V c (σ • p) c (σ • p) V -2mc 2 ; β = β -I 4 (16)
where the relativistic energy scale has been aligned with the non-relativistic one, as discussed above.

The so-called minimal substitution ( 15) is employed both in the relativistic and non-relativistic domain, but it is important to realize that it corresponds to a relativistic coupling of particles and elds [5,25]. Non-relativistic theory is notably unable to describe magnetic induction, which is a relativistic eect and forms the basis for the spin-orbit interaction. Figure 2 illustrates the basic mechanism of spin-orbit interaction generated by a clamped nucleus of charge +Ze. In the frame of the nucleus, which is the frame employed for the formulation of the electronic Hamiltonian, the nucleus is a source of a electrostatic potential φ only and the vector potential A is zero. However, carrying out a Lorentz transformation to the frame of the moving electron a non-zero vector potential appears, meaning that the electron in its own frame sees both an electric and magnetic eld due to the nucleus in relative motion. Spin-orbit coupling is the interaction of the electron spin with the magnetic eld induced by a charge in relative motion. A subtle feature of the spin-orbit interaction arises from the fact that the motion of the electron is not of constant velocity and thus has to be followed by a succession of Lorentz transformations. The product of two pure Lorentz transformations (boosts) is not a third boost, rather a boost combined with a rotation [26], so that the succession of boosts induces a rotational motion, denoted the Thomas precession, which reduces the spin-orbit interaction by a factor two.

B. Two-electron part

The derivation of the relativistic expression for the two-electron interaction can be seen as a continuation of the preceding discussion in that we replace the potentials of a clamped nucleus by those of a second moving electron, taking care of avoiding double counting of interactions [START_REF] Moss | Advanced Molecular Quantum Mechanics[END_REF]. In Coulomb gauge the general expression for the scalar potential is

φ(r 1 , t) = ˆρ(r 2 , t) r 12 dτ 2 ; r 12 = |r 1 -r 2 | ( 17
)
which corresponds to the instantaneous Coulomb interaction. The complexity of the relativistic two-electron interaction shows up in the expression for the corresponding vector potential

A(r 1 , t) = 4π c 2 ˆj⊥ (r 2 , t r ) r 12 dτ 2 ; t r = t - r 12 c , (18)
where j ⊥ is the solenoidal, or divergence-free, current density (∇ • j ⊥ = 0), compatible with the constraint of a purely solenoidal vector potential (A = A ⊥ ) in Coulomb gauge. Contrary to the expression for the scalar potential we now have retarded time t r , and not instantaneous time t, under the integral sign, reecting the nite speed of the propagation of interactions in the relativistic domain. For classical electrons we accordingly need to know the position r 2 of the second electron at retarded time. In fact, the full history of the interacting particles is required for a complete relativistic description of the two-electron interaction, and no closed expression is available for use in the electronic Hamiltonian. Rather a perturbation expansion of the full two-electron interaction can be used, starting from the instantaneous Coulomb term known from the non-relativistic domain

g C (1, 2) = e 2 r 12 . (19)
We can indeed consider the instantaneous Coulomb interaction as the non-relativistic limit of electrodynamics, as suggested by the observation the the vector potential [START_REF] Seijo | Computational chemistry: Review of current trends[END_REF] goes to zero in the non-relativistic limit. The question of the non-relativistic limit of electrodynamics is to some extent open, since there is no experiment which can probe this limit. It has been suggested that magnetic induction, e.g. the spin-orbit interaction, goes to zero in the non-relativistic limit, but not the magnetic elds themselves [START_REF] Kutzelnigg | Relativistic Electronic Structure Theory Part 1. Fundamentals[END_REF][START_REF] Kutzelnigg | [END_REF]. However, this leads to incoherent theory, as shown by a simple thought experiment: Suppose that we consider a system consisting of a molecule to which we apply an external magnetic eld. If we now extend the system to include the sources of the magnetic eld, then the magnetic eld will vanish in the non-relativistic limit for the extended system, but not for the original one. A more consistent proposal is that all magnetic interactions as well as retardation disappears in the non-relativistic limit [5,[START_REF] Saue | Relativistic Electronic Structure Theory. Part 1. Fundamentals[END_REF]. A non-relativistic framework can certainly accommodate magnetic elds, but this amounts to a non-relativistic description of particles and a relativistic description of their coupling to electromagnetic elds.

We can derive the rst-order relativistic correction to the two-electron interaction in a heuristic manner by starting from the classical Darwin approximation to the relativistic interaction between two moving charges [START_REF] Darwin | [END_REF]

E int = q 1 φ 2 -q 1 v 1 • A 2 ∼ q 1 q 2 r 12 - q 1 v 1 • q 2 v 2 2c 2 r 12 - (q 1 v 1 • r 12 ) (q 2 v 2 • r 12 ) 2c 2 r 3 12 . (20)
The rst term on the far right is, as expected, the Coulomb term, whereas the nal two terms involve the vector potential of the second electron and come in to order O c -2 . For quantization we shall need the relativistic velocity operator which can be obtained from the Heisenberg equation of motion

dr dt = -i [r, h D ] = cα. ( 21
)
The form of this operator is quite dierent from the non-relativistic velocity operator p m and reects Zitterbewegung [32], a highly oscillatory motion superimposed on the average motion of the electron. We obtain the Breit term

g B (1, 2) = -e 2 cα 1 • cα 2 2c 2 r 12 + (cα 1 • r 12 ) (cα 2 • r 12 ) 2c 2 r 3 12 (22)
by the heuristic substitution v i → cα i in [START_REF] Whien | [END_REF], which is indeed how Breit obtained this operator himself [33]. The Breit term can be re-arranged to [34] g where the rst term is known as the Gaunt term [35] and represents current-current interaction, contrary to the Coulomb term which represents charge-charge interaction.

B (1, 2) = g G (1, 2) + g gauge (1, 2) = - ecα 1 • ecα 2 c 2 r 12 - (ecα 1 • ∇ 1 ) (ecα 2 • ∇ 2 ) r 12 2c 2 , (23)
Both terms give rise to spin-orbit interaction. The separation of spin-orbit interaction due to relative motion of electrons into two terms, contrary to the single term associated with electron-nucleus interaction, is due to the choice of reference frame for the specication of the electronic Hamiltonian, as illustrated in gure 3. The second electron contributes only a scalar potential in its own frame, just as the the clamped nuclei in the Born-Oppenheimer frame.

However, Lorentz transforming to the frame of the nuclei a non-zero vector potential appears, which gives rise to the Gaunt (Breit) term. Both the scalar and vector potential of the second electron, as seen in the Born-Oppenheimer frame, contribute to the vector potential in the frame of the reference electron and give rise to the spin-same-orbit (SSO) and spin-other-orbit (SOO) interactions, respectively. The spin-other-orbit interaction arises solely from the Gaunt term and not from the gauge-dependent term g gauge (1, 2) [START_REF] Saue | Principles and Applications of Relativistic Molecular Calculations[END_REF]. From these considerations we see that the spinorbit interaction associated with nuclei is of the spin-same -orbit type. Finally, it should be noted that the Lorentz transformation of the instantaneous Coulomb term to some other frame does in general not give potentials satisfying

Coulomb gauge, but this can be achieved by carrying out the proper gauge transformation in the new frame. 16) is combined with the Coulomb two-electron operator [START_REF] Klobukowski | Computational chemistry: Review of current trends[END_REF]. Although the latter operator has essentially the same form as in the non-relativistic domain, its physical content is dierent in that it includes spin-same orbit (SSO) interaction, as discussed in section II B. The truncation of the full two-electron interaction means that the electronic Hamiltonian is not strictly relativistic, but it is sucient for most chemical purposes [START_REF] Visser | [END_REF]. For accurate studies of molecular spectra including ne structure it is recommended to include spinother orbit interaction through the Gaunt (Breit) term. The Gaunt term also carries the full spin-spin interaction, whereas the gauge-dependent term g gauge (23) must be included for the full orbit-orbit interaction [START_REF] Saue | Principles and Applications of Relativistic Molecular Calculations[END_REF].

C. The electronic Hamiltonian

The rst 4-component relativistic Hartree-Fock calculations based on the Dirac-Coulomb Hamiltonian and using nite basis sets [3840] were awed since the coupling of large and small components were not taken into account [4144]. From the Dirac equation for an electron in a molecular eld [START_REF] Dolg | [END_REF] the exact coupling is found to be

ψ S = Rψ L ; R(E) = 2mc 2 -V + E + -1 c (σ • p) = 1 2mc B(E) (σ • p) ; B(E) = 1 + E + -V 2mc 2 -1 . ( 24
)
For positive-energy solutions the energy-dependent factor B(E) is on the order of unity O c 0 , demonstrating that the small components are on average a factor c smaller than the large components, hence their name. For negative-energy solutions B(E) can be re-expressed as

B(E) = 2mc 2 E --V , (25)
and is of order O c 2 , such that the relative sizes of the large and small components switch. Due to the energy dependence of the exact coupling [START_REF] Gell-Mann | [END_REF] one rather employs its non-relativistic limit

lim c→∞ cψ S = 1 2m (σ • p) ψ L . (26)
Note, however, that this limit is only valid for |E + -V | 2mc 2 which restricts the energy to the positive-energy branch, for which E + ∼ O c 0 , and the potential φ nuc to non-singular ones, the latter obtained by replacing the point charge nuclear model of the non-relativistic domain by a nite charge distribution. A further strong motivation for the introduction of a nite nuclear charge distribution is that a point nucleus induces a singularity in the electronic wave function at the nucleus, and not a cusp as in the non-relativistic domain. Although the singularity is weak in the sense that it still allows normalization of the wave function, it can not be described in a nite basis using conventional Slater-or Gaussian-type functions and serves as a 'black hole' in basis set optimizations in that it can not be saturated by the introduction of tight functions. With the introduction of a nite nuclear charge model, typically a Gaussian, the singularity is replaced by a Gaussian shape at the origin, favoring the use of Gaussian-type basis functions [4547].

The implementation of (26) at the basis set level has been denoted kinetic balance [48], since it provides a correct representation of the kinetic energy operator in the non-relativistic limit. It is, however, important to have sucient exibility in the basis such that the exact coupling can be obtained [49]. The eect of the missing energy-dependent factor B(E) is illustrated in gure 4 which shows the small component radial function R S of the 1s 1/2 orbital of the mercury atom obtained by a numerical atomic 4-component relativistic Hartree-Fock calculation (exact coupling) or generated from the large component radial function R L using restricted kinetic balance. In this particular case B(E) can be expanded as

B (ε 1s ) = 1 + |ε 1s | 2mc 2 - eφ nuc 2mc 2 + . . . (27)
where ε 1s is the 1s 1/2 orbital energy. A nite (Gaussian) nucleus model was used in this calculation such that there is no singularity of the corresponding potential φ nuc at the origin. Yet the potential attains signicant values close to the nucleus and the B(E) factor will accordingly dampen the RKB radial functional in this region, as seen in the gure. However, bearing in mind that the radial expectation value of the 1s 1/2 orbital is 0.0017 a 0 it can be seen that this modication is extremely local, which implies that if the exact coupling can be obtained by a kinetically balanced basis for an isolated atom it will also work when this atom is placed in a molecule [START_REF] Saue | Theoretical Chemistry and Physics of Heavy and Superheavy Elements[END_REF]. Eq. ( 27) furthermore suggests that away from the nucleus, for suciently small φ nuc , the factor B (ε 1s ) will be slightly larger than unity.

In a nite basis calculation the modication of the coupling of large and small components from kinetic to exact balance is carried by the basis set expansion coecients, which commute with the operator (σ • p). It has therefore been pointed out [5254] that imposing kinetic balance is equivalent to transforming the Dirac equation ( 16) according to

ψ L ψ S = I 2 0 2 0 2 1 2mc (σ • p) ψ L φ L (28)
where φ L is denoted the pseudo-large component. The resulting equation

V T T W 0 -T ψ L φ L = I 2 0 2 0 2 1 2mc 2 T ψ L φ L E; T = p 2 2m ; W 0 = 1 4m 2 c 2 (σ • p) V (σ • p) (29)
was later revived and denoted the modied Dirac equation by Dyall [53], who used it to obtain a spin-orbit free form of the Dirac equation. This is accomplished by using the Dirac identity (6)

(σ • p) V (σ • p) = pV • p + iσ • (pV × p) (30)
and remove the second term. In a quaternion formalism [55] this corresponds to the deletion of imaginary parts of the quaternion modied Dirac equation. Visscher and van Lenthe suggested that the separation of spin-free and spindependent terms is not unique [56]. This is only partially true. Going back to gure 2 we see that a precise mechanism for the elimination of the spin-orbit interaction would be to eliminate the vector potential of the nucleus in the rest frame of the electron. The transformation back to the nuclear frame is equivalent to the forward transformation with the opposite sign of the velocity and would therefore result in the introduction of a compensating magnetic eld in the nuclear frame. However, there is no practical realization of such a procedure.

Another interesting feature of the modied Dirac equation is that, contrary to the parent equation ( 16), it contains only inverse powers of the speed of light such that the non-relativistic limit can be taken without the operator blowing up. The non-relativistic limit of the Dirac equation can also be obtained by the simple metric transformation

ψ L ψ S = I 2 0 2 0 2 I 2 c -1 ψ L ψS (31)
which leads to the equation

V (σ • p) (σ • p) -2m 1 -V 2mc 2 ψ L ψ S = I 2 0 2 0 2 c -2 I 2 ψ L ψ S E, (32)
for which the non-relativistic limit is known as the Lévy-Leblond equation [57]. Note, however, that in both cases it is assumed that E 2mc 2 and V 2mc 2 , which is only true for positive-energy solutions and extended nuclei, respectively. Equation ( 32) is the starting point for direct perturbation theory (DPT) [START_REF] Kutzelnigg | Relativistic Electronic Structure Theory Part 1. Fundamentals[END_REF]5860] which has recently been revived by Stopkowich and co-workers [6163].

The Dirac-Coulomb Hamiltonian has no bound solutions. Expanding the many-electron wave function in an Nparticle basis of Slater determinants one nds that a determinant of bound orbitals will be degenerate with an innite number of determinants containing orbitals from both the positive-and negative-energy continuum. These determinants will mix in and cause a continuum dissolution [64], referred to as the Brown-Ravenhall disease. It can be avoided by embedding the Dirac-Coulomb Hamiltonian by projection operators removing any Slater determinant containing negative-energy orbitals [65]. Such a no-pair Hamiltonian is simply obtained by expressing the electronic Hamiltonian in second quantization

Ĥ = pq h pq a † p a q + 1 2 pqrs g pqrs a † p a † r a s a q (33)
and removing all negative-energy orbitals from the summations. A crucial point and a source of much confusion in the literature is that the separation of the space spanned separately by positive-and negative-energy orbitals depends on the actual potential such that projectors can not be introduced without carefully specifying the potential on which they are based.

The one-electron basis for which the N-particle basis is constructed is generated by solving some eective oneelectron Dirac equation, typically the Hartree-Fock equation. Such an equation has, just as the Dirac equation itself, solutions of both positive and negative energy. The bound orbitals employed for the construction of the mean-eld potential are not selected according to an Aufbau principle, starting from the orbitals of lowest energy. Rather one starts from the lowest positive-energy orbital, which is readily identied due to the large energy gap, on the order of 2mc 2 , down to the negative-energy solutions. This procedure corresponds to the implicit use of projection operators updated in every iteration of the SCF procedure, as suggested by Mittleman [66]. It furthermore corresponds to the variational solution of the Dirac equation, replacing the minimization principle of the non-relativistic domain by a min-max principle [67,68]. The original proposal by Sucher [65] was to employ projection operators constructed from the solutions of the free-particle Dirac equation. However, as pointed out by Heully and co-workers [69], this will induce serious errors. A better approximation to the mean-eld projectors is to use the bare nucleus projectors, that is, using the solutions of the Dirac equation in the molecular eld [START_REF] Dolg | [END_REF], which, as we shall see in section III, are implicitly used in most schemes for generating 2-component relativistic Hamiltonians.

Starting from the Hartree-Fock solutions, electron correlation may be added through Conguration Interaction (CI) or Coupled Cluster (CC) procedures. The N-particle determinants are usually constructed from positive-energy orbitals only. This bars, however, the complete relaxation of the one-particle basis and thus the projectors to the full instantaneous potential seen by the electrons of the system and implies that a full CI in this determinantal space can not be considered as the exact solution of the optimal projected Dirac-Coulomb Hamiltonian. Various authors [7072] have reported CI calculations using an N-particle basis constructed from both positive-and negative-energy orbitals.

Pestka and co-workers [73] have reported a procedure in which the electronic ground state is treated as a resonance and extracted from the continuum by a complex-coordinate rotation technique. However, in all such procedures going beyond the no-pair approximation the electronic ground state is therefore at best treated as a metastable state.

Ultimately this would imply that matter is not stable, which is physically not very satisfying. In the early days of the Dirac equation the stability of matter caused considerable worry amongst theoreticians since calculations based on the Dirac equation suggested that the hydrogen atom would have a lifetime of about a nanosecond, after which the electron would descend down the negative-energy continuum, liberating an innite amount of energy. Dirac nally proposed that all negative-energy orbitals are occupied and thus not available. A more satisfying treatment of the negative-energy states is obtained in the framework of QED which treat electrons and positrons on an equal footing.

Present-day relativistic molecular electronic structure calculations treat the negative-energy orbitals as an orthogonal complement which is continuously updated at the SCF level and frozen at the correlated level. A complete relaxation of the electronic wave function is possible at the MCSCF level with the N-particle basis constructed according to the no-pair approximation [START_REF] Saue | Theoretical Chemistry and Physics of Heavy and Superheavy Elements[END_REF]. In the limit of a full CI rotations between occupied and virtual positiveenergy orbitals become redundant. This is not the case for rotations between occupied positive-energy orbitals and virtual negative-energy orbitals which assure complete update of the projection operators in which the Dirac-Coulomb Hamiltonian is embedded.

III. 2-COMPONENT HAMILTONIANS

A.

General considerations

We now turn to a discussion of 2-component relativistic Hamiltonians. The basic idea is to generate a 2-component

Hamiltonian which reproduces the positive-energy spectrum of the parent 4-component Hamiltonian. Foldy and

Wouthuysen [74] proposed to look for a unitary transformation U that will formally decouple the large and small components

U † h LL h LS h SL h SS U = h++ 0 0 h-- (34)
Another approach that has been discussed in the literature is the elimination of the small component (ESC). The two approaches can be shown to be equivalent [75] and we shall focus on the latter. Let us, however, begin with the elimination of the small components. We will write the 4-component relativistic eective one-electron equation, e.g.

the Dirac equation, as a system of coupled equations

h LL ψ L + h LS ψ S = Eψ L (35)
h SL ψ L + h SS ψ S = Eψ S (36)
Using the exact coupling R of the large and small component [START_REF] Gell-Mann | [END_REF] to eliminate ψ S in (35) we obtain a 2-component equation for the large component only

[h LL + h LS R] ψ L = Eψ L . (37)
However, the large component is not normalized to unity and for this purpose we seek a normalization operator N + . We set φ = N + ψ L and nd

φ | φ = ψ | ψ = ψ L | ψ L + ψ S | ψ S = ψ L 1 + R † R ψ L (38)
from which we deduce that the normalization operator is

N + = √ 1 + R † R.
h ++ φ = Eφ; h ++ = N + [h 11 + h 12 R] N -1 + = 1 + R † R [h 11 + h 12 R] 1 √ 1 + R † R (39)
An alternative route to the elimination of the small component was proposed by Dyall [76]: Summing the equations ( 35) and ( 36), the latter multiplied by R † from the left, gives

h LL + h LS R + R † h SL + R † h SS R ψ L = E 1 + R † R ψ L . ( 40
)
Comparing with (38) one sees that a metric has been introduced on the right hand side of ( 40) which provides correct normalization of ψ L . Equations ( 37) and ( 40) have therefore been dubbed unnormalized and normalized elimination of the small components (UESC and NESC), respectively.

Due to the explicit energy-dependence appearing in the exact coupling operator R in Eq. ( 24) early 2-component relativistic Hamiltonians were based on approximate couplings, as will be discussed in section III B. However, in a very insightful paper Heully and co-workers [77] demonstrated that it is possible to derive an expression for the coupling operator R without any explicit energy-dependence. Multiplying through (35) with R from the left it is possible to equate the left-sides of ( 35) and ( 36)

Rh LL ψ L + Rh LS Rψ L = Eψ S = h SL ψ L + h SS Rψ L . ( 41
)
Since the large component ψ L is now arbitrary we arrive at an operator equation for the exact coupling R without any explicit reference to energy

h SL + h SS R = Rh LL + Rh LS R. ( 42
)
This operator equation is second order in R and has accordingly two sets of solutions, R + and R -, corresponding to the positive-and negative-energy branch, respectively, and being of order O(c -1 ) and O(c 1 ). By an analogous procedure it is also possible to arrive at an operator equation for the inverse coupling ψ L = R -1 ψ S which takes the form

h LS + R -1 h SS = -h LL R -1 + R -1 h SL R -1 . ( 43
)
Taking the Hermitian conjugate of the coupling equation shows that R -1 and -R † fulll the same equation. Combined with order analysis Heully et al. [77] were able to make the identications

R -1 -= -R † + ; R -1 + = -R † - ( 44
)
and to provide a general form for the Foldy-Wouthuysen transformation (34)

U = Ω + -R † Ω - RΩ + Ω - ; Ω + = 1 √ 1+R † R = N -1 + Ω -= 1 √ 1+RR † (45)
It is convenient [START_REF] Kutzelnigg | Relativistic Electronic Structure Theory Part 1. Fundamentals[END_REF]78] to split the transformation into two parts U = W 1 W 2 . The rst transformation

W 1 = 1 -R † R 1 (46)
decouples the large and small components whereas the second transformation

W 2 = Ω + 0 0 Ω - ; Ω + = 1 √ 1+R † R Ω -= 1 √ 1+RR † (47)
re-establishes normalization.

It is instructive to rst consider the unitary decoupling transformation (45) applied to the Dirac 4-spinor

U † ψ L ψ S = Ω + ψ L + R † ψ S Ω -ψ S -Rψ L (48)
For positive-energy solutions we want the lower components to be zero which leads to the relation

ψ S + = Rψ L + , (49)
indeed showing that R can be identied with the coupling R + between positive-energy large-and small components.

Further manipulation of the positive-energy 2-component solution gives

ψ 2c + = 1 √ 1 + R † R ψ L + R † ψ S = 1 √ 1 + R † R ψ L + R † Rψ L = 1 + R † Rψ L , (50)
suggesting that it can be expanded in the large component basis only.

The decoupling transformation W 1 gives a Hamiltonian on the form

h LL + h LS R + R † h SL + R † h SS R -h LL R † + h LS -R † h SL R † + R † h SS -Rh LL -Rh LS R + h SL + h SS R Rh LL R † -R † h LS -h SL R † + h SS . (51)
It can be seen that the lower and upper o-diagonal blocks are identical to the coupling equation ( 42) and its Hermitian conjugate, respectively, and are accordingly zero for exact coupling R. The second transformation W 2 accomplishes renormalization, thus assuring unitarity of the total transformation, and generates a 2-component relativistic Hamiltonian for positive-energy solutions on the form

h ++= 1 √ 1 + R † R h 11 + h 12 R + R † (h 21 + h 22 R) 1 √ 1 + R † R . ( 52
)
In the case of exact coupling we can use the coupling equation [START_REF] Mclean | Current Aspects of Quantum Chemistry[END_REF] to simplify the Hamiltonian further to

h ++ = 1 + R † R [h 11 + h 12 R] 1 √ 1 + R † R . (53)
The nal 2-component equation now reads

h ++ ψ 2c + = Eψ. 2c + (54)
The 2-component Hamiltonian ( 53) is strictly identical to the one obtained by UESC (39) and shows the equivalence of the two decoupling approaches, but only for exact coupling. On the other hand, using ( 52) and ( 50) the 2-component equation ( 54) is straightforwardly re-arranged to the NESC equation (40), showing their equivalence also in the case of approximate coupling. NESC can indeed be shown to provide more accurate energies in the approximate coupling case and is as such to be preferred over UESC [76,[START_REF] Dyall | Theoretical Chemistry and Physics of Heavy and Superheavy Elements[END_REF].

B. Approximative 2-component relativistic Hamiltonians

A rst generation of 2-component relativistic Hamiltonians were developed based on approximations to the exact energy-dependent coupling [START_REF] Gell-Mann | [END_REF] of the large and small components and have found their way into various widely distributed computer codes. In this section I give a brief overview of the most important Hamiltonians of this kind.

The section is intentionally brief since our main focus is on exact 2-component relativistic Hamiltonians and since the approximate variants are amply described elsewhere [START_REF] Dyall | Introduction to Relativistic Quantum Chemistry[END_REF][START_REF] Reiher | Relativistic Quantum Chemistry: The Fundamental Theory of Molecular Science[END_REF].

(i) = V + T i - p 4 i 8m 3 c 2 mass-velocity + 1 4m 2 c 2 σ • [(∇V ) × p i ] spin-orbit + 1 8m 2 c 2 ∇ 2 V Darwin (55)
is based on the approximate coupling

R = 1 2mc 1 + E + -V 2mc 2 -1 (σ • p) ∼ 1 2mc (σ • p) (56)
and is obtained by applying the unitary decoupling transformation (45) and retaining terms only to O(c -2 ), which adds three relativistic correction terms to the non-relativistic one-electron Hamiltonian. The mass-velocity term is easily recognized as a rst-order correction to the kinetic energy by comparing with the relativistic free-particle energy (5). The second term is the generic form of the one-electron spin-orbit operator and is seen to not contain the orbital angular momentum operator. Upon insertion of the electrostatic potential of clamped point nuclei we obtain the more familiar form

V = - Z A e 2 r iA → h SO = A Ze 2 4m 2 c 2 r 3 iA σ • l Ai . (57)
However, it should be noted that the orbital angular momentum operators l Ai are dened with respect to each nucleus A and represents the relative motion of electron and nuclei.

The third term, the Darwin term, represents a correction to the electron-nucleus interaction, V = -eφ(r) → -eφ(r + δ), (58) due to Zitterbewegung. The Darwin term may be generated by Taylor-expanding the modied interaction to second order in the Zitterbewegung amplitude δ and carrying out suitable time averaging (indicated by brackets . . . T ) [START_REF] Moss | Advanced Molecular Quantum Mechanics[END_REF] V + h Darwin = -eφ(r) -e (δ • ∇) T φ(r) -

1 2 e (δ • ∇) 2 T φ(r) = -eφ(r) -e δ 2 T 6 ∇ 2 φ(r) (59)
Assuming isotropicity of the Zitterbewegung the rst-order term averages to zero and only diagonal terms δ 2

x T = δ 2 x T = δ 2

x T = 1 3 δ 2 T are retained in the second-order term. The precise form of the Darwin term is then obtained through the identication δ x = δ y = δ z = 1 2mc . Interestingly, this amplitude can be obtained from a hand-waving model of Zitterbewegung depicted in gure 5: In the vicinity of an electron its eld is suciently strong to allow the creation of a electron-positron pair. The positron annihilates the original electron and the new electron takes over.

The creation of an electron-positron pair requires an energy of at least 2mc 2 . Using the energy-time uncertainty relation we deduce that the maximum particle displacement is

∆E∆t ≥ 1 → ∆x = c∆t ≈ 1 2mc (60)
which corresponds exactly to the Zitterbewegung amplitude given above. This simple model accordingly suggests that the Dirac equation does not describe a single electron, rather a relay of electrons ! A nice feature of the Pauli Hamiltonian is that its underlying physics is reasonably transparent, yet it has limited applicability in relativistic molecular electronic structure calculations. One disadvantage is that the mass-velocity term has no lower bound, in contrast to the square-root expression of the relativistic free-particle energy (5) from which it is extracted. This feature bars the use of the Pauli Hamiltonian in variational calculations. The Pauli Hamiltonian has on the other hand been successfully employed to provide rst-order relativistic corrections in high-precision work on interaction potentials of light atoms such as helium (see for instance [START_REF] Przybytek | [END_REF]). A second disadvantage of the Pauli Hamiltonian is that the Darwin term and higher corrections contain derivatives of the nuclear potential and are thus highly singular terms which are dicult to describe in a nite basis approximation.

The regular approximation alleviates these complications by using a dierent approximation to the exact coupling

(24) R = c 2mc 2 -V 1 + E + 2mc 2 -V -1 (σ • p) ∼ c 2mc 2 -V (σ • p) . (61)
ĤZORA = V + T ZORA ; T ZORA = 1 2m (σ • p) K (r) (σ • p) ; K (r) = 2mc 2 2mc 2 -V (62)
contains an eective kinetic energy operator T ZORA which goes into the non-relativistic one when V → 0, in practice away from nuclei. As depicted in gure 6, for the mercury atom the relativistic factor K (r) is essentially unity at a distance of 0.2 a 0 , just outside the radial expectation value of the 3s 1/2 orbital. Expanding K (r) in orders of the potential and comparing with the Pauli Hamiltonian shows that the ZORA Hamiltonian contains no mass-velocity term, only parts of the Darwin term and the full one-electron spin-orbit term [START_REF] Dyall | Introduction to Relativistic Quantum Chemistry[END_REF]. Spin-dependent and -independent parts of the eective kinetic energy operator T ZORA can be separated using the Dirac identity (6). A nice feature of the ZORA Hamiltonian is that the insertion of the electronic mean-eld or the full Kohn-Sham potential in T ZORA automatically provides the inclusion of two-electron spin-orbit contributions. The potential term in T ZORA leads on the other hand to electrical gauge-dependence in that adding a constant ∆ to the potential, V → V + ∆, does not lead to the addition of the same constant to the energy, E → E + ∆, as is observed in both the Schrödinger and the Dirac equation. This gauge-dependence, which shows up in the calculation of molecular or electronic binding energies, can be curbed by using the so-called scaled ZORA equation [83], based on approximate renormalization, or by replacing the potential term V by a model potential constructed from a superposition of atomic potentials [85] or densities [86].

The direct analytical calculation of T ZORA in a nite basis is rather forbidding due to the presence of a potential term V in the denominator and is better handled by 3D numerical integration, which favors implementation in the framework of density functional theory (DFT) rather than wave function theory (WFT). The rst implementation of ZORA in WFT was nevertheless reported by Faas and co-workers [8789]. A fully analytic implementation of the ZORA Hamiltonian may be obtained by setting the SS block of the metric of the modied Dirac equation ( 29) to zero [55]. An unnormalized elimination of the small components then leads to an eective kinetic energy operator on the form

T ZORA = T (T -W 0 ) -1 T, (63)
to be constructed in matrix representations of T and W 0 . Further manipulation leads to the form

T ZORA = T (T -W 0 ) -1 (T -W 0 + W 0 ) = T + W ; W -1 = W -1 0 -T -1 (64)
employed in the analytic implementation of Filatov and Cremer [90,91] to which later a gauge-correction was added [92].

A third class of 2-component approximate relativistic Hamiltonians have been motivated by the observation that the variational instability of the Pauli Hamiltonian arises from the truncated expansion of the relativistic free-particle energy (5), whereas the complete square-root expression has a lower bound. The idea is therefore to keep the squareroot form for the kinetic energy in the nal 2-component relativistic Hamiltonian. This is achieved by rst carrying

U † 0 ĥU 0 = E p -mc 2 0 0 -E p -mc 2 + A [V + R 0 V R 0 ] A A [R 0 , V ] A -A [R 0 , V ] A A [V + R 0 V R 0 ] A , (65)
where

E p = m 2 c 4 + c 2 p 2 ; R 0 = c (σ • p) E p + mc 2 ∼ O c -1 ; A ≡ Q + 0 = E p + mc 2 2E p ∼ O c 0 . ( 66
)
The free-particle part now has the desired diagonal form and singularities arising from the potential V are regulated by the kinematical factors (66), but the o-diagonal blocks of the potential part are of order unity in V and c -1 , such that the total Hamiltonian U † 0 ĥU 0 ++ , denoted the rst-order Douglas-Kroll-Hess (DKH1) Hamiltonian and extracted from the (++) block of the transformed Hamiltonian, has only limited applicability in the relativistic domain [93,94]. Further decoupling is therefore necessary. Douglas and Kroll proposed decoupling through a sequence of unitary transformations

U = U 0 U 1 U 2 . . . ( 67
)
in orders of the external potential V [95]. Bernd Hess devised matrix techniques for the construction of the kinematical factors (66), starting from a matrix representation of the non-relativistic kinetic energy operator in uncontracted nite standard basis sets of quantum chemistry, and thereby opened up for widespread use of this approach in quantum chemistry. His contribution can not be underestimated, and this approach is therefore justly referred to as the Douglas-Kroll-Hess transformation. The exact decoupling equation ( 24) for the free-particle transformed Hamiltonian (65) can be expressed as

[R 1 , E p ] + = -A [R 0 , V ] A O(V 1 ) + [A [V + R 0 V R 0 ] A, R 1 ] O(V 2 ) -R 1 A [R 0 , V ] AR 1 O(V 3 ) (68)
From order analysis in terms of the external potential the coupling is approximated as

[R 1 , E p ] + ∼ -A [R 0 , V ] A, (69)
dening the approximate decoupling transformation U 1 which generates the second-order Douglas-Kroll-Hess Hamiltonian (DKH2). The original implementation by Hess stopped at this order, but the method has been extended to higher [9698] and even arbitrary order [99,100].

Barysz, Sadlej and Snijders carried out an order analysis of the exact coupling equation ( 68) in terms of c -1 rather than the external potential [101]. The lowest-order approximation to R 1 is then

R 1 ∼ - 1 2mc 2 A [R, V ] A O(c -3 )
giving a 2-component Hamiltonian correct through order O c -4 , but also containing higher-order terms. Rather than performing a sequence of unitary transformations, as in the DKH-approach, they instead proposed to solve the coupling equation ( 68) in an iterative manner to obtain the coupling R 1 correct through some odd order 2k -1 in c -1 and then perform a single unitary transformation, in addition to the free-particle one, to obtain a 2-component Hamiltonian correct through order 2k in c -1 [101103], denoted the BSS h 2k Hamiltonian. The innite-order 2component Hamiltonian generated in the limit k → ∞ were given the acronym IOTC by Kedziera and Barysz [104].

C. Exact 2-component relativistic Hamiltonians

The development of exact 2-component relativistic Hamiltonians (X2C), reproducing exactly the positive-energy spectrum of the parent 4-component Dirac Hamiltonian, arose from the conuence of two important realizations:

1. Obtaining the exact coupling between the large and small components requires the solution, explicit or implicit, of the Dirac equation. However, the problem at hand in a typical quantum chemical calculation is to nd approximate solutions to the many-electron Hamiltonian [START_REF] Hehre | Ab Initio Molecular Orbital Theory[END_REF]. The computational cost of such calculations is completely dominated by quantities associated with the two-electron term. Initially solving the one-electron problem is therefore meaningful and inexpensive.

2. The analytic expression of the exact coupling is not available in closed form. However, starting from a matrix representation of the parent 4-component Hamiltonian the entire decoupling transformation can be carried out using matrix algebra, without ever having to program integrals over additional complicated operator expressions.

The simplest algorithm for the generation of an exact 2-component relativistic Hamiltonian can then be stated in words:

1. Solve the parent 4-component one-electron equation on matrix form.

2. Extract the coupling R (24) from the solutions.

3. Construct the transformation matrix U (45), next h X2C .

The two key realizations leading to the simple construction of the X2C Hamiltonian matured over a number of years in the relativistic quantum chemistry community, but both are present in a 1997 paper by Dyall [76]. Dyall sought to interface relativistic and non-relativistic methods for molecular electronic structure calculations within the nite basis approximation. For this purpose he devised the NESC equation ( 40) which goes smoothly into the Schrödinger equation in the non-relativistic limit. He does not reference the work of Heully et al. [77] and works with an explicitly energy-dependent matrix expression for the coupling R (U in his notation). He does state, though: in principle, a single matrix U can be dened which simultaneously connects all the large and small components for the entire set of solutions, and need not to be given for each solution [76]. As a numerical example he considers a one-electron system, N e +9 , for which he obtains NESC solutions using the energy-dependent expression of U , but also extracting the coupling matrix directly from the corresponding 4-component eigenvectors.

The rst computer implementation that actually delivered the X2C Hamiltonian on matrix form for use in general molecular applications was reported by Hans Jørgen Aagaard Jensen and Miroslav Ilia² and employed in a relativistic Coupled Cluster study of lead oxide [105]. The X2C Hamiltonian was constructed in the spirit of the BSS Hamiltonian, but in a simpler manner. In a rst step the free-particle Dirac equation was solved on matrix form. The matrix representation of the Dirac Hamiltonian in the molecular eld ( 16) was then transformed using the free-particle eigenvectors, a procedure which is equivalent to the free-particle Foldy-Wouthusen transformation (65). In a second step, the Dirac equation, now in free-particle basis, was solved on matrix form, the coupling extracted and a matrix representation of h X2C constructed. The implementation and application of this algorithm was presented by Jensen at the 2005 Conference on Relativistic Eects in Heavy Elements (REHE), but unfortunately not followed up by any publication.

The initial free-particle transformation allowed Jensen and Ilia² to construct various nite-order 2-component Hamiltonians such as the second-order DKH Hamiltonian. However, for the construction of the innite-order 2component Hamiltonian this initial step is not needed. Kutzelnigg and Liu [54,106,107] therefore proposed a one-step procedure, at the time denoted exact quasirelativistic theory (XQR), for the generation of the exact 2-component oneelectron Hamiltonian, starting from a matrix representation of the modied Dirac Hamiltonian (29), and obtaining the coupling either by diagonalisation or by iterative solution of the coupling equation [START_REF] Mclean | Current Aspects of Quantum Chemistry[END_REF]. The scheme was further developed and extended to the Kohn-Sham framework by Liu and Peng [108,109]. A one-step procedure based on diagonalisation of the 4-component parent matrix equation was independently developed by Ilia² and Saue under the acronym IOTC [75]. To avoid a profusion of acronyms for the same object, albeit obtained by somewhat dierent algorithms, it was agreed [START_REF] Jensen | The generic acronym X2C was introduced following discussion at[END_REF] to use the generic acronym X2C for the eXact 2-Component Hamiltonian obtained by matrix algebra.

Let us briey go through the one-step algorithm for the construction of the X2C Hamiltonian as implemented in the dirac program package [START_REF] Dirac | a relativistic ab initio electronic structure program, Release DIRAC10[END_REF]:

1. We start from a matrix representation of the parent 4-component one-electron equation, e.g. the Dirac equation, in a nite non-orthogonal atomic orbital (AO) basis

h 0 c 0 = S 0 c 0 ε. (70)
The four components of the Dirac 4-spinor are individually expanded in Cartesian Gaussian basis functions, the large and small component bases being related by unrestricted kinetic balance. The equation is subsequently transformed to orthonormal basis, using canonical orthonormalisation [START_REF] Löwdin | [END_REF] which provides control of linear dependencies

h 1 c 1 = c 1 ε; h 1 = V † 1 h 0 V 1 ; V † 1 S 0 V 1 = I (71) 2.
A matrix representation W of the transformation (28) to the modied Dirac equation is constructed in the orthonormal basis and subsequently applied

h 2 c 2 = S 2 c 2 ε; S 2 = W † W, (72)
thus enforcing restricted kinetic balance [55]. Orthonormality is restored by a transformation V 2

h 3 c 3 = c 3 ε; V † 2 S 2 V 2 = I, (73)
thus avoiding pitfalls [109] of carrying out the decoupling transformation U in a non-orthogonal basis. The total transformation W V 2 is accordingly unitary.

3. The eigenvalue problem ( 73) is solved and the coupling R extracted. In the matrix representation the coupling relations are given as

C S + = RC L + ; C L -= -R † C S - (74)
where C

L(S) +(-) is the block of the eigenvector matrix corresponding to the large (small) components of the positive(negative)-energy solutions. Rather than solving either equation directly we follow the original implementation of Jensen and Ilia² and construct the equation

AR = B; A = C S -C S † -; B = -C S -C L †
-, where the matrix A is Hermitian and positive-denite such that the system can be solved by a Cholesky decomposition [75].

4. From the coupling R the decoupling transformation matrix U is constructed. Note that the renormalization transformation W 2 (47) can be constructed using canonical orthonormalisation. The matrix representation of h X2C is constructed and subsequently transformed back to the initial AO-basis.

It should be noted that spin-orbit interaction can be eliminated from the 2-component Hamiltonian in two dierent ways, either from the parent 4-component modied Dirac Hamiltonian, or from the nal 2-component Hamiltonian.

The rst route gives an exact 2-component Hamiltonian which reproduces exactly the positive-energy solutions of the 4-component spin-free Hamiltonian, whereas the latter route corresponds to the approach taken to eliminate spin-dependent terms of approximate Hamiltonians such as DKH2.

For large systems the diagonalisation of the Dirac Hamiltonian for the full molecule may become a bottleneck, which can be avoided by extracting the coupling at the level of individual atoms or fragments, as already proposed by Liu and co-workers [108,112].

D. Beyond one-electron Hamiltonians

So far our discussion of 2-component relativistic Hamiltonians has focused rather exclusively on one-electron Hamiltonians. However, construction of the full electronic Hamiltonian (1) also requires nding a suitable form for the two-electron operator at the 2-component level. The most consistent choice would be to transform the 4-component two-electron operator using the same decoupling transformation as for the one-electron part. On matrix form this can be expressed as [113] [U (1

) ⊗ U (2)] † G(1, 2) [U (1) ⊗ U (2)] =     G++ ++ G+- ++ G-+ ++ G-- ++ G++ +- G+- +- G-+ +- G-- +- G++ -+ G+- -+ G-+ -+ G-- -+ G++ -- G+- -- G-+ -- G-- --     ,
where G(1, 2) represents the supermatrix of two-electron integrals at the 4-component level. After transformation the G++ ++ block would be combined with the h++ block from the corresponding transformation (34) of the one-electron

Hamiltonian. However, the transformation involves the full set of two-electron integrals at the 4-component level and the generation of the 2-component two-electron integrals will be more expensive than the corresponding 4-component calculation. Such a scheme has nevertheless been reported by Seino and Hada [114] at the innite-order DKH level, and reductions of computational cost are presently being sought by exploiting locality [START_REF] Seino | Local unitary transformation in two-component relativistic scheme for large-scale molecular systems[END_REF]. A third-order DKH scheme including the mean-eld potential in the decoupling procedure was previously reported by Nakajima and Hirao [START_REF] Nakajima | [END_REF]117].

A straightforward approximation is to combine the 2-component one-electron Hamiltonian with the untransformed two-electron Coulomb operator. However, this translates into a neglect of two-electron spin-spin and spin-orbit contributions, the latter mandatory for a correct description of the ne structure of molecular spectra and for magnetic properties such as the g tensor. Approximate corrections to the untransformed Coulomb operator are most straightforwardly done in a DFT framework since they enter through an eective one-electron potential [86,112,118120].

In a WFT framework two-electron spin-orbit corrections may be introduced by an atomic mean-eld approach [121].

An example is the widely used AMFI code [122] in which the two-electron spin-orbit contributions at the rst-order DKH level are constructed in a mean-eld fashion for each constituent atom of the molecule and then added to the matrix representation of the molecular one-electron Hamiltonian. At the correlated CI or CC level corrections to the untransformed two-electron operator can be introduced by a molecular mean-eld approach [113]. To see how this works, we re-write the second-quantized electronic Hamiltonian (33) on normal-ordered form

Ĥ = pq F pq a † p a q + 1 2 pqrs g pqrs a † p a † r a s a q (75)
where now appears elements F pq of the Fock matrix. An initial SCF calculation is carried out at the 4-component level and the decoupling transformation is then dened and carried out with respect to the converged Fock matrix. This can be justied from the higher computational scaling of the correlated calculation. The 2-component Fock matrix, which now carries innite-order molecular mean-eld spin-orbit corrections is then combined with the untransformed two-electron operator and employed at the correlated level.

An important observation is that errors introduced by an approximate treatment of the two-electron interaction at the 2-component level can be assessed by 4-component calculations [113]. Going back to the unitary decoupling transformation (34) it is evident that an Hamiltonian which reproduces the positive-energy solutions only of the parent 4-component Hamiltonian is simply obtained by full diagonalisation of the parent 4-component Hamiltonian and then backtransformation using only the positive-energy solutions. This is not a very useful approach in the sense that the solutions must still be expanded in bases for both the large and small components, but it tells us that the generation of the 2-component exact Hamiltonian h++ is equivalent to projecting out the negative-energy solutions of the parent 4-

h 2c = U † h 4c U ++ . (76)
For consistency 4-component property operators Ω 4c must be subjected to the same decoupling transformation as the Hamiltonian, that is

Ω 2c = U † Ω 4c U ++ (77)
Use of the approximate expression

Ω 2c ≈ Ω 4c LL (78)
leads to picture change errors that may be larger than the relativistic eects [123126]. A nice feature of the simple algorithm for the construction of the X2C Hamiltonian is that it provides an explicit representation of the transformation matrix U such that property operators can be transformed on the y, thus avoiding picture change errors, as well as the programming of new integrals.

As an example we may consider the electron density in some point P which is formally to be calculated as the expectation value of the corresponding operator. However, at the 4-component SCF level the expectation value reduces to the familiar sum of products of the occupied orbitals

ρ 4c (P ) = -e i ψ 4c i |δ (r -P )| ψ 4c i = -e i ψ 4c † i (P ) ψ 4c i (P ) .
This is not the case at the 2-component level where the picture change transformation leads to a modication (smearing) of the Dirac delta function

ρ 2c (P ) = -e i ψ 2c i U † δ (r -P ) U ++ ψ 2c i = -e i ψ 2c † i (P ) ψ 2c i (P ) .
The dierence is illustrated in gure 7 where we trace the electron number density of the mercury atom, calculated as orbital products, as a function of radial distance at the 4-and 2-component level. Close to the origin the two curves deviate signicantly. The dierence is appreciable, though, only within the radial expectation value r 1s of the 1s orbital and basically disappears on the more chemical scale employed in gure 8, and so one may question its chemical relevance. There are, however, a number of molecular properties which probe the electron density near nuclei, such as NMR parameters, electric eld gradients, molecular gradients as well as Mössbauer isomer shifts, and which provide local information with great sensitivity to the chemical environment. For such properties proper picture change transformation is crucial in order to obtain meaningful results at the 2-component relativistic level, as for instance shown in our recent study of the Mössbauer isomer shifts of mercury uorides [127].

IV. NUMERICAL EXAMPLE

It is beyond the scope of the present mini-review to furnish detailed benchmark studies of the relativistic Hamiltonians described in the text. However, to provide a rst impression of their performance and to illustrate various features discussed in the text I give in Table I orbital energies of the mercury atom extracted from Hartree-Fock calculations using various relativistic Hamiltonians. All results have been obtained with the dirac code [START_REF] Dirac | a relativistic ab initio electronic structure program, Release DIRAC10[END_REF] with identical basis sets, such that all dierences can be unambiguously attributed to the dierence between Hamiltonians. in Table I. It can be seen that with the exception of the 1s 1/2 orbitals the deviation with respect to the reference is quite small. This also holds for the spin-orbit splittings ∆ SO .

We next consider the results obtained with the two lowest-order Douglas-Kroll-Hess Hamiltonians. The one-electron

Hamiltonians have been combined with the untransformed Coulomb term and two-electron spin-orbit corrections from the amfi package [122]. The rst-order Douglas-Kroll-Hess Hamiltonian DKH1 is obtained, as discussed in section III B, by a free-particle Foldy-Wouthuysen transformation. As this transformation provides an exact decoupling of the free-particle Dirac Hamiltonian, the DKH1 Hamiltonian could also be denoted X2C(f ree). The results obtained with the DKH1 Hamiltonian are accordingly comparable to the 4-component DC results using free-particle operators, denoted DC(P f ree ) in Table I, the dierences being due to picture-change errors from the two-electron term. It can indeed be seen that both Hamiltonians, DKH1 and DC(P f ree ), have comparable poor performance and can not be recommended for use in quantum chemical calculations. The second-order Douglas-Kroll-Hess Hamiltonian DKH2 has a signicant better performance, which explains why it has found such widespread use in quantum chemistry codes.

Interestingly, the ZORA Hamiltonian has larger errors in the inner core region than the DC(P f ree ) Hamiltonian, but reproduces the reference DC Hamiltonian quite well for outer core and valence orbitals. The ZORA numbers in

V. SUMMARY AND OUTLOOK

None of the many-electron Hamiltonians discussed in this mini-review are truly relativistic, in the sense of being Lorentz invariant. The time-dependent Schrödinger equation is manifestly non-relativistic since it involves secondderivatives in the spatial coordinates and a rst-derivative in the time coordinate. This is contrary to the theory of special theory in which time and space are treated on an equal footing to within a sign dictated by the metric. As such, the Dirac equation for an electron in a molecular eld ( 16) also appears to violate Lorentz invariance. This equation can, however, be derived from a time-dependent Lorentz invariant Dirac equation by xing the reference frame to the nuclear one and then factor o time by standard mathematical procedures. A fully relativistic approach to many-electron systems is on the other hand necessarily time-dependent due to the eects of retardation in the two-electron interaction, which in principle requires a complete specication of the history of the interacting particles.

Relativistic Hamiltonians for chemistry are therefore chosen according to their ability to capture essential relativistic eects combined with their ease of use.

The development of exact 2-component relativistic Hamiltonians undoubtedly constitute a major breakthrough of relativistic quantum chemistry due to the simplicity of their construction and the superior accuracy with respect to approximate 2-component relativistic Hamiltonians. Removing the need for a small component basis gives a speedup of an order of magnitude at the SCF and 4-index transformation level compared to 4-component calculations, although by suitable tricks it is possible to reduce the computational cost at the latter level to that of 2-component calculations, at least in a DFT framework [112]. At the correlated CI and CC level, the cost of 2-and 4-component calculations is strictly identical, since both levels employ the no-pair approximation. This is easily seen from the second quantized form of the electronic Hamiltonian (33): It is dened in terms of one-and two-electron integrals, so once the 4-index transformation has been carried out, it is not possible to deduce whether the orbitals came from a 2-or 4-component calculation.

What is lost with respect to 4-component theory is the real-space representation of quantities such as the charge and current density, since the X2C Hamiltonian and property operators are formulated in matrix algebra. One possible option is to backtransform 2-component orbitals to the 4-component level for constructing such objects. Another complication with respect to the 4-component level is that when calculating energy derivatives for molecular properties and gradients, such derivatives have to be taken with respect to the decoupling transformation as well [130,131]. Some consideration is also required for properties, such as NMR parameters, where the response contribution from negativeenergy orbitals can not be neglected [132134].

With respect to approximate 2-component Hamiltonians the X2C Hamiltonian further benets from the ease of picture change transformation, making the calculation of molecular properties much more simple to implement. One may ask what will be the future for approximate Hamiltonians such as ZORA and DKH2. Both are widely distributed in computer codes and will therefore continue to be used. I personally believe that further development of the DKH Hamiltonian will eventually cease, whereas the ZORA Hamiltonian may have a brighter future due to the more transparent decoupling transformation. All calculations were carried out by the dirac code [START_REF] Dirac | a relativistic ab initio electronic structure program, Release DIRAC10[END_REF] and employed an uncontracted 24s19p12d9f large component Gaussian basis [128,129] and, where relevant, a small component basis generated by restricted kinetic balance. For an explanation of acronyms, see text.

With the development of ecient 2-component methods for relativistic molecular electronic structure calculations, one may ask if relativity in chemistry has been solved, and so all that remains in this domain is the straightforward and possibly tedious transfer of the full repertoire of non-relativistic methodology to the relativistic domain. Personally I do not believe so. There are challenges ahead, notably exploring the interface with QED [START_REF] Kutzelnigg | [END_REF][START_REF] Saue | Theoretical Chemistry and Physics of Heavy and Superheavy Elements[END_REF], which may be of chemical relevance [? ], and the relativistic denition of certain molecular properties, such as the spin-rotation constant [135]. A major challenge is also to develop ecient electron correlation methods which incorporate spin-orbit interaction from the start [136]. Faced with bombastic statements about the completeness of relativistic molecular quantum mechanics it may be appropriate to end this mini-review with the following quote from a speech of Albert A.

Michelson at the dedication of Ryerson Physics Lab, University of Chicago, in 1894 : The more important fundamental laws and facts of the physical sciences have all been discovered and these are so rmly established that the possibility of ever being supplanted in consequence of new discoveries is exceedingly remote . . .Our future discoveries must be looked for in the sixth place of decimals. I believe that nobody knows what surprises relativistic molecular quantum mechanics may have in store for us.

Figure 1 .

 1 Figure 1. A 3D chart of theoretical model chemistries

  will give a numerical illustration, where all results have been obtained with the same computer code and the same basis set, to illustrate the performance of the various Hamiltonians discussed in the text. Finally, I will discuss what has been gained by the introduction of the X2C Hamiltonian, and also what are its possible disadvantages. I will conclude by a personal outlook on the future of 2-component relativistic Hamiltonians and of relativistic quantum chemistry in general.

Figure 2 .

 2 Figure 2. Basic mechanism of spin-orbit interaction generated by a clamped nucleus. The arrow represents a Lorentz transformation from the frame of clamped nuclei to the frame of the moving electron.

Figure 3 .

 3 Figure 3. Basic mechanism of spin-orbit interaction between electrons. The arrows represents Lorentz transformations between the frame of the moving electrons and the frame of clamped nuclei.

Most 4 -

 4 component relativistic molecular calculations are based on the Dirac-Coulomb Hamiltonian, where the Dirac Hamiltonian h D in the molecular eld (

Figure 4 .

 4 Figure 4. The small component 1s radial function of the mercury atom (nite Gaussian nucleus) generated from the large component using exact coupling or restricted kinetic balance. Radial expectation values of the 1s 1/2 and 2s 1/2 orbitals are indicated for reference. The data is based on a numerical 4-component relativistic Hartree-Fock calculation using the grasp code [51].

Figure 5 .

 5 Figure 5. A pictorial model of Zitterbewegung.

Figure 6 .

 6 Figure 6. The factor K(r) of the ZORA kinetic energy operator plotted as a function of radial distance (in a0) for the mercury atom with point charge nucleus.

component

  Hamiltonian. Applying the same projectors to the two-electron part in 4-component calculations therefore mimics 2-component calculations with the correctly transformed two-electron operator. The above observation also tells us that a atom-or fragment-based construction of the X2C Hamiltonian should preferably start from converged Fock (Kohn-Sham) matrix to allow maximally relaxed projectors. E. Picture change errors We have seen in section III A that 2-component Hamiltonians are obtained by taking the ++ block of the transformed 4-component parent Hamiltonian, that is

  Our reference will be the results obtained with the 4-component Dirac-Coulomb (DC) Hamiltonian, but I also display results obtained with the Dirac-Coulomb-Gaunt (DCG) Hamiltonian. When adding the Gaunt operator one notes an appreciable upward shift of the orbital energies of the inner core orbitals. Spin-orbit splittings ∆ SO are somewhat reduced due to the inclusion of spin-other orbit interaction. As discussed in section II C, SCF calculations based on the DC (or DCG) Hamiltonian excludes negative-energy orbitals from the construction of the mean-eld potential, which corresponds to the implicit use of projection operators, as suggested by Mittleman[66]. In TableII also report results obtained with other choices of projection operators. It can be seen that the use of projection operators P f ree constructed from solutions to the free-particle Dirac equation introduce signicant errors compared to the fully relaxed DC results. The errors associated with the use of projection operators P V ext constructed from solutions of the bare-nucleus Dirac Hamiltonian, which denes the Dirac equation for an electron in the molecular eld[START_REF] Dolg | [END_REF], are on the other hand quite acceptable.

Figure 7 .

 7 Figure 7. Comparison of 4-and 2-component electron number density, the latter calculated without picture change, of the mercury atom, as a function of radial distance. All quantities in atomic units.

Figure 8 .

 8 Figure 8. Comparison of 4-and 2-component electron number density, the latter calculated without picture change, of the mercury atom, as a function of radial distance. All quantities in atomic units.

  The 2-component equation of the correctly normalized function φ is accordingly

Table I

 I have been generated using an analytic implementation based on(63). The scaled ZORA Hamiltonian, denoted scZORA in TableI, outperforms the DKH2 Hamiltonian and approaches the X2C Hamiltonian in accuracy.

Table I .

 I -3062.411 -3074.239 -3218.372 -3074.248 -3069.951 -3214.349 -3061.460 -3380.366 -3067.397 Orbital energies (in E h ) from Hartree-Fock calculations on the mercury atom based on various relativistic Hamiltonians.

		DCG	DC	DC(P f ree ) DC(PV ext) X2C(Vext)	DKH1	DKH2	ZORA	scZORA
	1s 1/2
	2s 1/2	-548.905 -550.250	-566.498	-550.252	-549.615	-565.900 -548.677 -561.253 -549.292
	2p 1/2	-524.536 -526.846	-531.525	-526.850	-526.451	-531.271 -526.045 -537.952 -526.461
	2p 3/2	-453.605 -455.153	-455.092	-455.155	-454.824	-454.786 -454.829 -463.518 -454.746
	3s 1/2	-132.861 -133.110	-136.662	-133.110	-132.974	-136.534 -132.769 -133.861 -132.904
	3p 1/2	-122.195 -122.634	-123.707	-122.634	-122.557	-123.659 -122.463 -123.439 -122.564
	3p 3/2	-106.261 -106.541	-106.476	-106.541	-106.477	-106.416 -106.482 -107.145 -106.460
	3d 3/2	-89.242	-89.432	-89.228	-89.432	-89.418	-89.216	-89.430	-89.980	-89.418
	3d 5/2	-85.887	-86.016	-85.817	-86.016	-85.989	-85.788	-86.000	-86.505	-85.981
	4s 1/2	-30.593	-30.644	-31.527	-30.644	-30.613	-31.498	-30.562	-30.670	-30.596
	4p 1/2	-26.026	-26.119	-26.371	-26.120	-26.104	-26.363	-26.082	-26.168	-26.107
	4p 3/2	-22.131	-22.184	-22.153	-22.184	-22.172	-22.141	-22.174	-22.215	-22.169
	4d 3/2	-14.767	-14.792	-14.733	-14.792	-14.792	-14.733	-14.796	-14.821	-14.793
	4d 5/2	-14.036	-14.048	-13.990	-14.048	-14.045	-13.987	-14.049	-14.070	-14.044
	4f 5/2	-4.476	-4.469	-4.424	-4.469	-4.472	-4.428	-4.475	-4.480	-4.474
	4f 7/2	-4.321	-4.307	-4.264	-4.307	-4.310	-4.266	-4.313	-4.316	-4.310
	5s 1/2	-5.091	-5.099	-5.267	-5.099	-5.094	-5.262	-5.084	-5.093	-5.090
	5p 1/2	-3.520	-3.533	-3.568	-3.533	-3.532	-3.568	-3.528	-3.534	-3.532
	5p 3/2	-2.832	-2.838	-2.822	-2.838	-2.837	-2.821	-2.838	-2.838	-2.836
	5d 3/2	-0.646	-0.646	-0.630	-0.646	-0.647	-0.631	-0.648	-0.648	-0.647
	5d 5/2	-0.572	-0.571	-0.556	-0.571	-0.572	-0.557	-0.573	-0.572	-0.572
	6s 1/2	-0.326	-0.326	-0.341	-0.326	-0.326	-0.340	-0.325	-0.325	-0.325
	∆SO(2p) 70.931	71.693	76.433	71.695	71.628	76.485	71.216	74.434	71.715
	∆SO(3p) 15.934	16.093	17.230	16.093	16.080	17.242	15.980	16.294	16.104
	∆SO(4p) 3.895	3.935	4.218	3.935	3.932	4.221	3.908	3.953	3.939
	∆SO(5p) 0.688	0.695	0.746	0.695	0.695	0.747	0.691	0.697	0.696
	∆SO(3d) 3.354	3.416	3.412	3.416	3.429	3.428	3.429	3.476	3.438
	∆SO(4d) 0.730	0.744	0.743	0.744	0.747	0.747	0.747	0.751	0.749
	∆SO(5d) 0.074	0.075	0.074	0.075	0.075	0.074	0.075	0.076	0.076
	∆SO(4f)	0.155	0.161	0.160	0.161	0.162	0.162	0.163	0.164	0.164