Ayaki Sunaga

Maen Salman

Trond Saue

eective QED 4-component relativistic Hamiltonian with eective QED potentials for molecular calculations

We report the implementation of eective QED potentials for all-electron 4-component relativistic molecular calculations using the DIRAC code. The potentials are also available for 2-component calculations, being properly picture-change transformed. The latter point is important; we demonstrate through atomic calculations that picture-change errors are sizable. Specicially, we have implemented the Uehling potential [E.

I. INTRODUCTION

Relativistic quantum chemistry is the proper framework for the theoretical study of heavy elements.

15 For example, the yellow color of gold, 6,7 as well as the cell potential of the lead-acid battery 8 cannot be explained without relativistic eects. Even for light elements, the ne structure of spectra is essentially due to spin-orbit (SO) interaction (e.g. Refs. 911).

Improvements in both computational power and methodology nowadays allow highly accurate electronic structure calculations including both relativistic and electron correlation eects. A next challenge for increased accuracy is the inclusion of the eects of quantum electrodynamics (QED), which in principle means going beyond the no-pair approximation. 5,12,13 We focus on QED eects generating the Lamb shift, roughly described as follows:

vacuum polarization (VP): a charge in space is surrounded by virtual electron-positron pairs and this a) euthor to whom orrespondene should e ddressedX trondFsuedirsmFupsEtlseFfrFY httpXGGdirFupsEtlseFfrGsue contributes to its observed charge the electron self-energy (SE): the electron drags along its electromagnetic eld and this contributes to its observed mass.

For hydrogen the splitting between the 2 S 1/2 and 2 P 1/2 states is a mere 4 meV, 14 but for U 91+ it has grown to a whopping 468 eV.

15 QED eects would possibly con- stitute the nal correction to chemistry concerning the fundamental inter-particle interactions because the next contribution, parity non-conservation (PNC) associated with the weak force, is typically ten orders of magnitude smaller.

16 The magnitude of QED eects has been esti- mated based on the ionization potential of alkali atoms, and the rule of thumb is that QED eects reduce relativistic eects by about one percent. 17 The rigorous QED approach for few-electron systems cannot be extended to many-electron systems because of the high computational cost. A more practical, but approximate approach is the introduction of eective QED potentials (eQED). 2327 In the atomic case, some codes for the calculation with eective potentials have been reported (e.g., GRASP, 28 QEDMOD, 29,30 and AMBiT 31 ).

A nice illustration is the recent work by Pa²teka and coworkers, 32 which was nally able to bring the calculated ionization potential (IP) and the electron anity (EA)

of the gold atom into meV agreement with experiment, high-order electron correlation being the missing crucial ingredient.

For the case of molecules in chemistry, pioneering works have been done by Kirk Peterson's group. PPs are widely used for the inclusion of relativistic effects, and generally give accurate results for valence properties compared with all-electron calculations.

42,43 How- ever, the PP approach cannot be applied to molecular core-properties such as NMR and Mössbauer parameters, which bars the possibility to investigate the eect of QED in the nuclear region where such eects are generated. 44

The eQED approach promoted by the Peterson group can in principle be applied to core-properties, but it should be noted that eQED potentials were added to ap- 53 In addition to the Pb system, we have also calculated the heavier analogue, Fl hydrides.

Very recently, Leonid Skripnikov reported the implementation of eective QED potentials for 4-component all-electron molecular calculations, so far with a focus on transition energies.

54 The initial report has been followed by applications to Ba + , BaF, RaF and E120F 55 as well as the ve low-lying excited states of RaF.

56 The implemen- tation is to some extent complementary to ours in that it uses the eective SE potential proposed by Shabaev and co-workers.

26,27,29 Interestingly, the implementation is based on the DIRAC code as well.

The paper is organized as follows: in Sec. II we review the eective QED potentials that we have implemented, and in Sec. III we discuss the numerical integration of these potentials. This is followed by Sec. IV which gives the computational details of our calculations. Our results are presented in Sec. V, followed by conclusions in Sec. VI. We also provide an appendix with more extensive theory and reading suggestions. SI-units are used throughout this paper.

II. THEORY

The starting point for our work is an electronic Hamiltonian on the generic form

H = V NN + i H D (x i ) + 1 2 i̸ =j g(x i , x j ) (1)
where V NN is the classical repulsion of xed nuclei. The one-electron part is the Dirac Hamiltonian

H D (x i ) = (β -1 4 )m e c 2 -iℏcα • ∇ i -eφ N (x i ), (2)
in the electric potential φ N of the xed nuclei and shifted by -m e c 2 to align energies with the non-relativistic scale. In the present work, the two-electron interaction g will be the instantaneous Coulomb term supplemented with the Gaunt term. 57 Further discussion of the result- ing DiracCoulombGaunt (DCG) Hamiltonian is for instance found in Ref. 3.

Our goal is to introduce QED eects, notably electron self-energy (SE) and vacuum polarization (VP), by extending the one-electron Hamiltonian by the corresponding eective QED potentials

H D → H D -eφ effQED ; φ effQED = A φ SE A + φ VP A .
(3)

Note that the eective QED potentials are formulated as a sum over atomic contributions due to their expected short-range nature (on the order of a reduced Compton wavelength λ = ℏ/m e c). 44

In the following we shall present the eective QED potentials selected for our implementation with some remarks on their construction which may provide indications on their expected performance. We shall proceed within the Ŝ-matrix (scattering matrix) formalism of QED. Since we hope to address a wider audience than QED specialists, we provide a more extensive theoretical background in Appendix A.

QED is the relativistic quantum eld theory that describes the interaction of electromagnetic radiation with relativistic matter (Dirac electrons). The interaction between electrons and photons is given by an interaction Hamiltonian density

ĤI (x) = -ec Ψ (x) γ µ Ψ (x) µ (x) . (4)
Here, Ψ (x) and Ψ (x) are the quantized Dirac eld operator and its corresponding adjoint, whereas µ (x) is the quantized photon eld operator. The job of these operators is to create and annihilate, at the spacetime point x = (ct, x), electrons and photons, respectively.

This last expression accounts (explicitly) for the coupling between electron and photon elds, and is obtained through minimal substitution of the four-gradient of the ordered ones with all possible contractions, which in turn can be translated into the iconic Feynman diagrams. 63

We limit attention to systems of n electrons and zero photons (photon vacuum). The latter implies that any string of normal-ordered photon operators µ (x) that is not fully contracted will vanish upon taking expectation values, such that the Ŝ-matrix expansion is ef-fectively limited to even-ordered contributions, associated with the ne-structure constant α = e 2 /4πε 0 ℏc as expansion parameter. To lowest order in α appears ve Feynman diagrams, shown in Fig. 1 ), rather than free-particle ones. In the atomic case, this provides us with a second perturbation expansion parameter Zα, as will be seen in the next section.

A. Eective QED potentials for vacuum polarization

The four-potential associated with the vacuum polarization eect can be written as

φ µ VP (x 1 ) = e 2πi C F dz d 3 x 2 Tr γ µ G A e (x 2 , x 2 ; z) γ 0 4πϵ 0 |x 1 -x 2 | , (5)
where the complex z-integral is to be evaluated along the Feynman contour C F that goes above and below positive-and negative-energy poles, respectively, of the bound electron Green's function G A e . This function is related to the bound electron propagator S F A e by Eq. (A42).

This VP potential leads to the following vacuum polarization energy-shift

∆E α,2 VP = -e i d 3 x ψi (x) γ µ ψ i (x) φ µ VP (x) . (6)
From consideration of time-reversal symmetry, one can show that in the case of a purely scalar external potential: A e = (φ e /c, 0), vector components of the vacuum polarization four-potential vanish

φ µ VP (x) = 0 for µ = 1, 2, 3. (7)
The bound Green's function G A e can be written in terms of the free Green's function G 0 and expanded in powers of the time-independent external potential A e (hence @A ingleEphoton exhngeF @A uum polriztionF @A elfEenergyF @dA hiret ule digrmF @eA ixhnge ule digrmF FIG. 1. The lowest-order QED corrections for a many-boundelectron system: of order α.

Zα in the atomic case) as shown in Eq. (A44). As discussed in Sec. A 6 b, the rst non-vanishing term of this expansion is the one that is linear in the external potential A e (x) (the one-potential term). The potential of Eq. ( 5) is divergent (as seen in Sec

φ point Ueh. (x) = Ze 4πε 0 r x 2α 3π K 1 2r x λ , (8)
that corrects the classical Coulomb potential. Here, r x ≡ |x| is the radial distance and expressed in terms of the function 66 (10)

K 1 (x) = ∞ 1 dζe -xζ 1 ζ 2 + 1 2ζ 4 ζ 2 -1. ( 9
In the case of a spherically symmetric nuclear charge distribution, one obtains, after angular integration 66

φ nuc. Ueh. (x) = Ze 4πε 0 r x λ 2α 3 ∞ 0 r y dr y ρ nuc. (r y ) × K 0 2 λ |r x -r y | -K 0 2 λ |r x + r y | , (11)
where appears the function

K 0 (x) = ∞ 1 dζe -xζ 1 ζ 3 + 1 2ζ 5 ζ 2 -1. ( 12
)
The integral functions K 0 and K 1 are related through

K 1 (x) = - d dx K 0 (x).
(

) 13
The Uehling potential generally represents the dominant vacuum polarization eect.

19,69

The Feynman di- agram associated with this process is presented in Fig. 3b, and associated with the α (Zα) perturbation order.

The higher-order vacuum polarization potentials, associated with the WichmannKroll: 70 α (Zα)

B. Eective QED potentials for self-energy

The energy-shift associated with the self-energy process, in which an electron emits and absorbs a virtual photon, is given by the following expression

∆E α,2 SE = -e i d 3 x 1 d 3 x 2 × ψ † i (x 2 ) φ SE (x 2 , x 1 ; E i ) ψ i (x 1 ) (14)
φ SE (x 2 , x 1 ; E i ) = - e 2πi C F dzα µ G A e (x 2 , x 1 ; z) α µ × exp g (x 2 , x 1 ; z -E i ) 4πϵ 0 |x 1 -x 2 | , g (x 2 , x 1 ; z) = + i ℏ |x 1 -x 2 | z 2 /c 2 + iϵ. (15)
Here, ϵ is a small positive number, and the z-integral is again to be evaluated along the Feynman contour C F . This expression is obtained using the covariant Feynman gauge photon propagator.

The corresponding expression obtained using Coulomb gauge photon propagator is given by Lindgren in Ref. 73 (section 4.6.1.2) (See also Malenfant in Ref. 74). As in the vacuum polarization case, the self-energy potential of Eq. ( 15) is divergent (as seen in Sec. A 6 c), and calls for a regularization and renormalization treatment in order extract the physical (nite) correction; see Sec. A 6 d.

In the next two sections, we shall assume that the nonlocal potential of Eq. ( 15) can be written in terms of a local eective potential φ SE (x 1 ) as

φ SE (x 2 , x 1 ; E i ) ≈ φ SE (x 1 )δ(x 2 -x 1 ), (16)
and discuss some choices of φ SE (x 1 ) that are designed to reproduce some precise self-energy correction calcula- In the evaluation of matrix elements over the operator of Eq. (A68), Flambaum and Ginges employ free-particle q p 1 p 2 @A glssil sttering proessF q p 1 p 2 @A pirst rditive orretionF FIG. 2. Momentum-space Feynman diagrams for the lowest-order scattering processes through (the exchange of) momentum-transfer q = p2 -p1.

solutions rather than atomic bound orbitals. This replacement yields the free-electron vertex-correction (VC) problem. This terminology can be understood from consideration of the scattering of a free electron due to the interaction with a classical external potential (the vertex process). In terms of momentum-space quantities (cf. Eq.(A65)), including free electron eld operators, the corresponding (non-radiative) S-matrix is given by

Ŝ(1) scattering = - e iℏ d 4 p 2 (2πℏ) 4 d 4 p 1 (2πℏ) 4 × : Ψ (p 2 ) γ µ A e µ (p 2 -p 1 ) Ψ (p 1 ) : (20)
(see for instance section 8.7 in Ref. 78). This process is represented in the left panel of Fig. 2, where the wiggly line ending with a cross × describes an interaction of a free electron with the classical external potential source through the exchange of a four-momentum q = p 2 -p 1 .

We note that in general a factor (-eγ µ ) is associated with each spacetime point (vertex). The right panel of Fig. 2 represents one of the four lowest-order radiative corrections to the left panel process. The corresponding S-matrix can be combined with the one of Eq. ( 20)

through the substitution

γ µ → Γ µ = γ µ + Λ µ (p 2 , p 1 ), (21)
where the vertex-correction function Λ µ (p 2 , p 1 ) is given by Eq.(A69). After a careful treatment of the divergence when q = 0, as done in Refs. 60 (section 6.3) and 79 (section 117), one obtains the regularized (physical) vertex-correction function Λ µ R (p 2 , p 1 ). Furthermore, using the fact that the vertex function is sandwiched between free-electron (on-mass-shell) eld operators, one can show that this function can be written as

Λ µ R (q) = γ µ F 1 q 2 + i 2m e c σ µν q ν F 2 q 2 , ( 22
)
where Since the free-electron vertex function of Eq. ( 22) only depends on the momentum-transfer q = p 2 -p 1 , it conveniently yields a local potential in real space. This can be clearly seen from the following relation

σ µν = i 2 [γ µ , γ ν ],
d 4 p 2 (2πℏ) 4 d 4 p 1 (2πℏ) 4 Ψ (p 2 ) Λ µ R (q) A e µ (q) Ψ (p 1 ) = 1 c d 4 x Ψ † (x) φ VC (x) Ψ (x) . ( 23
)
When the nucleus is described as a point charge the corresponding Coulomb potential,

A e 0 (q) = δ(q 0 ) 2πℏ 3 cϵ 0 Ze q 2 , ( 24
)
generates a vertex-correction potential of the form

φ point VC (x) = ℏ 2 ϵ 0 d 3 q (2πℏ) 3 e + i ℏ q•x Ze q 2 × F 1 -q 2 + 1 2m e c γ • qF 2 -q 2 = φ point elec (x) + φ point mag (x) , (25)
which splits into electric and magnetic scalar potentials. We note that due to the time-independence of the Coulomb potential, the time-like part of the 4momentum transfer q = p 2 -p 1 vanishes. In terms of the variable t = q 2 = -q 2 , the form factors are Hermitian analytic functions, 81 that is

F (t) = F * (t * ). (26)
This feature, combined with these functions being radial in terms of q and the clever use of complex analysis techniques, allowed Berestetskii et al. to express such functions in coordinate-space using only their imaginary parts in momentum-space

F (x) = 1 (2πℏ) 2 r x ∞ 4m 2 e c 2 dt ℑ [F (t)] exp - 1 ℏ r x √ t . (27)
(see eq.(114.4) of Ref. 79). It may be noted that the lower limit of integration over t is 4m 2 e c 2 , corresponding to the threshold of pair creation.

82 Expressions for the imagi- nary part of the form factors can be found in Refs. 79 (eq.(117.14-15)) and 83 (eq.(2.12))

ℑ [F 1 (t)] = α t (t -4m 2 e c 2 )
2m 2 e c 2 -3t/4

(28)

+ t/2 -m 2 e c 2 log t -4m 2 e c 2 λ 2 , ℑ [F 2 (t)] = αm 2 c 2 t (t -4m 2 e c 2 ) . ( 29
)
Building on the work of Berestetskii et al.,

79 Flambaum and Ginges have evaluated the integral of Eq. ( 25), and obtained the associated real-space potentials. After the variable substitution t = 4m 2 e c 2 ζ 2 , the magnetic potential was found to be

φ point mag (x) = αℏ 4πm e c iγ • ∇ x Ze 4πε 0 r x K m 2r x λ -1 , (30)
where we have introduced the function

K m (x) = ∞ 1 dζ e -xζ ζ 2 ζ 2 -1 , (31)
which can be recognized as the 2nd BickleyNaylor function Ki 2 (cf. Ref. 68). Note that the same variable ζ is employed in the Uehling potential (cf. Eqs. ( 9) and ( 12)).

The magnetic contribution gives the rst-order correction to the magnetic moment: the anomalous magnetic moment of the electron, rst calculated by Schwinger, see for instance Mandl in Ref. 78 (section 10.5).

On the other hand, the electric form factor yields the electric eective potential

φ point elec (x, λ) = - α π Ze 4πε 0 r x K e 2r x λ , (32)
where we have introduced the function

K e (x) = ∞ 1 dζ e -xζ ζ 2 -1 - 3 2 + 1 ζ 2 (33) + 1 - 1 2ζ 2 ln ζ 2 -1 + 2 ln 2m e c 2 λ .
These self-energy eective potentials where rst derived with respect to a point nucleus (Coulomb potential), and the corresponding generalized expressions for an arbitrary normalized nuclear distribution ρ nuc. are obtained by convolution, 84 as in Eq. ( 10).

The potential of Eq. ( 32) is called the high-frequency term, because it contains an energy parameter λ, already present in Eq. ( 28), that prevents the obtention of low-frequency contributions. This parameter is associated with the introduction of a small ctitious photon mass, which needs to be plugged in the photon propagator denominator in order to make the divergent (at small momenta) momentum-space integral, associated

φ point HF (x) = A (Z, x) φ point elec (x, λ) , (34)
where A (Z, x) is a tting function and choose a λ-value that will minimize the low-frequency contribution. They argue that λ should be on the order of electron binding energies, that is (Zα) 2 m e c 2 . They nally dene it through

ln 2m e c 2 λ = 2 ln 1 Zα + 1 2 , (35)
though, for better performance. Flambaum and Ginges next argue that the low-frequency (LF) potential should have the range of a 1s orbital of hydrogen-like atoms and therefore choose the functional form

φ point LF (x) = - B (Z) e Z 4 α 5 m e c 2 e -Zrx/a0 , ( 36
)
where a 0 = λ/α is the Bohr radius and

B (Z) = 0.074 + 0.35 × Zα, (37)
is a second tting function.

The tting function of the high-frequency contribution is written as 1), an operator on the form 30) and ( 32) by the potentials of nite nuclear charge distributions, and we have so far followed this approach which appears to be a reasonable approximation, as can be inferred from Table IV of Ref. 84.

A (Z, x) = Θ (Z, x) 1.071 -1.976y 2 -2.128y 3 + 0.169y 4 , ( 38
V ASHIFT = i |ψ i ⟩ω i ⟨ψ i |, ω i = ⟨ψ i | -eφ effQED |ψ i ⟩, (40)
We have adapted the GRASP eective QED potential routines to molecular calculations by using the numerical integration scheme implemented for relativistic KohnSham calculations in the DIRAC molecular code.

49

The scheme is based on the Becke partitioning 95 of the molecular volume into atomic ones for which numerical integration is carried out in spherical coordinates. Specifically, we use Lebedev angular quadrature, 96 by default setting ℓ = 15, combined with the basis-set adaptive radial grid proposed by Lindh and co-workers. 97 It may be noted that the eective QED potentials presented in the previous section are all radial, with the exception of the magnetic contribution to the FlambaumGinges SE potential, Eq. ( 30).

Due to the very local nature of the eective QED potentials 44 one-electron integrals over a potential asso- ciated with atomic center A can be well approximated by

V A µν ≈ R A 0 dr A Ω dΩ A χ µ v A χ ν (r A )r 2 A , (41)
where {χ µ } are Gaussian-type basis functions. The most delocal potential is the low-frequency contribution to the electric form factor of the FlambaumGinges SE potential, Eq. ( 36), since it has been designed to have the range of the 1s orbital of a hydrogen-like atom. For low Z the potential may thereby overlap signicantly with neighbor centers. By default, we therefore deactivate the eective QED potentials for Z < 19. We also determine the value of the upper limit of radial integration R A based on the convergence of the low-frequency term to a very conservative 10 -50 . For the atomic calculations reported in 114 We used the counterpoise correction 115 to minimize basis set superposition errors (BSSE).

IV. COMPUTATIONAL DETAILS

For the calculations of gold cyanide, we used the CCSD(T) method for comparison with experiment. In the CCSD(T) calculation, 4f 5s5p5d6s V. RESULTS

A. Atomic calculations

In Table I Au C N ω ⟨ε⟩ 5d 5/2 6s 1/2 2s 1/2 2p 1/2 2p 3/2 2s 1/2 2p 1/2 2p 3/2 3/2 -0.342 -0.01 0.00 0.00 0.00 0.92 0.00 0.00 1.09 1/2 -0.345 0.00 0.00 0.00 0.65 0.28 0.00 0.64 0.43 1/2 -0.586 0.15 0.39 0.92 0.17 0.34 0.01 -0.01 -0.02 1/2 -0.781 0.00 0.01 0.24 0.12 0.38 0.15 0.44 0.65

B = ℏ 4πI ⊥ ; I ⊥ = I xx = I yy = A m A z 2 A , (42)
when the molecular axis is aligned with the z-axis. z A is the distance of atom A from the center of mass. Eective r 0 and substitution r s structures are both obtained by assuming identical structures for all isotopomers of

|z A | = ℏ 4πµ 1 B A ′ - 1 B A ; 1 µ = 1 M + 1 ∆m A . ( 43
)
where M is the total mass of the parent isotopomer.

In

B ξ ν = B ξ e - i α ξ i ν i + d i 2 + 1 2 i,j γ ξ i,j ν i + d i 2 ν j + d j 2 + • • • . (44)
Here, ξ is the axis of rotation, α ξ and γ ξ are vibrationrotation interaction constants of dierent orders and d i is the degeneracy of vibration mode i.

C. van der Waals dimers

As a second molecular application of our implementation we consider spectroscopic constants of dimers with van der Waals bonding (M 2 , M = Hg, Rn, Cn, Og). In XV andXVI, the ratio of VP and SE is ∼ 1:-4.2 for the Pb system, while it is ∼ 1:-1.8 for the Fl system. Discussion along these lines is also found in Refs. 55,89. Finally, we note from Tables XV andXVI that the the atomic shift operator (ASHIFT), either using atomic ground state occupations or the eective atomic conguration in the molecules given in Table XIII, is not reliable for describing QED eects in the molecules.

VI. CONCLUSIONS

We have implemented eective QED potentials for relativistic molecular calculations by grafting code from the numerical atomic code GRASP onto the DFT grid of DIRAC. A general disadvantage of numerical integration is higher computational cost than analytical evaluation, to the extent that such expressions are available, although the implementation itself is easier and considerable savings are achieved by the locality of the eective QED potentials.

We report several applications of the new code, mostly using the the molecular mean-eld approximation Hamiltonian (X2Cmmf ). We demonstrate (Table II For the rare-gas dimers Hg 2 and Rn 2 we nd that QED increases bond lengths by about 0.15 pm. For the superheavy homologues the bond length increase is on the order of 0.30 pm; the eect on dissociation energies is quite small (∼0.4 %).

We have also investigated the eect of QED on the reaction energy of XH 4 --→ XH 2 + H 2 , (X=Pb,Fl). From projection analysis we do nd that there is a signicant change of valence s population of the metals during the reaction, in line with the proposition of Dyall and coworkers.

53 Interestingly, though, we also nd that in the tetrahydrides the valence s population essentially resides in bonding orbitals, but in non-bonding ones in the dihydrides. We nd for the dissociation of lead tetrahydride that QED reduces the magnitude of the reaction energy by 0.32 kcal/mol (-1.27 %); for the superheavy homologue the magnitude of the QED eect is basically the same (0.35 kcal/mol). This possibly surprising observation is explained by the reduction of the (negative) SE/VP ratio with increasing nuclear charge.

For these metal hydrides, and also AuCN, we have also tried a simpler approach for the incorporation of QED effects in molecular calculations in the form of an atomic shift operator, but we nd that this is not a reliable approach.

We would like to stress that our implementation of effective QED potential is general in the sense that they are available in all parts of the code. A natural continuation of our project will therefore be to explore the impact of these potentials on molecular properties probing electron density in the vicinity of nuclei, where the QED eects are generated. Our results so far indicate that QED eects may be more important than for the valence properties reported in the present work. For instance, the QED eect on the parity violation energy of

γ µ γ ν + γ ν γ µ = 2g µν 1 4 . (A1)
In Dirac basis they are represented by γ 0 = β and γ = βα. Following Lindgren 73 we shall complement the Dirac α matrices with α 0 = 1 4 to form a pseudo-4-vector.

We nally note that we put hats (ˆ) on quantities that contain creation/annihilation operators acting on occupation number states. Contrary to conventional QED sources, we have decided to express the formalism in full SI units.

Electron eld operator

The electron eld operator is given by the following annihilation expansion over all solutions of the Dirac equation

Ψ (x) = i ψ i (x) c i , with ψ i (x) = ψ i (x) e -i ℏ Eit . (A2)
In this expression, c i is the electron annihilation operator obeying the fermionic algebra relations

{c i , c † j } = δ ij , and {c i , c j } = {c † i , c † j } = 0, (A3)
and associated with the i-th spatial wavefunction ψ i (x) and energy-level E i that solve the time-independent Dirac equation in the presence of a time-independent external four-potential A e = (φ e /c, A e )

H D ψ i (x) = E i ψ i (x) ; with H D = cα • (-iℏ∇ + eA e (x)) -eφ e (x) + βm e c 2 . (A4)
The electron vacuum state is dened to be the one that vanishes after any annihilation: c i 0 e = 0, ∀i.

c i = a i , for E i > 0, c i = b † i , for E i < 0. (A6)
Here, operators a i and b i are introduced to distinguish between the particle (electron) and its hole (positron), (A17)

Ψ (x) = Ei>0 ψ i (x) a i + Ei<0 ψ i (x) b † i . ( A7
Finally, we note that the photon vacuum state is dened to be the state that satises the following relation

a (k, r) |0 p ⟩ = 0, ∀ k, r → Â+ µ (x) |0 p ⟩ = 0, ∀ µ, x. (A18)
We shall now consider the interaction between the noninteracting electron and photon elds and show how one can derive QED corrections using perturbation theory.

Perturbation theory

As in conventional perturbation theory, one wants to get the eigensolutions of the following total Hamiltonian

ĤS = Ĥ0

S + λ Ĥ1 S .

(A19)

The zeroth-order Hamiltonian

Ĥ0 S = Ĥ0 electron + Ĥ0 photon , (A20)
represents the free electron and photon elds. The electronic part is given by a spatial integral over the normalordered Dirac Hamiltonian density

Ĥ0 electron = d 3 x : Ψ † (x) H D (x) Ψ (x) : = Ei>0 E i a † i a i - Ei<0 E i b † i b i , (A21)
where in BSQED, H D is the Dirac Hamiltonian in the presence of the external four-potential A e = (φ e /c, A e ),

given in Eq. (A4), and where normal-ordering is indicated by double dots. The free photon Hamiltonian is written as an integral of the electromagnetic Hamiltonian density The perturbation Hamiltonian Ĥ1

Ĥ0 photon = 1 µ 0 d 3 x : -(∂ 0 µ (x))(∂ 0 µ (x)) + 1 2 (∂ ν µ (x))(∂ ν µ (x)) : = k 3 r=0 ℏω k ζ r a † (k, r) a (k, r) .
S complicates the problem, and prevents us from obtaining eigensolutions of the full Hamiltonian ĤS . λ is a dimensionless parameter that can be varied between 0 and 1, and which keeps track of the perturbation-order. This parameter is to be taken to 1 in order to account for the full perturbation by the end of the calculation. Notice that so far our

Hamiltonians have an S subscript; this is made to indicate that they are in the Schrödinger picture of quantum mechanics. Assuming that we know the eigensolutions of the unperturbed time-independent problem equation

Ĥ0 S |Φ α 0 ⟩ S = E α 0 |Φ α 0 ⟩ S ; |Φ α 0 (t)⟩ S = e -iE α 0 t/ℏ |Φ α 0 ⟩ S , (A23)
where the α superscript labels solutions (states and associated energy-levels), the ultimate goal is to nd eigensolutions of the perturbed problem

ĤS |Φ α ⟩ S = E α |Φ α ⟩ S ; |Φ α (t)⟩ S = e -iE α t/ℏ |Φ α ⟩ S . (A24)
GellMann and Low provided a closed form of the perturbed eigensolutions (E α , |Φ α ⟩ S ) in terms of the unperturbed ones (E α 0 , |Φ α 0 ⟩ S ) and the time-evolution operator;

154 see also Refs. 155 (pages 61-64) and 59 (sec- tion 11f.). A few years later, Sucher 156 provided an ex- pression of the perturbation energy-shift that is more symmetric in time

∆E α = E α -E α 0 = lim ϵ→0 λ→1 iϵλ 2 ∂ ∂λ log Φ α 0 Ŝ (ϵ, λ) Φ α 0 , (A25)
where ϵ is an energy-parameter, to be shortly discussed. This energy-shift expression contains the S-matrix operator that is dened to be the time-evolution operator that takes the interaction state from the very past t = -∞ to the very future t = +∞, and can be written as (see Dyson in Ref. 157 eq.( 4))

Ŝ (ϵ, λ) = T exp λ iℏc d 4 xe -ϵ ℏ |t| ĤI (x) . (A26)
In this expression, T stands for time-ordering, i.e., it re-orders the inside operators such that those associated with earlier times act rst. In the simplest case of two operators, the time-ordering operation is dened to be

T Â (x 1 ) B (x 2 ) ≡ Θ (t 1 -t 2 ) Â (x 1 ) B (x 2 ) ± Θ (t 2 -t 1 ) B (x 2 ) Â (x 1 ) , (A27)
where the minus sign applies when both operators  and B are of fermionic nature. Furthermore, the S-matrix in Eq. ( A26) is a functional of the interaction-Hamiltonian density ĤI (x), that is related to the interaction Hamiltonian Ĥ1 I (t) by the following integral

Ĥ1 I (t) = d 3 x ĤI (x) . (A28)
hold.

These relations show that the only non-zero contractions are between electron eld operators and their adjoints. The free Feynman electron propagator S F 0 (x, y), corresponding to the case A µ (x) = 0 µ , can be written as

S F 0 (y, x) = lim ϵ→0 d 4 p (2πℏ) 4 e -i ℏ p•(y-x) S F 0 (p) , with S F 0 (p) = γ µ p µ + m e c p 2 -m e 2 c 2 + iϵ , (A37)
where S F 0 (p) is the Fourier transformed free-electron propagator. The role of the small positive number ϵ is to shift energy-poles (at the energy-momentum relation) with respect to the Feynman prescription.

Similarly, the contraction of two photon operators (of Eq. ( A9)) is dened by the following expression

µ (x 1 ) Âν (x 2 ) ≡ 0 p T µ (x 1 ) Âν (x 2 ) 0 p = iℏD F µν (x 1 , x 2 ) , (A38)
where D F µν (x, y) is the photon propagator in the Feynman gauge, is given by the following expression

D F µν (x, y) = lim ϵ→0 d 4 p (2πℏ) 4 e -i ℏ p•(x-y) D F µν (p) , with D F µν (p) = g µν D F (p) = - ℏ 2 cϵ 0 g µν p 2 + iϵ (A39)
and satises the Maxwell Green's-type equation

∂ σ ∂ σ D F νθ (x, y) = g νθ ϵ 0 c δ (x -y) . ( A40
)
This equation is obtained after imposing the Lorenz gauge condition, otherwise this propagator will not be invertible; see Schwartz in Ref. 163 (section 8.5). We should nally note that the F superscript on both propagators is added to indicate that these are Feynman propagators. This means that when writing the propagators as Fourier transforms, the energy-integrals are to be taken along the Feynman contour. Dierent choices of paths (contours) lead to dierent propagators (retarded and advanced), but they all satisfy the corresponding Dirac and Maxwell equations.

Bound electron propagator expansion

The bound Feynman propagator S F A e (x 2 , x 1 ) of Eq. (A34) can be expanded in powers of the external potential as (Refs. 86 eq.(2-119) and 164 eq.( 16)):

S F A e (x 2 , x 1 ) = S F 0 (x 2 , x 1 ) -d 4 x 3 S F 0 (x 2 , x 3 ) eA e µ (x 3 ) γ µ S F 0 (x 3 , x 1 ) + . . . , (A41)
and written in terms of the free-electron propagator S F 0 (x 2 , x 1 ) of Eq. (A37). It is worth noting that the bound Feynman propagator can be related to the bound Dirac Green's function G A e by the relation of Ref. 162 (eq.( 32)):

S F A e (x 2 , x 1 ) = 1 iℏ 1 2πi C F dz G A e (x 2 , x 1 ; z) γ 0 e -i ℏ z(t2-t1) , (A42)
This Green's function satises

[H D (x 2 ) -z] G A e (x 2 , x 1 ; z) = 1 4 δ (x 2 -x 1 ) ; (A43) cf.
Eqs. (A4) and (A35). Using Eqs. (A41) and (A42), and integrating over time variables, one obtains the potential expansion associated with the Green's function

G A e (x 2 , x 1 ; z) = G 0 (x 2 , x 1 ; z) + ec d 3 x 3 G 0 (x 2 , x 3 ; z) A e µ (x 3 ) α µ G 0 (x 3 , x 1 ; z) + . . . (A44)
where the free Dirac Green's function is given by: G 0 = lim A e →0 G A e . These two expansions are known as the potential expansion, where consecutive terms are known as the zero-one-and many-potential terms. The main utility of this expansion is that it allows the isolation of ultraviolet divergent integrals encountered when evaluating loop integrals, as done by Baranger et al., 72 and later by many authors working within BSQED theory.

No-photon BSQED energy-shifts

Using the S-matrix expansion of Eq. (A29), one can expand Sucher's energy-shift expression of Eq. (A25) in powers of the interaction-Hamiltonian density and write, following Mohr in Ref. 165 (eqs.( 18) and ( 31)),

∆E α = lim ϵ→0 λ→1 iϵλ 2 Φ α 0 Ŝ(1) (ϵ, λ) Φ α 0 + 2 Φ α 0 Ŝ(2) (ϵ, λ) Φ α 0 -Φ α 0 Ŝ(1) (ϵ, λ) Φ α 0 2 + O λ 3 , (A45)
where Ŝ(n) is given in Eq. (A29). We shall now con- sider a system of n electrons and zero photons (photon vacuum), represented by the following electron-photon state, labeled by α:

|Φ α 0 ⟩ = n α e , 0 α p . (A46)
We remind the reader that the electron eld operators entering in our expressions describe non-interacting electrons, in the presence of an external potential, as also seen in the zeroth-order electron Hamiltonian of Eq. (A21). As already pointed out in Section II, the electron-electron interaction arises from terms describing exchange of virtual photons between bound electrons.

Since the photon state is chosen to be the vacuum one, this means that any string of photon operators that is not fully contracted, will vanish under the photon vacuum expectation value of Eq. (A25). Following this reasoning one concludes that the rst non-vanishing QED correction comes from the second-order Ŝ(2) -matrix

Ŝ(2) (ϵ, λ) = - λ 2 2ℏ 2 c 2 d 4 x 1 d 4 x 2 × e -ϵ ℏ (|t1|+|t2|) T ĤI (x 1 ) ĤI (x 2 ) . (A47)
Using Wick's theorem, we expand the electron and photon time-ordered products, and replace operator contractions by corresponding propagators, following the contraction denitions of Eqs. (A34) and (A38). Fur-thermore, using the symmetry properties of the photon propagator of Eq. (A39):

D F µν (x, y) = D F µν (y, x) = D F νµ (x, y) , (A48)
the second-order S-matrix of Eq. (A47) can be shown to reduce to the following expression

Ŝ(2) (ϵ, λ) = - λ 2 e 2 2ℏ 2 d 4 x 1 d 4 x 2 × e -ϵ ℏ (|t1|+|t2|) F (x 1 , x 2 ) , (A49)
where the operator F (x 1 , x 2 ) contains the following ve QED corrections to the non-interacting problem:

F (x 1 , x 2 ) = iℏ D F µ1µ2 (x 1 , x 2 ) : Ψ(x 1 )γ µ1 Ψ(x 1 ) Ψ(x 2 )γ µ2 Ψ(x 2 ) : SP + 2ℏ 2 D F µ1µ2 (x 1 , x 2 ) Tr S F A e (x 2 , x 2 )γ µ2 : Ψ(x 1 )γ µ1 Ψ(x 1 ) : VP -2ℏ 2 D F µ1µ2 (x 1 , x 2 ) : Ψ(x 1 )γ µ1 S F A e (x 1 , x 2 )γ µ2 Ψ(x 2 ) : SE -iℏ 3 D F µ1µ2 (x 1 , x 2 ) Tr S F A e (x 1 , x 1 ) γ µ1 Tr S F A e (x 2 , x 2 ) γ µ2 D1 + iℏ 3 D F µ1µ2 (x 1 , x 2 ) Tr S F A e (x 2 , x 1 ) γ µ1 S F A e (x 1 , x 2 ) γ µ2 D2 (A50)
Finally, using Sucher's energy expression of Eq. (A45), the second-order energy-shift becomes

∆E α,2 = lim ϵ→0 λ→1 iϵλ Φ α 0 Ŝ(2) (ϵ, λ) Φ α 0 = ∆E α,2 SP + ∆E α,2 VP + ∆E α,2 SE + ∆E α,2 D1 + ∆E α,2 D2 . (A51)
Each of these terms will be discussed in the next sec- and represented in Fig. 1a, describes electron-electron interaction in its lowest-order, where an electron feels the existence of the other electron through the exchange of a single virtual-photon. After integrating over times t 1 and t 2 in Eq. (A49), and taking limits in Eq. (A51), this correction yields an instantaneous direct interaction-term, in addition to a retarded exchange interaction-term, analogous to the direct and exchange terms in the Hartree-Fock theory: 205)):

∆E α,2 SP = e 2 2 i,j d 3 x 1 d 3 x 2 ψi (x 1 ) γ µ ψ i (x 1 ) 1 4πϵ 0 |x 1 -x 2 | ψj (x 2 ) γ µ ψ j (x 2 ) Direct - e 2 2 i,j d 3 x 1 d 3 x 2 ψi (x 1 ) γ µ ψ j (x 1 ) e + i cℏ |Ei-Ej ||x1-x2| 4πϵ 0 |x 1 -x 2 | ψj (x 2 ) γ µ ψ i (x 2 ) Exchange (
∆E α VP = -e i d 3 x 1 ψi (x 1 ) γ µ ψ i (x 1 ) φ µ VP (x 1 ) . (A53)
We note that the vacuum polarization eect is local, i.e., it can be written as an expectation value of a local vacuum polarization four-potential

φ µ VP (x 1 ) = ieℏ d 3 x 2 Tr γ µ S F A e (x 2 , x 2 ) 4πϵ 0 |x 1 -x 2 | . (A54)
The energy expression of Eq. (A53) (as well as the last potential) is divergent due to the fact that This naive estimation usually overestimates the eective (divergence), that we shall call E, and this can be seen after further analysis of the integral in question. As a consequence, some supercially divergent integral can be eectively less divergent, or hopefully convergent. In the next Table XVII, we list the supercial and eective divergences of the rst vacuum polarization terms. The reader should notice that with higher-order terms, more propagators are included in the momentum-space integral, meaning that more denominator powers are added, and as a consequence, the integral becomes less divergent.

lim x1→x2 S F A e (x 2 , x 1 ) = ∞, (
We shall now focus on the rst non-vanishing vacuum polarization contribution ∆E α VP,1 . As seen in VP,3 and presented in Fig. 3d, is known as the WichmannKroll eect. 70 As seen in Ta- 9)). We nally note that the latter authors provided a good approximation of the corresponding potential, in order to render the numerical evaluation more practical.

@A α (Zα) 0 F @A α (Zα) 1 F @A α (Zα) 2 F @dA α (Zα) 3 F FIG.
3. First four bound-state vacuum polarization processes, obtained after expanding the bound propagator in powers of the external potential. A wiggly line ending with a cross × indicates an interaction with the external eld.

∆E α,2 SE = -e i d 3 x 1 d 3 x 2 × ψ † i (x 2 ) φ SE (x 2 , x 1 ; E i ) ψ i (x 1 ) (A60) @A α (Zα) 0 F @A α (Zα) 1 F @A α (Zα) 2 F @dA α (Zα) 3 F FIG.
4. First four bound-state self-energy processes, obtained after expanding the bound propagator in powers of the external potential.

where the self-energy potential is given by

φ SE (x 2 , x 1 ; E i ) = - e 2πi C F dzα µ G A e (x 2 , x 1 ; z) α µ × exp + i ℏ |x 1 -x 2 | (z -E i ) 2 /c 2 + iϵ 4πϵ 0 |x 1 -x 2 | . (A61)
Similar forms of this equation are provided by Schweber in Ref. 59 (eq.( 205)) and Mohr in Ref. 175 (eq.(2.6)).

Notice at this point that unlike the vacuum polarization case, the self-energy is a non-local eect, as seen from Eq. (A60), and this is the reason behind the complexity of its analytical and numerical evaluation. As in the vacuum polarization case, the self-energy potential is a divergent quantity and needs to be regularized. In order to isolate divergent terms, one can use the Green's function (propagator) expansion of Eq. (A44) and write the total energy-shift as

∆E α,2 SE = ∆E α,2 SE,0 + ∆E α,2 SE,1 + ∆E α,2 SE,2 + . . . , (A62)
where ∆E α,2 SE,i represents the process in which the inter- nal electron interacts i times with the external potential, it is thus associated with the α (Zα) SE,1 are known as the (free) self-energy and the vertex-correction processes. These two contributions are @PHPPAD IHFIHTQGSFHHUWTIVF 57 tF eF quntD roeedings of the oyl oiety of vondonF eries e A122D SIQ @IWPWAF 58 wF qellEwnnD sl xuovo gimento @IWSSEIWTSA 4D VRV @IWSTAF 59 F F hweerD An Introduction to Relativistic Quantum Field Theory @owD eterson nd gompnyD ivnstonD sllinoisD IWTIAF 60 wF eskin nd hF hroederD An Introduction to Quantum Field TheoryD prontiers in hysis @estview ressD IWWSAF 61 F qreiner nd tF einhrdtD Field Quantization @pringerE erlg ferlin reidelergD IWWTAF 62 qF gF ikD hysF evF 80D PTV @IWSHAF 63 hF uiserD emerin ientist 93D IST @PHHSAF 64 F eF wF hirD roeedings of the oyl oiety of vondonF eries e 136D RSQ @IWQPAF 65 F rF purryD hysF evF 81D IIS @IWSIAF 66 vF yne pullerton nd qF eF inkerD hysF evF e 13D IPVQ @IWUTAF 67 F ulrsfeldD hysis vetters f 66D VT @IWUUAF 68 eF wF prolov nd hF wF rdlwD iuropen hysil tournl f 85D IR @PHIPAF 69 F feierD F tF wohrD rF erssonD nd qF o'D hysF evF e 58D WSR @IWWVAF 70 iF rF ihmnn nd xF wF urollD hysF evF 101D VRQ @IWSTAF 71 qF uällén nd eF ryD elskF wtF pysF wedd 29 @IWSSAF 72 wF frngerD rF eF fetheD nd F F peynmnD hysF evF 92D RVP @IWSQAF 73 sF vindgrenD Relativistic many-body theory: a new eldtheoretical approachD Pnd edFD olF TQ @pringerD PHITAF 74 tF wlenfntD hysF evF h 35D ISPS @IWVUAF 75 F sndelito nd F tF wohrD hysF evF e 58D ITS @IWWVAF 76 he -δ in le I of efF PR orresponds to 100×SE Li /∆E Hf s Li with vlues tken from le P of efF UUD where SE Li in turn omes from efF IVSF 77 F fourd nd F sndelitoD he iuropen hysil tournl h 8D SW @PHHHAF 78 pF wndl nd qF hwD Quantum eld theory @tohn iley 8 onsD PHIHAF 79 F fF ferestetskiiD iF wF vifshitzD nd vF F itevskiiD Quantum Electrodynamics (Course of Theoretical Physics,4) @ergE mon ressD yxfordD IWVPAF 80 F rofstdterD evF wodF hysF 28D PIR @IWSTAF 81 gF fF vng nd xF ukerD Mathematische Methoden in der Physik @pektrum ekdemisher erlgD PHITAF 82 F tF idenD roeedings of the oyl oiety of vondonF eries eF wthemtil nd hysil ienes 210D QVV @IWSPAF 83 F frieriD tF eF wignoD nd iF emiddiD sl xuovo gimento e @IWTSEIWUHA 11D VPR @IWUPAF 84 tF qinges nd tF ferengutD hysil eview e 93D HSPSHW @PHITAF 85 F qreiner nd tF einhrdtD Quantum ElectrodynamicsD Rth edF @pringerEerlg ferlin reidelergD PHHWAF 86 gF stzykson nd tF uerD Quantum Field TheoryD hover fooks on hysis @hover ulitionsD IWVHAF 87 F tF wohrD hysF evF e 46D RRPI @IWWPAF 88 F tF wohr nd FEuF uimD hysF evF e 45D PUPU @IWWPAF 89 gF hierfelder nd F hwerdtfegerD hysil eview e 82D HTPSHQ @PHIHAF 90 yF possgrdD yF qropenD iF ilivD nd F ueD tF gompF hysF 119D WQSS @PHHQAF 91 F huillrdD tFEfF otD F ueD nd uF pegriD tF ghemF hysF 124D ISRQHU @PHHTAF 92 uF qF hyllD tF ghemF hysF 139D PIIHQ @PHIQAF 93 fF wuenzieD sF qrntD nd F xorringtonD 101 tF ikkemD vF issherD F ueD nd wF slisD he tournl of ghemil hysis 131D IPRIIT @PHHWAF 102 tF ipek nd F qF wezeyD tF ghemF hysF 90D RWIT @IWVWAF 103 qF uniziD tF ghemF heory gomputF 9D RVQR @PHIQAF 104 uF qF hyllD tF hysF ghemF e 113D IPTQV @PHHWAF 105 uF qF hyllD heorF ghemF eF 117D RVQ @PHHUAF 106 uF qF hyllD heoretil ghemistry eounts 112D RHQ!RHW @PHHRAF 107 uF qF hyll nd eF F qomesD heoretil ghemistry eounts 125D WU!IHH @PHIHAF 108 uF qF hyllD heoretil ghemistry eounts 129D THQ!TIQ @PHIIAF 109 uF qF hyllD heoretil ghemistry eounts 135D IPV @PHITAF 110 uF qF hyllD privte ommunitionF 111 uF qF hyllD heoretil ghemistry eounts 108D QQS!QRH @PHHPAF 112 uF qF hyllD heoretil ghemistry eounts 115D RRI!RRU @PHHTAF 113 uF qF hyllD heoretil ghemistry eounts 131D I!PH @PHIPAF 114 vF issherD F tF veeD nd uF qF hyllD tournl of ghemil hysis 105D VUTW!VUUT @IWWTAF 115 F pF foys nd pF fernrdiD woleulr hysis 19D SSQ @IWUHAF 116 tF qrnt rillD eF yF witrushhenkovD nd uF eF etersonD he tournl of hemil physis 138D IQRQIR @PHIQAF 117 tF enekowitshD Spektroskopische Eigenschaften aus elektronischen WellenfunktionenD hFhF thesisD niversität prnkfurtD qermny @IWVVAD he ps progrm is ville upon reE quest from uirk eterson @shington tteAF 118 uF qF hyllD he tournl of hemil physis 100D PIIV @IWWRAF 119 vF issher nd F ueD tF ghemF hysF 113D QWWT @PHHHAF 120 F tF tephensD pF tF hevlinD gF pF ghlowskiD nd wF tF prishD tF hysF ghemF 98D IITPQ @IWWRAF 121 eF hF fekeD he tournl of ghemil hysis 98D IQUP @IWWQAF 122 tF wF vévyEvelondD gommunF wthF hysF 6D PVT @IWTUAF

Terms S E α (Zα) 0 1 0 α (Zα) 1 0 0 α (Zα) 2 < 0 < 0 . . . . .

  Calculations within the rigorous QED framework have been reported for few-electron systems, and they are in excellent agreement with experiment. Examples are the Lamb shift of Li-like uranium, 18,19 the hyperne cou- pling constant (HFCC) of few-electron atoms, 20,21 and the anomalous g factor, 22 that provides stringent tests of the accuracy of QED.

  , in Ref. 41, the QED eect was found to shorten the bond length of TsH + , LvH and OgH + . The reason for this opposite trend may be that the valence orbitals have p-rather than s-orbital contributions from the heavy atom.

  proximate one-component relativistic Hamiltonians without picture-change. 4547 To include QED eects in a more rigorous manner, it seems more appropriate to include eective QED potentials in 4-component relativistic allelectron calculations. In this work, we report the implementation of eective QED potentials in the DIRAC code for relativistic molecular calculations. 48 Three potentials have been im- plemented: the Uehling potential 23 for vacuum polariza- tion, Pyykkö and Zhao's model SE potential, 24 as well as the eective SE potential of Flambaum and Ginges. 25 Our implementation is based on numerical routines from the GRASP atomic code 28 that have been grafted onto the DFT grid of DIRAC. 49 As rst molecular applications of our implementation we have chosen three case studies: the AuCN molecule for which Pekka Pyykkö has suggested QED eects on the bond length. 50 the van der Waals dimers M 2 (M = Hg, Rn, Cn, Og) for which one might suspect QED eects to be on par with interaction energies. Interestingly, van der Waals forces have been described in terms of vacuum uctuations. 51,52 the reaction energy of Pb hydrides, PbH 4 --→ PbH 2 + H 2 , suggested by Dyall et al. as a possible candidate for a signicant QED eect in chemistry.

  : two of them give state-independent energy-shifts and are usually ignored within a perturbative setting, whereas the remaining three represent electron self-energy, vacuum polarization and single-photon exchange. The latter diagram describes the relativistic electron-electron interaction, mediated by photons, to lowest order and is in line with the statement of Dirac: Classical electrodynamics, in its accurate (restricted) relativistic form, teaches us that the idea of an interaction energy between particles is only an approximation and should be replaced by the idea of each particle emitting waves, which travel outward with anite velocity and inuence the other particles in passing over them. 64 In the diagrams of Fig. 1 double electron lines appear to indicate that we are working within the Bound-State QED (BSQED) framework in which the Dirac eld operators are expanded in solutions of the Dirac equation in some external (contravariant) four-potential: A e = (φ e /c, A e ) (Furry picture 65

  )(see alsoRefs. 67 and 68). This potential is named after Uehling who rst calculated it in 1935 for a point charge nuclear distribution (as indicated by the superscript point). The corresponding potential for an arbitrary nuclear distribution ρ nuc. , normalized to one, is obtained by the following convolution 66 φ nuc.Ueh. (x) = d 3 y ρ nuc. (y) φ point Ueh. (x -y) .

  (Zα) processes, are briey discussed at the end of Sec. A 6 b.

  ) in terms of the variable x = (Z -80) α and a cutothe contribution of φ point elec at short distances where the the locality of the eective SE potential breaks down. The coecients of the A and B tting functions above were adjusted to reproduce the self-energy corrections to high sand p-states, respectively, calculated accurately in Refs. 87,88 for Coulombic hydrogen-like atoms of 5 ≤ Z ≤ 110. It should be added that Thierfelder and Schwerdtfeger 89 later modied the tting function to A n (Z, x), that is, making it dependent of the principal quantum number n. These potentials with A n (Z, x) instead of A (Z, x) were used by Pa²teka et al. to calculate the electron anity and ionization potential of gold. 32 Ginges and Berengut, 84 on the other hand, made both tting functions A and B dependent on orbital angular momentum ℓ and further suggest to introduce a κ-dependence as well. The downside of making the eective QED potentials dependent on atomic orbital quantum numbers is that it complicates the extension of these potentials to the molecular regime. This is in fact what motivated us to implement the eective SE potentials of Pyykkö/Zhao 24 and Flambaum/Ginges, 25 rather the one proposed by Shabaev and co-workers. 26,27,29 C. Atomic shift operator With the above eective QED potentials available in an atomic code (see Section III), we have investigated their extension to molecular calculations by adding to the electronic Hamiltonian, Eq. (

  where {ω i } are expectation values of the eective QED potentials taken from atomic calculations and {ψ i } are pre-calculated atomic orbitals, in practice limited to those that are occupied in the electronic ground state of the atoms constituting the molecule under study, calculated in their proper basis. The import of atomic orbitals into molecular calculations is straightforward in the case of the DIRAC code, since such functionality is already available through projection analysis. 90,91 There is some overlap between the spectral representation of the self-energy proposed by Dyall 92 as well as the eective SE operator proposed by Shabaev and co-workers, 26 but those approaches are based on hydrogenic orbitals. III. IMPLEMENTATION Routines for the radiative potentials used in this work are available in the GRASP atomic code. 28 Routines for calculating the Uehling potential were reported as early as 1980. 93 McKenzie et al. follow the approach suggested by Wayne Fullerton and Rinker. 66 More precisely, they employ Eq. (11) for the inner grid points until a more approximate form, Eq. (6) of Ref. 66, becomes numerically valid. The latter form is then used until the magnitude of the potential falls below a threshold value. The eective SE potential of Flambaum and Ginges 25 was imple- mented more recently, 89 as is also the case 94 of the eec- tive SE potential of Pyykkö and Zhao. 24 As already men- tioned, the FG potential is in principle that associated with a point nucleus, although tting parameters have been optimized also to calculations with nite nuclear charge distributions. Thierfelder and Schwerdtfeger 89 adapted these potentials to nite nuclei by replacing the Coulomb potentials of Eqs. (

For

  all calculations we used a development version of DIRAC code; 48,98 precise version and build informa- tion is found in output les, see Ref. 99. A Gaussian model 100 for the nuclear charge distribution was em- ployed throughout our calculations. Unless otherwise stated, we applied the Uehling VP potential 23 and the SE potential of Flambaum and Ginges, 25 added to the DiracCoulombGaunt (DCG) Hamiltonian. For correlated calculations we employed the molecular mean-eld approximation Hamiltonian (X2Cmmf ) 101 based on the DCG Hamiltonian, which we denote as 2 DCG M . In this approach, the converged Fock matrix obtained with the DCG Hamiltonian, with the eective QED potentials included, is exactly transformed to two-component form, that is, without any picture-change errors. 4547 All basis sets were employed in uncontracted form with the small component generated by restricted kinetic balance (see Ref. 48 for details). Electronic structure analysis was carwas done at the molecular geometries optimized with respect to the employed Hamiltonian, except for DCG with eQED, where the DCG structures were employed.

  forbid the transition of positive-energy electrons to the negative-energy continuum by the Pauli exclusion principle, and obtain a stable atomic theory, Dirac 140 postulated that this continuum should be to- tally lled with electrons that are not observed (Dirac sea). This means that the vacuum state is redened to be the state containing no positive-energy electrons and a fully occupied negative-energy electron sea. Dirac then argued that when a negative-energy electron absorbs enough energy (E ≥ 2m e c 2 ), it becomes real (observable), and leaves, for mass-and charge-conservation reasons, a positron behind (Dirac hole theory). 141 This last reasoning allows one to dene 142

)

  Despite its experimental success in predicting the existence of the positron, 143 the hole theory (its physical implications) was, shortly after its introduction, abandoned. Many physicists including Pauli, Bohr, Weisskopf, Heisenberg and Majorana, opposed this theory, as clearly indicated in Refs. 144 (section 1.6), 145, 146 (section 4.4) and 147. This opposition came mainly from the following aws of the Dirac hole theory: 1) the existence of a non-observable innite negative energy and charge and 2) for massive boson systems, whose wavefunctions satisfy the Klein-Gordon equation, the Dirac argument would not hold, and the existence of these bosons is not justied. Modern quantum eld theory reached the same mathematical expressions derived with respect to Dirac's hole theory, but provided a more symmetric picture between electrons and positrons, in which 1) one only sees electrons and positrons with positive energies, 2) the innite negative-energy electron sea assumption is no longer necessary, and 3) operators such as the Hamiltonian, and charge are replaced by their normal-ordered forms. This physical interpretation leads to the modern denition of the vacuum state, that obeys a i 0 e = b i 0 e = 0, ∀i, (A8) and contains zero positive-energy electrons and positrons. To get a wider and more detailed vision of the historical development of the quantum eld theory, the reader may consult Weinberg in Ref. 148 (section 1.2 and chapter 5) and Ref. 149, Mehra in Ref. 150 (chapter 29), Schweber in Ref. 144, Kragh in Ref. 151 and Weisskopf in Ref. 152. and the polarization vectors satisfy the following completeness relation 3 r=0 ζ r ϵ µ (k, r) ϵ ν (k, r) = -g µν .

  and discussions on the photon Hamiltonian, the reader may consult Greiner and Reinhardt in Ref. 61 (section 7.3) and Mandl and Shaw in Ref. 78 (chapter 5).

  tions, and is represented by a Feynman diagram in Fig. 1. The elements of these diagrams are the following: 1. Double external lines represent bound-electrons, i.e., with wavefunctions and energies satisfying the interacting Dirac equation of Eq. (A4), in the presence of a classical time-independent external potential A e (x).

2.

  Double internal lines represent a virtual boundelectron propagation between the two vertices S F A e (x 2 , x 1 ) and arise from a single contraction of two electron eld operators. 3. Internal wiggly-lines connecting two vertices represent propagations of virtual-photons D F µ2µ1 (x 2 , x 1 ) and come from a single contraction of two photon eld operators. The last two contributions from Eq. (A50) correspond to fully contracted products, and they are thus free of creation and annihilation operators. This means that their corresponding energy-shifts ∆E α,2 D1 and ∆E α,2 D2 are state-independent and hence do not contribute to energydierences. They are therefore discarded from further consideration; see for instance Mohr in Ref. 165. On the other hand, the rst three contributions correspond to partially contracted products, associated with the follwing physical processes: a. SP: Single-photon exchange This process, coming from the SP term in Eq. (A50)

  A52) as noted by Mohr in Ref. 166 (section IV). Notice that for µ = 0 , and µ = 1, 2, 3 these integrals account for the Coulomb and Gaunt interaction, respectively. On the other hand, if we used the Coulomb gauge photon propagator instead of the Feynman one, we would get the retarded Breit interaction, as noted by Lindgren in Ref. 167 (page 262). b. VP: Vacuum polarization This process, presented in Fig. 1b, accounts for the instantaneous interaction of a bound-electron with the electron-positron pair cloud, polarized by the presence of a classical potential source. After plugging the VP term of Eq. (A50) in the second-order scattering matrix expression, one can use Sucher's formula of Eq. (A45) to write the energy-shift associated with the vacuum polarization process as Ref. 59 (chapter 15 eq.(

  A55) as mentioned in Ref. 168. The isolation of the divergent terms in this expression can be done by expanding the propagator inside the trace using Eq. (A41), and write the energy as∆E α,2 VP = ∆E α,2 VP,0 + ∆E α,2 VP,1 + ∆E α,2 VP,2 + . . . (A56)where ∆E α,2VP,i represents the term that corresponds to an i number of interactions with the external potential (Zα) i . The rst four terms are presented in Figs. 3a to 3d. Notice that the double-line loop is replaced by a single-line one. This is made to indicate that the propagators between these vertices are the free ones S F 0 , in- stead of the bound-ones S F A e . Using Furry's theorem, 169 that is based on a charge conjugation symmetry argument, one can show that any diagram containing a freeelectron loop with an odd number of vertices does not contribute. This means that the above energy expression reduces to ∆E α,2 VP = ∆E α,2 VP,1 + ∆E α,2 VP,3 + . . . (A57) A naive estimation of the degree of divergence of a QED integral can be done by calculating the supercial degree of divergence S that simply counts overall momentum powers of the integral in question (in momentumspace): S ≡ 4 -N e -2N p , (A58) where 4 are the spacetime dimensions and N e and N p are the number of electron and photon propagators, respectively, in the loop in question; see for instance Refs. 60 (section 10.1) and 86 (sections 7-1-4 and 8-1-3). The integral is said to be supercially divergent if S ≥ 0. The possible cases are: S ≤ 0 convergence S = 0 logarithmic divergence S = 1 linear divergence S = 2 quadratic divergence (A59)

  ble XVII and noted by Gyulassy, 171 this contribution is free of divergences. Wichmann and Kroll calculated the eective potential associated with the ∆E α,2 VP,3 cor- rection in Laplace space. On the other hand, in Ref. 172 (section 4), Blomqvist has evaluated the inverse Laplacetransform, and obtained the real-space potential expression for a point nuclear charge distribution. The last reference presents a relatively complex analytical expression for this α (Zα) 3 potential, and this motivated Fainshtein et al. 173 to provide an approximation that facilitates the numerical computation, yet conserving precision. We nally note that in the fourth-order BSQED correction, one nds the KällénSabry potentials 71 of order α 2 (Zα) that can be obtained by expanding the bound propagators of Ref. 162 (g. 25 diagrams b,c VPVP). In order to make this momentum-space potential usable in practical calculations, in Ref. 172 (section 3) Blomqvist derived its real-space version for a point nucleus distribution, whereas Wayne Fullerton and Rinker generalized this result to account for an extended nuclear charge distribution; see Ref. 66 (eq.(

  c. SE: Self-energyThe self-energy process, presented in Fig.1c, is the dominant radiative QED correction in electronic atoms, as seen in the work of Johnson and So ofRef. 123 (g.2).This process describes the interaction of the bound-electron with itself, by emitting and absorbing a virtual-photon. The rst calculation for this correction was made in 1947 by Hans Bethe in a purely nonrelativistic framework, 174 where he used a renormaliza- tion technique (by subtracting the free self-energy) to render the integral less divergent, and introduced a reasonable virtual-photon energy cuto at E = m e c 2 . This simple calculation gave hope in digging for the physical Lamb shift in the frustrating non-physical divergences in the QED theory. Using Sucher's energy formula of Eq. (A45), the SE term of Eq. (A50) leads to the following energy-shift

  the last expansion are represented in Figs. 4a to 4d. The zero-and one-potential terms: ∆E α,2 SE,0 and ∆E α,2

  123 F tohnson nd qF o'D etomi ht nd xuler ht les 33D RHS @IWVSAF 124 tF hyssenD Development and Applications of Methods for Correlated Relativistic Calculations of Molecular PropertiesD hFhF thesisD niversity of outhern henmrk @PHHIAF 125 sn these lultions the hseg keyword yixpegy ws set to oneD suh tht oritl eigenvlues stis(es uoopmns9 theoremF 126 F eF erokhin nd F wF hevD tournl of hysil nd ghemil eferene ht 44D HQQIHQ @PHISAF 127 vF vzowskyD sF qoidenkoD wF okmnD nd F yykköD hysF evF e 59D PUHU @IWWWAF 128 F leskiEijgierdD wF tzshkeD nd F yykköD he tournl of hemil physis 128D PPRQHQ @PHHVAF 129 hF fF qrotjhnD wF eF frewsterD nd vF wF iurysD tournl of the emerin ghemil oiety 124D SVWS @PHHPAF 130 F ykyshiD iF F ykyshiD pF uotoD F sshidD nd wF nimotoD tournl of the emerin ghemil oiety 131D IIUIP @PHHWAF 131 hF piggenD qF uhutD wF holgD nd rF tollD ghemil hysis

  where the initial t 0 and nal times t are at ∓∞, to ensure Lorentz invariance. Upon expansion of the Ŝ-matrix operator in the fundamental charge e, the nth-order term Ŝ(n) contains a time-ordered string of n interaction Hamiltonian densities H I , as seen

	in Eq. (A29). Using Wick's theorem, 62 a time-ordered
	string is converted into a linear combination of normal-

Dirac Lagrangian density, in accordance with the principle of minimal electromagnetic interaction (term coined by GellMann 58 ). For detailed derivations and discus- sions, the reader may consult Schweber in Ref. 59 (chapter 10), Peskin and Schroeder in Ref. 60 (chapter 4), as well as Greiner and Reinhardt in Ref. 61 (section 8.6).

The scattering matrix is the special case of the timeevolution operator Û (t, t 0 ),

  and F 1 and F 2 are known as the

	term form factor comes from diraction physics; see for
	instance Ref. 80.
	electric and magnetic form-factors, respectively (corre-
	sponding to f -1 and g in Eq. (116.6) of Ref. 79). The

  We used an virtual energy cuto of 40 E h . Dyall cv3z basis sets, 106108 designed for core-valence correlation, were employed for the Hg and Cn species, whereas Dyall acv3z basis sets, 111113 where the Dyall cv3z basis sets have been augmented by diuse functions, were employed for Rn and Og species. Electronic structure calculations were done at the level of coupled-

Table I we employed Dyall v3z basis sets; 104109 the basis set for Uue was specially optimized by Dyall for this work. 110

For van der Waals dimers, the following orbitals were correlated: 5d6s for Hg, 5d6s6p for Rn, 6d7s for Cn, and 6d7s7p for Og.

cluster singles-and-doubles with approximate triples correction (CCSD(T)) using the RELCCSD module.

  for Au, and all electrons of C and N were correlated, which is the same level as the previous work. 116 Dyall ae3z and ae4z ba-

	sis sets, 106,107,109 designed for correlation of all electrons,
	were employed in the calculations. We employed a vir-
	tual energy cuto of about 50 E h and 80 E h for dyall.ae3z
	and dyall.ae4z, respectively, which assures that correlat-
	or better on the gradient.	In addition, to estimate
	the relativistic eects we employed the two-component
	non-relativistic (by using .NONREL keyword), 4c-scalar-
	relativistic, 118,119 and the DiracCoulomb Hamiltonians
	at the density functional theory (DFT) level. In these cal-
	culations we employed the B3LYP functional 120,121 and
	the dyall.3zp basis sets. 106,107,109

ing h and i orbitals, respectively, are included. Eective QED potentials for C and N atoms were not used, as explained in Section III. The potential energy surface (PES) was calculated in the vicinity of the equilibrium structure with a total of 49 points for each basis set, using internal coordinates r 1 (Au-C distance), r 2 (CN distance). The bond angle was xed at θ = 180 • . The step size for bond distances was 0.1 a 0 . The surface tting and determination of the equilibrium structure was carried out using the SURFIT program, 117 with convergence 3.2 × 10 -10

For the calculation of Pb and Fl hydrides, the DCG Hamiltonian with and without eective QED potentials, as well as the Lévy-Leblond (LL) 119,122 Hamiltonian were employed. The dyall.3zp basis sets were used for all of the elements. The B3LYP functional was used for both the projection analysis and the geometry optimization.

TABLE I .

 I Calculated ns orbital energies in eV of group 1 and 11 elements from AOC-HF/dyall.v3z calculations based on the NR and DC Hamiltonians. The VP (Uehling) and SE (FlambaumGinges) corrections have been calculated as expectation values.

		NR	DC	VP	SE	∆QED	SE/VP	SE/VP 123	∆QED/∆R[%]
	Li	-5.342	-5.343	-1.373E-06	4.092E-05	3.955E-05	-29.7949	-29.7058	-9.01
	Na	-4.955	-4.962	-1.536E-05	2.950E-04	2.796E-04	-19.2057	-18.7963	-4.37
	K	-4.013	-4.028	-3.423E-05	5.155E-04	4.813E-04	-15.0615	-14.7030	-3.19
	Rb	-3.752	-3.811	-1.309E-04	1.361E-03	1.231E-03	-10.3981	-10.0783	-2.08
	Cs	-3.365	-3.490	-2.989E-04	2.304E-03	2.005E-03	-7.7089	-7.4266	-1.61
	Fr	-1.740	-3.611	-1.438E-03	6.333E-03	4.895E-03	-4.4038	-4.3351	-0.26
	Uue	-2.993	-4.327	-1.034E-02	2.157E-02	1.123E-02	-2.0859	-4.3351	-0.84
	Cu	-6.480	-6.649	-2.355E-04	2.840E-03	2.604E-03	-12.0606	-11.7316	-1.54
	Ag	-5.985	-6.452	-7.342E-04	6.448E-03	5.714E-03	-8.7825	-8.4755	-1.22
	Au	-6.003	-7.923	-4.635E-03	2.374E-02	1.910E-02	-5.1219	-4.9912	-1.00
	Rg	-5.441	-11.425	-3.251E-02	8.408E-02	5.157E-02	-2.5863	-2.7223	-0.86

TABLE II .

 II Relativistic and QED eects on the orbital energies ε(E h ) of the Au atom at the B3LYP/dyall.3zp level. The Uehling VP potential has been combined with two dierent SE potentials: FG (FlambaumGinges) and PZ (PyykköZhao) in variational calculations. Numbers in parentheses shows the percentage-wise ratio ∆QED/∆R for each combination of eective QED potentials.

		NR	DCG	∆(U+FG)		∆(U+PZ)
	1s 1/2	-2689.451	-2955.841	6.377E00	(-2.39)	6.243E00	(-2.34)
	2s 1/2	-449.932	-523.020	8.448E-01	(-1.16)	8.586E-01	(-1.17)
	2p 1/2	-432.492	-500.523	5.895E-02	(-0.09)	5.955E-02	(-0.09)
	2p 3/2	-432.492	-433.755	1.194E-01	(-9.45)	-3.055E-02	( 2.42)
	3s 1/2	-105.753	-123.999	1.862E-01	(-1.02)	1.924E-01	(-1.05)
	3p 1/2	-97.515	-114.096	7.320E-03	(-0.04)	1.462E-02	(-0.09)
	3p 3/2	-97.515	-99.311	2.213E-02	(-1.23)	-8.516E-03	( 0.47)
	3d 3/2	-82.165	-83.236	-1.430E-02	( 1.34)	-9.090E-03	( 0.85)
	3d 5/2	-82.165	-80.090	-5.817E-03	(-0.28)	-8.587E-03	(-0.41)
	4s 1/2	-22.559	-27.178	4.650E-02	(-1.01)	4.833E-02	(-1.05)
	4p 1/2	-18.993	-22.979	1.124E-03	(-0.03)	3.336E-03	(-0.08)
	4p 3/2	-18.993	-19.424	4.651E-03	(-1.08)	-2.427E-03	( 0.56)
	4d 3/2	-12.440	-12.630	-3.439E-03	( 1.82)	-2.307E-03	( 1.22)
	4d 5/2	-12.440	-11.971	-1.648E-03	(-0.35)	-2.185E-03	(-0.47)
	4f 5/2	-3.648	-3.228	-2.383E-03	(-0.57)	-1.461E-03	(-0.35)
	4f 7/2	-3.648	-3.091	-1.824E-03	(-0.33)	-1.425E-03	(-0.26)
	5s 1/2	-3.253	-4.116	9.003E-03	(-1.04)	9.430E-03	(-1.09)
	5p 1/2	-2.108	-2.745	-4.569E-05	( 0.01)	4.105E-04	(-0.06)
	5p 3/2	-2.108	-2.139	5.619E-04	(-1.82)	-5.703E-04	( 1.85)
	5d 3/2	-0.346	-0.333	-4.662E-04	(-3.61)	-3.486E-04	(-2.70)
	5d 5/2	-0.346	-0.276	-2.881E-04	(-0.41)	-3.206E-04	(-0.46)
	6s 1/2	-0.148	-0.205	6.519E-04	(-1.14)	6.726E-04	(-1.18)
	come with opposite sign and are dominated by the
	latter. 123 However, the ratio SE/VP decreases signi-
	cantly with increasing nuclear charge and, indeed, VP is predicted to eventually overtake SE at very high nu-clear charges. 89 Pyykkö et al. calculated the SE contri-bution as ⟨V U ⟩ * (SE/VP) where (SE/VP) is the ratio for 2s 1/2 of the corresponding hydrogen-like systems, includ-ing the nuclear-size eect, tabulated for 1 ≤ Z ≤ 100 by Johnson and So 123 II we show the eect of relativity and QED on all orbital energies of the gold atom. Two combinations of eective QED potentials have been used in variational calcula-tions: The Uehling (U) VP potential has been combined either with the FlambaumGinges (FG) or PyykköZhao SE potentials. One sees that for both combinations of ef-

(a more recent compilation is pro- vided by Yerokhin and Shabaev 126 ). As conrmed by later calculations 127 and the numbers in Table

I

, this is a quite reasonable approximation.

Comparing relativistic and QED eects, one sees that the latter corrects the former by about -1% for the heavier atoms. For the gold atom it is exactly so. In Table

fective

QED potentials the relativistic eects is, with very few exceptions, reduced with a few percent. For s 1/2 orbitals the dierence in QED shift between the U+FG and U+PZ combinations is below 5%; for other orbitals the dierence is generally larger. We note in particu-

TABLE III .

 III First-order QED eects on the orbital energies ε(E h ) of the Au atom at the B3LYP/dyall.3zp level using the DCG Hamiltonian or the X2C Hamiltonian, the latter with or without picture-change (PC) transformation. The Uehling VP potential has been combined with two dierent SE potentials: FG (FlambaumGinges) and PZ (PyykköZhao). Numbers in parentheses shows the percentage-wise ratio ∆QED/∆ for each combination of eective QED potentials. orbitals; the very largest relative deviation concerns 5p 1/2 , but this can probably be attributed to noise, since the QED shift on the energy of this orbital is particularly small with both combinations of eective QED potentials.

		DCG		X2C -PC	X2C -noPC
		∆(U+FG)	∆(U+PZ)	∆(U+FG)	∆(U+PZ)	∆(U+FG)	∆(U+PZ)
	1s 1/2	6.620E00	6.465E00	6.632E00	6.477E00	7.522E00	8.253E00
	2s 1/2	9.003E-01	9.022E-01	9.012E-01	9.030E-01	1.079E00	1.175E00
	2p 1/2	1.184E-01	1.051E-01	1.187E-01	1.047E-01	1.971E-01	3.608E-02
	2p 3/2	1.671E-01	4.922E-03	1.682E-01	5.005E-03	1.279E-01	6.342E-03
	3s 1/2	2.016E-01	2.040E-01	2.020E-01	2.043E-01	2.441E-01	2.665E-01
	3p 1/2	2.357E-02	2.674E-02	3.574E-02	1.357E-03	2.671E-02	1.727E-03
	3p 3/2	3.554E-02	1.335E-03	3.574E-02	1.357E-03	2.671E-02	1.727E-03
	3d 3/2	-1.234E-03	1.064E-05	-1.220E-03	1.066E-05	3.906E-03	3.605E-07
	3d 5/2	6.494E-03	-1.680E-07	6.528E-03	-1.675E-07	3.296E-03	-1.132E-07
	4s 1/2	5.082E-02	5.151E-02	5.097E-02	5.166E-02	6.173E-02	6.741E-02
	4p 1/2	5.631E-03	6.658E-03	5.638E-03	6.637E-03	9.363E-03	2.339E-03
	4p 3/2	8.392E-03	3.315E-04	8.444E-03	3.370E-04	6.292E-03	4.295E-04
	4d 3/2	-1.331E-04	2.846E-06	-1.299E-04	2.851E-06	9.447E-04	9.962E-08
	4d 5/2	1.461E-03	-4.294E-08	1.468E-03	-4.277E-08	7.944E-04	-2.822E-08
	4f 5/2	-2.785E-04	-1.968E-10	-2.782E-04	-1.974E-10	1.089E-05	-1.184E-10
	4f 7/2	2.190E-04	-9.047E-11	2.194E-04	-9.111E-11	9.735E-06	-9.557E-11
	5s 1/2	1.001E-02	1.015E-02	1.004E-02	1.018E-02	1.217E-02	1.329E-02
	5p 1/2	9.658E-04	1.152E-03	9.668E-04	1.149E-03	1.605E-03	4.052E-04
	5p 3/2	1.373E-03	5.478E-05	1.382E-03	5.572E-05	1.029E-03	7.103E-05
	5d 3/2	-1.078E-05	2.580E-07	-1.048E-05	2.582E-07	8.410E-05	9.077E-09
	5d 5/2	1.216E-04	-3.635E-09	1.222E-04	-3.618E-09	6.669E-05	-2.375E-09
	6s 1/2	7.888E-04	8.001E-04	7.908E-04	8.020E-04	9.586E-04	1.047E-03
	Table III also shows perturbative QED shifts of or-
	bital energies obtained with the X2C Hamiltonian. When
	the eective QED potentials have been correctly picture-
	changed transformed, deviations from the parent 4c cal-
	culation are below 3 %, which clearly validates the use
	of these potentials in 2-component relativistic calcula-
	tions. On the other hand, without picture-change, sig-
	nicant errors are observed; the average unsigned error

lar that the shifts have systematically opposite sign for p 3/2 orbitals. Not surprisingly the largest absolute shifts are observed for inner core orbitals, whereas the largest relative shift 0.33% is seen for the 6s 1/2 orbital.

In Table

III we

show QED shifts of orbital energies, this time obtained perturbatively as expectation values. Compared to the shifts obtained from variational inclusion of the eective QED potentials, the largest absolute deviations concern the inner core orbitals. The smallest relative deviations are observed for s 1/2 orbitals and decreasing towards core. The largest relative deviations, on the other hand, are seen for p for U+FG and U+PZ is 130 % and 47 %, respectively. This is possibly worrisome since the U+PZ combination, expressed in terms of Gaussians, have been used without picture-change in scalar DKH calculations by Peterson and co-workers.

3336

TABLE IV

 IV

		CuCN	AgCN	AuCN
	r0	183.231(7) 203.324(45)	191.251(16)	MW a
	rs	183.284(4) 203.4182(27) 191.22519(84) MW a
	re	182.36	202.42	191.05	Calc. 128
	re	182.65	202.99	190.71	Calc. 116
		a gugxX efF IPWF AgCN,AuCN: Ref. 130
	B. Gold cyanide
		In 2008 Pyykkö and co-workers reported CCSD(T)/cc-
	pVQZ calculations on the noble metal cyanides (MCN,
	M=Cu, Ag, Au). 128 Small-core scalar-relativistic eec-
	tive core potentials (SRECP) 131 were used for the metal
	atoms and spin-orbit corrections added at the PBE-
	ZORA/QZ4P level. In 2013 Peterson and co-workers re-
	ported CCSD(T)-F12/cc-pV5Z calculations on the same
	compounds, using the same SRECPs as the previous au-
	thors and adding a number of corrections. 116 As seen
	from Table IV the newer calculations brought the M-C
	bond lengths of CuCN and AgCN in better agreement

. M-C bond lengths (in pm) in MCN (M=Cu, Ag, Au) from microwave (MW) spectroscopy and calculations.

with experiment, but increased the gap between theory and experiment for AuCN. This led Pyykkö to conjecture that this could be the rst evidence of the eect of QED on molecular structure.

50

To possibly verify this conjecture we rst carried out

TABLE V .

 V Relativistic and QED eects on the rms radius ⟨r 2 ⟩ 1/2 (pm) of the Au atom at the B3LYP/dyall.3zp level.We next turn to the AuCN molecule. We rst, in TableVI, report bonding analysis in localized orbitals. One may note that the QED eect is of the same order as the eect of adding the Gaunt term 89In passing we note from Table VII that incorporation of QED eects through the atomic shift operator (ASHIFT) described in Section II C also leads to a bond extension, albeit only capturing half of the full QED eect. We then added the triples ∆T and quadruples ∆Q corrections reported by Hill et al. 116 to obtain a

	Eective QED potentials: VP(Uehling)+SE(Flambaum
	Ginges).
	NR DCG DCG+QED ∆R ∆QED ∆QED/∆R(%)
	5s 1/2 57.71 52.23	52.28	-5.48 0.05	-0.86
	5p 1/2 63.16 56.78	56.78	-6.38 0.00	-0.02
	5p 3/2 63.16 62.37	62.38	-0.80 0.01	-1.04
	5d 3/2 91.07 90.81	90.78	-0.26 -0.03	10.46
	5d 5/2 91.07 95.75	95.73	4.68 -0.02	-0.36
	6s 1/2 196.07 167.54	167.79	-28.53 0.25	-0.88
	exploratory calculations at the B3LYP/dyall.3zp level.
	Table V shows the eects of relativity and QED on or-
	bital sizes of the gold atom. For the valence 6s 1/2 we ob-
	serve an impressive relativistic contraction of 28.53 pm,
	whereas QED leads to an orbital expansion of 0.25 pm,
	roughly -1 % of the relativistic eect.
					116 How-
	ever, this contraction is almost canceled when adding the
	Gaunt term, which brings spin-other-orbit interaction 3
	and which was not considered by Hill and co-workers. 116
	To obtain more accurate bond lengths, we proceeded as
	indicated in Table VIII:	2 DCG	M -CCSD(T) calculations
	were carried out in the Dyall ae3z and ae4z basis sets
	and then extrapolated to the basis-set limit, 132 indicated
	by ae∞z.

91 One nds a single σ-type Au-C bond, dominated by carbon 2s 1/2 and gold 6s 1/2 , as well as a triple C-N bond. Equi- librium bond lengths r e with respect to dierent Hamiltonians are reported in Table

VII

. We see a very signicant scalar relativistic bond-length contraction of 25.31 pm, on par with the 6s 1/2 orbital contraction observed in Ta- ble V. When going from a spin-free (SF) Hamiltonian to the DiracCoulomb one, one nds a further contraction of 0.29 pm, which agrees very well with the spin-orbit correction of -0.28 pm obtained by Hill et al. taking the same dierence, albeit at the CCSD(T) level. At this level of theory, the total relativistic eect on the bond length is thereby -25.38 pm; in future work we hope to assert the eect of replacing the Gaunt term by the full Breit one. Finally, adding QED eects, we observe a bond-length extension of 0.19 pm, -0.75 % of the relativistic eect. Au-C bond length of 190.75 pm, very close to the value 190.71 pm reported by Peterson and co-workers. Finally, we add a QED correction of 0.19 pm, identical to what we obtained at the B3LYP/dyall.3zp level, to obtain our nal value of 190.99 pm.

TABLE VI .

 VI Gross populations obtained by projection analysis using PipekMezey localized orbitals at the DCG/B3LYP/dyall.3zp level. ⟨ε⟩ is the expectation value in E h with respect to the KohnSham matrix.

TABLE VII .

 VII Equilibrium bond lengths re (in pm) of AuCN calculated at the B3LYP/dyall.3zp level using various Hamiltonians. Numbers in parenthesis indicate the change with respect to the previous line, except ASHIFT, which refers to DCG. SF refers to a spin-free 4-component relativistic Hamiltonian.

	Hamiltonian	Au-C	C-N
	NR	218.54	115.71
	SF	193.23 (-25.31)	115.54 (-0.17)
	DC	192.94 (-0.29)	115.56 (+0.02)
	DCG	193.16 (+0.22)	115.58 (+0.01)
	QED	193.35 (+0.19)	115.57 (+0.00)
	ASHIFT	193.25 (+0.09)	115.58 (+0.00)

The devil is, however, in the details: Our Born Oppenheimer equilibrium bond lengths r e are not directly comparable to the structural parameters extracted from the rotational spectra recorded by Okabayashi and co-workers. 130 Experiment gives access to rotational con- stants B ν for individual vibrational states. For a linear molecule like AuCN the rotational constant, in units of frequency, is expressed as

TABLE VIII .

 VIII Final, recommended equilibrium Au-C bondlength re (pm) at the 2 DCG M -CCSD(T) level for the AuCN molecule. ∆QED is the dierence between the extrapolated basis set limit ae∞z with QED and without QED. target molecule observed in experiment. 133,134 Ef- fective structures r 0 are obtained by least-square tting of experimental ground-state inertial moments, whereas substitution structures r s are obtained from observation of how rotational constants (and center of mass) change upon single isotope substitution A → A ′ . For a linear molecule one has

		Au-C	C-N
	ae3z	190.89	116.66
	ae4z	190.70	116.28
	ae∞z	190.58	116.07
	∆T 116	0.26	-0.10
	∆Q 116	-0.09	0.19
	Final w/o QED	190.75	116.16
	∆QED	0.19	0.00
	Final	190.94	116.16

the

  the case of AuCN |z C | and |z N | could be estimated

	from corresponding single isotope substitutions. How-
	ever, since gold has a single naturally occurring isotope, 197 Au, |z Au | was obtained from the denition of center of mass. 135
	Empirically one typically nds r 0 ≥ r s ≥ r e , 134
	which suggests that we should rather compare our recom-
	mended r e =190.99 pm for Au-C with the corresponding
	substitution bond length r s =191.22519(84) pm reported by Okabayashi and co-workers. 130 However, a better com-
	parison is provided by calculating the ground-state rota-

tional constant B 0 from B e . From perturbation theory, excluding Fermi resonances, the rotational constant for a given vibrational state ν of a general molecule is related to B e as follows:

134

TABLE IX .

 IX Calculated and experimental rotational constants (in MHz) for AuCN. 197 Au 12 C 14 N 197 Au 13 C 14 N 197 Au 12 C 15 N

	α1	14.55	13.53	14.06
	α2	-10.98	-10.25	-10.59
	α3	12.17	11.91	11.30
			w/o QED
	Be	3237.5	3184.5	3086.4
	B0	3235.1	3182.1	3084.3
			with QED
	Be	3232.8	3179.9	3082.0
	B0	3230.4	3177.5	3079.9
	B0(exp.) 130 3230.21115(18) 3177.20793(13) 3079.73540(12)
	TABLE X. Calculated and experimental substitution struc-
	ture (in pm) for AuCN
		rs(Au-C)		rs(C-N)
	w/o QED	190.991		115.910
	with QED	191.184		115.909
	Exp.	191.22519(84)	115.86545(97)
	the potential surface to improve agreement with experi-
	ment; we note for instance that our calculated vibration-
	rotation interaction constant associated with bending for the most abundant isotopomer 197 Au 12 C 14 N is -10.98
	MHz, compared to -11.9781 MHz when extracted from
	experiment.

Table

  XI we report our calculated the equilibrium bond lengths r e , harmonic frequencies ω e , anharmonic constants ω e x e and dissociation energies D e for these species.

	We see that the QED eect on bond length is on the
	order of 0.15 pm for row-6 dimers and approximately
	doubles when going to the superheavy elements; for Og 2
	the QED bond length extension is in line with what was
	reported by Hangele and Dolg using relativistic eective
	core potentials. 40 The QED eect on dissociation energies
	is rather small: on the order of 0.4 % for the superheavy
	dimers.
	D. Reaction energies: Pb and Fl hydrides
	As a nal case study, we consider the reaction energy

of XH 4 --→ XH 2 + H 2 to which Dyall et al. proposed that the Lamb shift could make a chemically signicant contribution.

53 Their argument was based on the observation that QED eects are most important for s orbitals, as seen in Tables

II and V

, and that this is

TABLE XI .

 XI Spectroscopic constants of heavy group 12 and 18 dimers obtained at the 2 DCG M -CCSD(T) level, using the U+FG combination of eective QED potentials and either dyall.cv3z (Hg 2 ,Cn 2 ) or dyall.acv3z (Rn 2 ,Og 2 ) basis sets. Numbers in parentheses indicate the QED eect. re/pm ωe/cm -1 ωexe/cm -1 De/cm -1

	Hg 2	385.71	16.65	0.232	277.7
		(0.15)	(-0.03)	(-0.002)	(-0.02)
	Cn 2	354.75	22.95	0.255	532.7
		(0.36)	(-0.11)	(0.001)	(-2.78)
	Rn 2	463.60	13.79	0.286	174.9
		(0.14)	(-0.02)	(-1.E-04)	(-0.41)
	Og 2	449.97	17.10	0.210	391.1
		(0.28)	(-0.04)	(0.001)	(-1.32)
	a reaction with a signicant change of the valence s
	population of a heavy element. We have investigated
	this at the B3LYP/dyall.3zp level and also included the
	corresponding reaction involving the heavier homologue
	erovium. Optimized equilibrium structures are given
	in Table XII. For the tetrahydrides we assumed T d sym-metry, in line with experiment 137 (PbH 4 ) and previous calculation 138 (FlH 4 ).
	Turning next to Tables XV and XVI we see that both
	reactions are endothermic at the non-relativistic level,
	but becomes clearly exothermic when adding relativity.
	For Pb QED reduces the relativistic eect by 1.25%. Its
	value is 0.32 kcal/mol, which is at the lower end of the
	perturbation estimate of Dyall et al. 53 For Fl QED re-
	duces the relativistic eect by 0.50%. Interestingly, its
	value is very close to that for the Pb reaction, despite Fl
	being a much heavier atom. The reason for this unex-
	pected result is the cancellation between the SE and VP
	eects. From Tables

To monitor valence s populations we carried out bonding analysis in PipekMezey localized MOs.

91,102 From

Table XIII, the change of the valence s population from XH 4 to XH 2 is 0.45 and 0.30 for Pb and Fl systems, re- spectively. From Table XIV, one sees that in the tetrahydrides the valence s population is contained in the four σ XH bonds. In contrast, in the dihydrides the two σ XH bonds are mediated by the valence p orbitals of the metals, and most of the valence s population is found in a non-bonding (nb) orbital.

TABLE XII .

 XII Optimized equilibrium structures of Pb and Fl hydrides at the B3LYP/dyall.3zp level based on the DCG Hamiltonian. re and αe refer to the X-H bond length (Å) and H-X-H angle (degree), respectively.

		PbH 2		FlH 2		PbH 4	FlH 4
		re	αe	re	αe	re	re
	NR	1.879	90.83	2.017	90.85	1.816	1.959
	DCG	1.845	91.18	1.920	93.35	1.756	1.825

TABLE XIII .

 XIII Charge Q and electronic congurations of Pb and Fl atoms in the title compounds obtained by projection analysis at the B3LYP/dyall.3zp level.

		Q	atomic conguration
	PbH 2	0.39	5d 4.00 3/2 5d 5.99 5/2 6s 1.86 1/2 6p 0.90 1/2 6p 0.86 3/2
	PbH 4	0.66	5d 3.99 3/2 5d 5.98 5/2 6s 1.41 1/2 6p 0.85 1/2 6p 1.10 3/2
	FlH 2	0.32	6d 3.99 3/2 6d 5.97 5/2 7s 1.91 1/2 7p 1.26 1/2 7p 0.55 3/2
	FlH 4	0.46	6d 3.98 3/2 6d 5.94 5/2 7s 1.61 1/2 7p 1.17 1/2 7p 0.84 3/2
	H 2 Po 2 is 2.38 %, although it depends on the choice of eective QED potentials. 139

TABLE XIV .

 XIV Gross population obtained by projection analysis of the localized bonding orbitals in the title compounds at the B3LYP level based on the DCG Hamiltonian. < ε > refers to the expectation value with respect to the converged KohnSham matrix (in E h ).

						X				Hi
			< ε >	5s 1/2	5d 3/2	5d 5/2	6s 1/2	6p 1/2	6p 3/2	1s 1/2
	PbH 4	σXH i	-0.4192		0.00	0.00	0.35	0.21	0.27	1.19
	PbH 2	σXH i	-0.3384	0.00	0.00	0.00	0.03	0.37	0.37	1.24
		nb	-1.0182	0.13	0.00	0.00	1.62	0.14	0.11	-0.03
			< ε >	6s 1/2	6d 3/2	6d 5/2	7s 1/2	7p 1/2	7p 3/2	1s 1/2
	FlH 4	σXH i	-0.4524		0.00	0.01	0.40	0.27	0.20	1.11
	FlH 2	σXH i	-0.3469	0.00	0.00	0.01	0.03	0.52	0.25	1.18
		nb	-1.2235	0.11	0.04	0.06	1.57	0.15	0.04	-0.03

TABLE XV .

 XV Reaction energy of PbH 4 --→ PbH 2 + H 2 (in kcal/mol). ∆DCG refers to the dierence between DCG and NR. Other ∆ refers to the dierence from the DCG value.

		QED eect	reac. energy	∆(kcal/mol)	∆(%)
	NR	none	16.47
	DCG	none	-8.99	-25.46
		VP	-9.09	-0.10	0.41
		SE	-8.56	0.42	-1.67
		VP+SE	-8.66	0.32	-1.27
		VP+SE(ASHIFT a )	-8.97	0.02	-0.07
		VP+SE(ASHIFT b )	-9.06	-0.08	0.31

a yuption of tomi frgment ws 6s 2 6p 2 b sing the ouptions of le sss through their anti-commutation relation

TABLE XVI .

 XVI Reaction energy of FlH 4 --→ FlH 2 + H 2 (in kcal/mol). ∆DCG refers to the dierence between DCG and NR. Other ∆ refers to the dierence from the DCG value.

		QED eect	reac. energy	∆(kcal/mol)	∆(%)
	NR	none	9.52
	DCG	none	-60.02	-69.54
		VP	-60.43	-0.41	0.59
		SE	-59.27	0.75	-1.08
		VP+SE	-59.67	0.35	-0.50
		VP+SE(ASHIFT a )	-60.10	-0.07	0.10
		VP+SE(ASHIFT b )	-60.23	-0.21	0.30
	a yuption of the tomi frgment ws 7s 2 7p 2
	b sing the ouptions of le sss
	and the second line indicates that the annihilation of a
	negative-energy electron with c i is equivalent to the cre-
	ation of a (positive-energy) hole (positron) with b † i . The
	electron eld operator of Eq. (A2) can be written, with
	respect to these denitions, as

Table XVII

 XVII

	Terms S E
	α (Zα) 1 2	0
	α (Zα) 3 0 < 0
	α (Zα) 5 < 0 < 0
	. . .	. . .	. . .
	TABLE XVII. Supercial and eective degrees of divergence
	for the bound-state vacuum polarization contributions.
	i.e., if the external vector potential A e (x) vanishes, cf.
	Greiner in Ref. 170 (eqs.(12.52-53)), then only the time-
	component potential φ 0 VP,1 survives, and one obtains the Uehling potential, 23 given in Eq.(8).
	For detailed discussions and derivations of the one-
	potential bound-state vacuum polarization correction,
	the reader may consult the calculation of Greiner and
	Reinhardt in Ref. 85 (section 5.2) where the authors used
	PauliVillars regularization, in addition to Peskin and
	Schroeder in Ref. 60 (section 7.5) who used dimensional-
	regularization to treat the occurring divergences; see also
	Mandl and Shaw in Ref. 78 (section 10.4), in addition
	to Schwartz in Ref. 163 (section 16.2.2). Contrary to
	the conventional momentum-space approach to evaluate
	QED corrections, Indelicato and Mohr in Ref. 168 (sec-
	tion B) considered the vacuum polarization problem in
	coordinate space, and derived the physical Uehling con-
	tribution using coordinate-space PauliVillars regulariza-
	tion.
	The second non-vanishing vacuum polarization eect,
	associated with ∆E α,2
			, this term is of supercial quadratic divergence, but
			it is, eectively, only logarithmic. This can be shown
			using the Ward identity, as mentioned by Peskin and
			Schroeder in Ref. 60 (section 7.5). After the employment
			of regularization, followed by renormalization (discussed
			in Sec. A 6 d), one may extract the physical contribu-
			tion out of the divergent ∆E α,2 VP,1 . In the case where the
			Hamiltonian is invariant under time-reversal symmetry,

  F rngeleD wF holgD nd F hwerdtfegerD he tournl of ghemil hysis 138D IURIIQ @PHIQAF 40 F rngele nd wF holgD ghemil hysis vetters 616D PPP @PHIRAF 41 F rngele nd wF holgD he tournl of hemil physis 138D HRRIHR @PHIQAF 42 F hwerdtfegerD ghemhysghem 12D QIRQ @PHIIAF 43 wF holg nd F goD ghemF evF 112D RHQ @PHIPAF 44 eF xF ertemyevD i'etive ih rmiltoninsD in Handbook of Relativistic Quantum ChemistryD edited y F viu @pringer ferlin reidelergD PHITA pF I!IWF 45 F hwerdtfeger nd tF qF nijdersD tournl of hysis fX etomiD woleulr nd yptil hysis 23D QPPS!QPRH @IWWHAF 46 F uellö nd eF tF dlejD snterntionl journl of quntum hemistry 68D ISW @IWWVAF 47 uF qF hyllD snterntionl tournl of untum ghemistry 78D RIP!RPI @PHHHAF 48 F ueD F fstD eF F F qomesD rF tF eF tensenD vF issherD sF eF eurD F hi emigioD uF qF hyllD iF ilivD iF psshuerD et al.D he tournl of ghemil hysis 152D PHRIHR @PHPHAF 49 F ue nd F relgkerD tF gomputF ghemF 23D VIR @PHHPAF 50 F yykköD hvidson vetureD delivered for niversity of xorth exs for PS ytoer PHIQF eville tX http://www.chem. helsinki.fi/~pyykko/Videos/UNT.mp4F eugx is disussed t RUXRR of the videoF 51 eF vF gyrique qenetD prneso sntrvi nd F eynudD ennF pondF vF de froglie 29D QQI @PHHRAD elso ville from the pulisher9s wesiteX https://fondationlouisdebroglie.org/ AFLB-Web/fldb-annales-index.htmF 52 F impson nd F veonhrdtD edsFD Forces of the Quantum Vacuum An Introduction to Casimir Physics @orld ienE ti(D xew terseyD PHISAF 53 uF qF hyllD gF F fushliherD hF F hwenkeD nd F yykköD ghemF hysF vettF 348D RWU @PHHIAF 54 vF F kripnikovD he tournl of ghemil hysis 154D PHIIHI @PHPIAF 55 vF F kripnikovD hF F ghuukovD nd F wF hkhovD he tournl of ghemil hysis 155D IRRIHQ @PHPIAF 56 eF itsevskiiD vF F kripnikovD xF F wosyginD F ssevD F fergerD eF eF freierD nd F qiesenD tF ghemF hysF 224303

	taking it into consideration) should be the bare ones, in-	28 uF hyllD sF qrntD gF tohnsonD pF rpiD nd iF lummerD
	stead of the measured ones. This awareness played a crucial role in formulating the renormalization theory. Since we do not have access to bare quantities, and since	gomputer hysis gommunitions 55D RPS @IWVWAF 29 F hevD sF upitsynD nd F erokhinD gomputF hysF gomE munF 189D IUS @PHISAF 30 F hevD sF upitsynD nd F erokhinD gomputF hysF gomE
	innity is not natural (not measurable), the renormal-ization theory says that we are allowed to redene our physical constants such that the bare ones absorb the emerging divergences and lead to overall values of the physical constants that correspond to the experimentally	munF 223D TW @PHIVAF 31 iF F uhl nd tF gF ferengutD gomputer hysis gommuniE tions 238D PQP @PHIWAF 32 vF pF ²tekD iF ilivD eF forshevskyD F uldorD nd F hwerdtfegerD hysil eview vetters 118D HPQHHP @PHIUAF 33 fF gF heplerD xF fF flnovD nd uF eF etersonD he tourE
	observed ones. Detailed discussions on renormalization in the quantum eld theory are provided by Collins in Ref. 182, Greiner and Reinhardt in Ref. 85 (chapter 5), Peskin and Schroeder in Ref. 60 (chapter 7), Itzykson and Zuber in Ref. 86 (section 7.1), and Huang in Ref. 183	nl of hysil ghemistry e 109D IHQTQ @PHHSAF 34 fF gF heplerD xF fF flnovD nd uF eF etersonD he tourE nl of hemil physis 127D ITRQHR @PHHUAF 35 uF eF etersonD he tournl of hemil physis 142D HURIHS @PHISAF 36 F wF goxD wF gitirD F fF ermentroutD F F ftteyD nd uF eF
	(chapter 13) and 184. 1 uF qF hyll nd uF paegri trD Introduction to relativistic quan-tum chemistry @yxford niversity ressD PHHUAF 2 wF eiher nd eF olfD Relativistic quantum chemistry: the fundamental theory of molecular science @tohn iley 8 onsD	etersonD tF ghemF hysF 144D IVRQHW @PHITAF 37 F peng nd uF eF etersonD he tournl of ghemil hysis 147D HVRIHV @PHIUAF 38 F rngeleD wF holgD wF rnrthD F goD nd F hwerdtE fegerD he tournl of hemil physis 136D PIRIHS @PHIPAF
	PHIRAF
	3 F ueD ghemhysghem 12D QHUU @PHIIAF 4 F yykköD ennul review of physil hemistry 63D RS @PHIPAF 5 F hwerdtfegerD vF pF ²tekD eF unnettD nd F yF fowE
	mnD xuF hysF e 944D SSI @PHISAF 6 F yykko nd tF F hesluxD eounts of ghemil eserh
	12D PUT @IWUWAF 7 F yykkoD ghemil eviews 88D STQ @IWVVAF 8 F ehujD eF flomqvistD F vrssonD F yykköD nd F leskiE ijgierdD hysil review letters 106D HIVQHI @PHIIAF 9 iF ilivD F uldorD nd F sshikwD ghemil physis letters	.	. . .
	222D VP @IWWRAF 10 xF xF hutt nd F wjumderD hysil eview e 85D HQPSIP
	@PHIPAF
	11 eF heeD F ueD vF issherD nd eF evero ereir qomesD
	he tournl of hemil physis 149D IURIIQ @PHIVAF 12 F ue nd vF issherD in Theoretical Chemistry and Physics
	of Heavy and Superheavy Elements @pringerD PHHQA ppF PII!
	PTUF
	13 F viu nd sF vindgrenD he tournl of ghemil hysis 139D
	HIRIHV @PHIQAF
	14 F iF vm nd F gF etherfordD hysF evF 72D PRI @IWRUAF 15 F töhlkerD F woklerD pF foshD F hunfordD pF prnzkeD
	yF ulepperD gF uozhuhrovD F vudziejewskiD pF xoldenD rF eE
	ihD et al.D hysil review letters 85D QIHW @PHHHAF 16 F yykkoD ghemil reviews 112D QUI @PHIPAF 17 F yykköD wF okmnD nd vF vzowskyD hysil eview e
	57D TVW @IWWVAF 18 rF erssonD sF vindgrenD nd F lomonsonD hysi ript
	1993D IPS @IWWQAF 19 rF erssonD sF vindgrenD F lomonsonD nd F unnergrenD
	hysil eview e 48D PUUP @IWWQAF 20 wF uhlski nd uF hukiD hysil review letters 111D
	PRQHHI @PHIQAF
	21 wF ridrD FEF hongD F uoroovD nd tFEF urrD hysil
	eview e 101D HPPSHI @PHPHAF 22 F eoymD wF rykwD F uinoshitD nd wF xioD hysil
	eview vetters 109D IIIVHU @PHIPAF 23 iF eF ehlingD hysF evF 48D SS @IWQSAF 24 F yykkö nd vFEfF hoD tournl of hysis fX etomiD woleE
	ulr nd yptil hysis 36D IRTW @PHHQAF 25 F F plmum nd tF F wF qingesD hysF evF e 72D HSPIIS
	@PHHSAF
	26 F wF hevD sF sF upitsynD nd F eF erokhinD hysF evF
	e 88D HIPSIQ @PHIQAF 27 eF F wlyshevD hF eF qlzovD F wF hevD sF sF upitsynD
	F eF erokhinD nd F eF ytsevD hysF evF e 106D HIPVHT
	@PHPPAF

39

  gomputer hysis gommunitions 21D PQQ @IWVHAF 94 gF hierfelderD privte ommunitionF 95 eF hF fekeD he tournl of ghemil hysis 88D PSRU @IWVVAF 96 F sF veedev nd hF xF vikovD in Dokl. Math.D olF SW @wesu xukGsnterperiodiD IWWWA ppF RUU! RVID engulr qudrture prmeters ville fromX http://server.ccl.net/cca/software/SOURCES/FORTRAN/ Lebedev-Laikov-Grids/index.shtml. 97 F vindhD F ÅF wlmqvistD nd vF qglirdiD heoretil ghemistry eounts 106D IUV @PHHIAF 98 hsegD reltivisti initio eletroni struture progrmD eE lese hsegPP @PHPPAD written y rF tF eF tensenD F fstD eF F F qomesD F ue nd vF issherD with ontriutions from sF eF eurD F fkkenD gF ghiuezeD tF greutzergD uF qF hyllD F huillrdD F ikströmD iF ilivD F inevoldsenD iF pÿhuerD F pleigD yF possgrdD vF rlertD iF hF redegårdD F relgkerD fF relmih!risD tF renrikssonD wF vn rornD wF sli²D ghF F toD F unehtD F uomorovskýD yF uulE lieD tF uF vaerdhlD gF F vrsenD F F veeD xF rF vistD rF F xtrjD wF uF xykD F xormnD qF ylejnizkD tF ylsenD tF wF rF ylsenD eF pdopoulosD F gF rkD tF uF edE ersenD wF ernpointnerD tF F ototshnigD F di emigioD wF episkyD uF uudD F ªekD fF himmelpfennigD fF enjenD eF heeD tF ikkemD eF ungD eF tF horvldsenD tF hyssenD tF vn trlenD wF vF idlD F illumeD yF isserD F intherD F mmoto nd F un @ville t http://dx.doi.org/10. 5281/zenodo.6010450D see lso http://www.diracprogram.orgAF 99 eF ungD wF lmnD nd F ueD ixyhy @PHPPAD IHFSPVIGzenodoFTVURUPVF 100 vF issher nd uF qF hyllD etomi ht nd xuler ht les 67D PHU @IWWUAF

ACKNOWLEDGMENTS

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 101019907 HAMP-vQED), as well as from the Agence Nationale de la Recherche (ANR-17-CE29-0004-01 molQED). AS acknowledges nancial support from the Japan Society for the Promotion of Science (JSPS) KAKENHI Grant No. 17J02767, 20K22553, and 21K1464, and JSPS Overseas Challenge Program for Young Researchers, Grant No. 201880193. Computing time from CALMIP (Calcul en Midi-Pyrenes), supercomputer of ACCMS (Kyoto University) and Research Institute for Information Technology, Kyushu University (General Projects) are gratefully acknowledged.

We would like to thank Kirk Peterson (Washington State), Toshiaki Okabayashi (Shizuoka), Pekka Pyykkö (Helsinki), Jacinda Ginges (Brisbane) and Radovan Bast (Tromsø), for valuable discussions.

AUTHOR'S CONTRIBUTIONS

The overall project was conceived and supervised by TS. All programming and calculations were carried out by AS, whereas MS is the main contributor to the theory sections, notably the appendix.

DATA AVAILABILITY

The data that support the ndings of this study are openly available in ZENODO at https://doi.org/10. 5281/zenodo.6874728, see also Ref. 99. The DIRAC outputs contain information about the precise build of the corresponding executable and, most importantly, the git commit hash, which provides a unique identier of the precise version of DIRAC generating the output, hence allowing reproduction of the results. We have also included the potential surface generated by Hill et al.

116

for AuCN, kindly made available by Kirk Peterson.

CONFLICTS OF INTEREST

The authors have no conicts to disclose.

Appendix A: Theory background Since we hope to reach a wider audience than QED specialists, we provide in this Appendix a compact, yet accessible introduction (crash course) to QED that would otherwise have necessitated consulting disparate sources.

More precisely, in this Appendix, we shall discuss the lowest-order BSQED corrections, and show how the eective potentials associated with these QED processes can be derived within the scattering matrix (S -matrix) formalism. These eective potentials are to be used in practical relativistic calculations in order to account for the physics that is missing from the Dirac theory. In conventional QED one studies how the free-electron eld interacts with the free quantized electromagnetic eld and/or with an potential source (the scattering problem). On the other hand, BSQED theory studies the same problem but with electrons that are already interacting with some time-independent external eld, i.e. their wavefunctions are solutions to the bound-state Dirac equation instead of the free one. This is known as the Furry picture of quantum mechanics; see for instance Refs. 65 and 59 (section 15g).

We shall use unbold symbols for four-quantities such as spacetime points (events):

x = (ct, x), here in contravariant coordinates, where x is the spatial position vector, and the contravariant metric tensor g µν = diag (+1, -1, -1, -1). The gamma matrices are dened

Photon eld operator

The photon eld operator is written as a sum over positive and negative plane-wave Fourier modes

where N k = ℏ/(2ϵ 0 ω k V ) is the normalization constant, the zeroth component four-wave vector is k 0 = |k| = ω k /c, ε µ (k, r) are the four polarization vectors, and a (k, r) (and a † (k, r)) is the annihilation (creation) operator that annihilates (creates) a photon with wave vector k and polarization r, respectively (see Refs. 78 (eqs.(5.16a-c)) and 153 (section 8.4)). The choice of k 0 is imposed by the fact that the photon eld operator must satisfy the Maxwell equation

The boson creation and annihilation operators do satisfy the following bosonic commutation relations :

Here, ζ r is a function dened by the following relation

Recall that Ĥ1 I (t) is the interaction-picture version of the Schrödinger-picture interaction-Hamiltonian Ĥ1 S of Eq. (A19).

We shall note that the interaction density is multiplied by a damping factor e -ϵ ℏ |t| , cf. Eq. ( A26), where ϵ is a small positive quantity that has energy dimensions.

This term is known as the adiabatic switch that al-lows the interpolation between the perturbed and unperturbed problems (t = 0, ±∞), and was rst introduced by Gell-Mann and Low in Ref. 154 (Appendix A) while extending the S-matrix formalism to cover the boundelectron problem (see also Ref. 158 (section 1.3)). The scattering matrix of Eq. ( A26) may be expanded in powers of the perturbation parameter λ as

This form of the Ŝ-matrix + : Ô (x 1 ) Ô (x 2 ) Ô (x 3 ) Ô (x 4 ) . . . : + . . .

Contracted operators are moved next to each other, noting that under normal-ordering (fermionic) bosonic operators can be permuted as if they (anti)commuted. A contraction is represented by a line that links two operators and is dened to be the vacuum expectation value of the time-ordered product Ô (x 1 ) Ô (x 2 ) ≡ 0 T Ô (x 1 ) Ô (x 2 ) 0 .

(A33)

Since our QED interaction-Hamiltonian density contains electron and photon operators, the time-ordered product in our S-matrix of Eq. (A29) will be expanded with two types of contractions: electronic and photonic. The contraction of two electron eld operators (of Eq. (A7))

components α and β is dened with respect to the last formula by

The next step is to transform these two real-space integral S-matrices into Fourier-space ones. We rst use the electron and photon propagators of Eqs. (A37) and (A39) and write electron eld operators and the external (classical) potential in their Fourier-integral forms

We nally note that variable dependence indicates in which space the corresponding physical quantity is: We use x variables for spacetime points (coordinate-space), and p and q variables for four-momentum points (in momentum-space). The rst S-matrix Ŝ(2) SE,0 (ϵ, λ) (the zero-potential bound-state self-energy) becomes Ŝ(2) SE,0 (0, 1) =

Here Σ (q) is the so-called self-energy matrix function (see, for instance, Mandl and Shaw in Ref. 78 (eq.(9.20)).

Notice that in the limit of large momentum p, the integrand behaves as ∝ 1 γ µ (q-p)µ 1 p 3 , which indicates a su- percial linear divergence. However, with further investigation, one can show that this divergence is reduced by one degree, as noted by Schweber in Ref. 59 (section 15a), and presented in Table XVIII. Following the same steps, one can show that the second scattering matrix, associated with the one-potential bound-state self-energy process Ŝ(2) SE,0 (0, 1) can be written as

4 A e µ (p 2 -p 1 ) : Ψ(p 2 )Λ µ (p 2 , p 1 ) Ψ(p 1 ) :