Four-component relativistic range-separated density-functional theory: Short-range exchange local-density approximation

Julien Paquier * and Julien Toulouse † Laboratoire de Chimie Théorique (LCT), Sorbonne Université and CNRS, F-75005 Paris, France (Dated: Octobre 19, 2018) We lay out the extension of range-separated density-functional theory to a four-component relativistic framework using a Dirac-Coulomb-Breit Hamiltonian in the no-pair approximation. This formalism combines a wavefunction method for the long-range part of the electron-electron interaction with a density(-current) functional for the short-range part of the interaction. We construct for this formalism a short-range exchange local-density approximation based on calculations on a relativistic homogeneous electron gas with a modified Coulomb-Breit electron-electron interaction. More specifically, we provide the relativistic short-range Coulomb and Breit exchange energies per particle of the relativistic homogeneous electron gas in the form of Padé approximants which are systematically improvable to arbitrary accuracy. These quantities, as well as the associated effective Coulomb-Breit exchange hole, show the important impact of relativity on short-range exchange effects for high densities.

I. INTRODUCTION

Range-separated density-functional theory (RS-DFT) (see, e.g., Refs. 1 and 2) is an alternative to Kohn-Sham densityfunctional theory (DFT) [3] for electronic-structure calculations of atoms, molecules, and solids. It permits to rigorously combine an explicit wave-function calculation for the longrange part of the electron-electron interaction with a compact density functional for the complement short-range part of the electron-electron interaction. RS-DFT leads to a faster basis convergence than standard wave-function methods [4] and can provide improvement over usual Kohn-Sham DFT approximations for the description of strong-correlation effects (see, e.g., Refs. 5 and 6), weak intermolecular interactions (see, e.g., Refs. [7][8][9], fractionally charged subsystems (see, e.g., Refs. 10 and 11), and electronic excitation energies (see, e.g., Refs. [12][13][14].

For the description of compounds with heavy elements, relativistic effects have to be incorporated into RS-DFT. A simple approach that has been used consists in using standard scalar-relativistic effective-core potentials in RS-DFT [15,16]. A more sophisticated approach was developed by Kullie and Saue [17] who extended RS-DFT to a four-component relativistic Dirac-Coulomb Hamiltonian, using second-order Møller-Plesset (MP2) perturbation theory for the long-range part of the calculation and usual non-relativistic short-range semi-local exchange-correlation density-functional approximations. To make this approach more rigorous and possibly more accurate, especially for core properties, relativistic short-range density-functional approximations should be used. A number of relativistic density-functional approximations have been proposed for four-component relativistic Kohn-Sham DFT (see, e.g., reviews in Refs. 18 and 19), but no relativistic short-range density-functional approximations has been developed yet. This is unfortunate since relativistic effects are most important in spatial regions of high density whose contribution to the energy in RS-DFT mainly comes from the short-range exchange-correlation density functional.

In the present work, we remedy to this limitation by constructing a relativistic short-range exchange local-density approximation (LDA) based on a relativistic homogeneous electron gas (RHEG) with a modified electron-electron interaction. The choice to focus on the exchange energy and not on the correlation energy is motivated by the fact that, at high densities where relativistic effects are important, exchange largely dominates correlation, at least in not very strongly correlated systems. The choice to target a LDA, as opposed to a generalized-gradient approximation (GGA), is motivated by the fact that LDA is the standard first-level approximation to consider in DFT. Moreover, even though corrections to the relativistic LDA were shown to be important for the case of the full-range electron-electron interaction [START_REF] Engel | [END_REF]21], such corrections beyond LDA are expected to be much smaller for the case of a short-range electron-electron interaction, as known in non-relativistic RS-DFT [START_REF] Toulouse | [END_REF]22]. Beyond the goal of using this relativistic short-range exchange LDA in RS-DFT calculations of molecular and solid-state systems, the present work also aims at analyzing the importance of relativistic effects on the short-range exchange energy.

The paper is organized as follows.

In Section II, we briefly lay out the formalism of RS-DFT for a fourcomponent relativistic Dirac-Coulomb-Breit Hamiltonian. In Sections III A, III B, and III C, we review the calculation of the full-range exchange energy per particle of the RHEG with the standard Dirac-Coulomb-Breit Hamiltonian in a way that will prepare for the extension to the short-range case. We discuss the importance of the separate Coulomb and Breit contributions to the exchange energy per particle and to the exchange hole. In Section III D, we derive the exchange energy per particle of a RHEG with the short-range version of the Coulomb-Breit electron-electron interaction. Whereas the exchange energy per particle of a non-relativistic homogeneous electron gas with a short-range interaction can be obtained analytically quite easily [START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF]23,24], the calculation of the relativistic analogue turned out to be quite a formidable task. We did not manage to obtain a closed-form expression for the relativistic short-range exchange energy per particle, but we could ex-press it as a divergent series in inverse powers of the speed of light that can be summed to high accuracy using Padé approximants, which is suitable for all practical purposes. In Section III E, we derive expansions of the relativistic shortrange exchange energy per particle for small and large values of the range-separation parameter, and use them to construct a simple approximation to the relativistic short-range exchange energy per particle which can be used as an alternative to the Padé approximants. Finally, Section IV contains our conclusions and prospects for future work. Details of the calculations and the complete expressions obtained are given in the Supplementary Information [START_REF]See Supplementary Information for details on the calculation of the sums over spins in Eqs. (A1) and (A2), details on the calculation of the large-c expansions in Eqs. (36) and (41), and for a Mathematica notebook containing the complete expressions of the large-c expansions[END_REF].

II. RELATIVISTIC RANGE-SEPARATED DENSITY-FUNCTIONAL THEORY

A. No-pair Dirac-Coulomb-Breit wave-function theory

For relativistic electronic-structure calculations, a good starting point is the Dirac-Coulomb-Breit (DCB) electronic Hamiltonian in the no-pair approximation [START_REF] Sucher | [END_REF] ĤDCB

+ = ĤD + + ŴCB + , (1)
where ĤD + is the one-electron Dirac Hamiltonian and ŴCB + is the two-electron Coulomb-Breit (CB) interaction. In secondquantized position representation, ĤD + can be written as, in atomic units,

ĤD + = dr ψ † + (r) c (α • p) + β mc 2 + v ne (r)I 4 ψ+ (r), (2)
where p = -i∇ r is the momentum operator, c = 137.036 a.u. is the speed of light, and m = 1 a.u. is the electron mass that has been kept for clarity, v ne (r) is the electron-nuclei scalar potential energy operator, I 4 is the 4 × 4 identity matrix, and α and β are the 4 × 4 Dirac matrices

α = 0 2 σ σ 0 2 and β = I 2 0 2 0 2 -I 2 , (3)
where σ = (σ x , σ y , σ z ) is the 3-dimensional vector of the 2×2 Pauli matrices, and 0 2 and I 2 are the 2 × 2 zero and identity matrices, respectively. In Eq. ( 2), ψ+ (r) and ψ † + (r) are the (projected) annihilation and creation field operators

ψ+ (r) = + p ψ p (r) âp and ψ † + (r) = + p ψ † p (r) â † p , (4)
where the sum is over a set of orthonormal 4-componentspinor orbitals {ψ p (r)} which are positive-energy eigenfunctions of a Dirac Hamiltonian with some potential [27], and âp and â † p are the corresponding annihilation and creation operators of these orbitals. The restriction in the sum to positive-energy states corresponds to the no-pair approximation, in which negative-energy states corresponding to positronic states are disregarded. The two-electron Coulomb-Breit interaction term is written as

ŴCB + = 1 2 dr 1 dr 2 ψ † + (r 1 ) ψ † + (r 2 )
w CB (r 12 ) ψ+ (r 2 ) ψ+ (r 1 ), (5) where w CB (r 12 ) is the sum of the Coulomb and Breit contributions

w CB (r 12 ) = w C (r 12 ) + w B (r 12 ), (6)
where r 12 = r 1r 2 and r 12 = |r 12 |. The Coulomb interaction is w C (r 12 ) = w ee (r 12 ) (I 4 ) 1 (I 4 ) 2 ,

where w ee (r 12 ) = 1/r 12 , and (I 4 ) 1 and (I 4 ) 2 are the 4×4 identity matrices acting on electron 1 and 2, respectively. The Breit interaction is

w B (r 12 ) = - 1 2 w ee (r 12 )       α 1 • α 2 + (α 1 • r 12 ) (α 2 • r 12 ) r 2 12       ,( 8
)
where the Dirac matrices α 1 and α 2 act on electron 1 and 2, respectively. The Coulomb-Breit interaction corresponds to the single-photon exchange electron-electron scattering amplitude in quantum-electrodynamics (QED) evaluated with the zero-frequency limit of the photon propagator in the Coulomb electromagnetic gauge. More specifically, the instantaneous Coulomb interaction corresponds to the longitudinal component of the photon propagator, whereas the Breit interaction is obtained from to the zero-frequency transverse component of the photon propagator. The Breit interaction comprises the instantaneous magnetic Gaunt interaction, -w ee (r 12 )α 1 •α 2 , and the remaining lowest-order retardation correction (see, e.g., Ref. 28).

The no-pair Hamiltonian is not unique since it depends on the choice of the set of orbitals {ψ p }. It has been proposed [START_REF] Saue | Theoretical Chemistry and Physics ofHeavy and Superheavy Elements[END_REF][START_REF] Almoukhalalati | [END_REF] to define the no-pair relativistic ground-state energy of a N-electron system using a minmax principle, that we will formally write as

E = min Ψ max {ψ p } Ψ| ĤD + + ŴCB + |Ψ , (9)
where the maximization is done with respect to the set of positive-energy orbitals {ψ p } [on which the Hamiltonian is projected via Eq. ( 4)] by rotations with its complement set of negative-energy orbitals, and the minimization is done with respect to normalized multideterminant wave functions Ψ within the N-electron space generated by the set of positiveenergy orbitals {ψ p }. In practice, Eq. ( 9) can be realized with a multiconfiguration self-consistent-field (MCSCF) algorithm which targets a saddle point in the parameter space [START_REF] Almoukhalalati | [END_REF][31][32].

For one-electron Dirac Hamiltonians, this type of minmax principle appears well founded (see Refs. [33][34][35][36]. For the nopair Dirac-Coulomb-Breit Hamiltonian with a correlated wave function, to the best of our knowledge this minmax principle has not been rigorously examined mathematically but the results of Ref. 30 on two-electron systems suggests that it is indeed a reasonable definition of the no-pair ground-state energy.

B. No-pair Dirac-Coulomb-Breit range-separated density-functional theory

The Dirac-Coulomb RS-DFT introduced by Kullie and Saue [17] can be readily extended to a Dirac-Coulomb-Breit Hamiltonian. We will do so using the general formalism of relativistic current-density-functional theory [37][38][39], even though we do not consider any external vector potential.

The starting point is the following decomposition of the Coulomb-form electron-electron potential w ee (r 12 ) = 1/r 12 appearing in Eqs. (7) and (8) w ee (r 12 ) = w lr,µ ee (r 12 ) + w sr,µ ee (r 12 ), (10) where w lr,µ ee (r 12 ) = erf(µr 12 )/r 12 and w sr,µ ee (r 12 ) = erfc(µr 12 )/r 12 are long-range and short-range potentials, respectively. Here, erf(x) = (2/ √ π)

x 0 e -t 2 dt is the error function, erfc(x) = 1erf(x) is the complementary error function, and µ is the range-separation parameter acting like an inverse smooth cutoff radius. We then assume that the no-pair relativistic groundstate energy can be expressed as

E = min Ψ max {ψ p } Ψ| ĤD + + ŴCB,lr,µ + |Ψ + ĒCB,sr,µ Hxc [n Ψ , j Ψ ] , (11)
where, as before, the maximization is over the set of projecting positive-energy orbitals {ψ p } by rotations with its complement set of negative-energy orbitals and the minimization is over normalized multideterminant wave functions Ψ within the N-electron space generated by the set of positiveenergy orbitals. In Eq. ( 11), ŴCB,lr,µ + is the long-range version of the two-electron Coulomb-Breit interaction as defined by Eqs. ( 5)-( 8) but with the substitution w ee (r 12 ) -→ w lr,µ ee (r 12 ), and ĒCB,sr,µ Hxc [n Ψ , j Ψ ] is a complement short-range Coulomb-Breit Hartree-exchange-correlation functional of the electron density n Ψ (r) = Ψ|n(r)|Ψ and current electron density j Ψ (r) = Ψ| ĵ(r)|Ψ , where n(r) = ψ † + (r) ψ+ (r) and ĵ(r) = ψ † + (r)cα ψ+ (r) are the density and current density operators, respectively. The term "complement short-range" means that the functional contains not only a pure short-range contribution but also a mixed long-range/short-range contribution [40]. We note that, even though Eq. ( 11) seems a natural extension of non-relativistic RS-DFT [START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF][START_REF] Toulouse | [END_REF], the existence of a universal density-current functional ĒCB,sr,µ Hxc [n, j] giving the no-pair relativistic ground-state energy via the minmax procedure of Eq. ( 11) is not established. Indeed, similarly to the problem of defining static density functionals for excited-state energies (see, e.g., Refs. 41 and 42), it seems a priori only possible to define a no-pair relativistic functional which either is universal but satisfies only a stationary principle (instead of a minmax principle) or satisfies a minmax principle but is not universal (i.e., depending on an external potential via the projecting orbitals {ψ p }). It might be necessary to go at the QED level to rigorously formulate a RS-DFT approach, as it was done for Kohn-Sham DFT [START_REF] Engel | Relativistic Electronic Structure Theory, Part 1: Fundamentals[END_REF]. These aspects are beyond the scope of the present work.

The limiting cases of the relativistic RS-DFT approach of Eq. ( 11) should be mentioned. For µ → ∞, the long-range interaction reduces to the Coulomb form, w lr,µ→∞ ee (r 12 ) = 1/r 12 , and the complement short-range functional vanishes, ĒCB,sr,µ→∞ Hxc

[n] = 0, so Eq. ( 11) reduces to the wave-function theory of Eq. ( 9). For µ = 0, the long-range interaction vanishes, w lr,µ=0 ee (r 12 ) = 0, and the complement shortrange functional reduces to a full-range density functional, ĒCB,sr,µ=0

Hxc [n] = E CB Hxc [n]
, so it reduces to a no-pair relativistic Kohn-Sham DFT method (see, e.g., Ref. 43), and the corresponding minimizing wave function in Eq. ( 11) is a single Slater determinant Ψ µ=0 = Φ.

We decompose now the short-range density-current functional into Hartree, exchange, and correlation contributions ĒCB,sr,µ

Hxc [n, j] = E CB,sr,µ H [n, j] + E CB,sr,µ x [n, j] + ĒCB,sr,µ c [n, j]. (12)
The short-range Hartree contribution, containing Coulomb and Breit contributions, is an explicit functional of the density n and the current j

E CB,sr,µ H [n, j] = 1 2 w sr,µ ee (r 12 ) n(r 1 )n(r 2 ) dr 1 dr 2 - 1 4c 2 w sr,µ ee (r 12 ) j(r 1 ) • j(r 2 ) dr 1 dr 2 + w sr,µ ee (r 12 ) j(r 1 ) • r 12 j(r 2 ) • r 12 r 2 12 dr 1 dr 2 . ( 13
)
Note that, in relativistic Kohn-Sham DFT, the Hartree energy is usually defined with the full QED photon propagator in the Feynman gauge [START_REF] Engel | Relativistic Electronic Structure Theory, Part 1: Fundamentals[END_REF]. Here, instead, we define the short-range Hartree energy with the two-electron Coulomb-Breit interaction in the Coulomb gauge, in order to be consistent with the corresponding long-range wave-function contribution. The short-range exchange density-current functional is consequently defined by

E CB,sr,µ x [n, j] = Φ[n, j]| ŴCB,sr,µ + |Φ[n, j] -E CB,sr,µ H [n, j], (14)
where Φ[n, j] is the Kohn-Sham determinant associated with density n and current j, and ŴCB,sr,µ + is the short-range version of the two-electron Coulomb-Breit interaction obtained as in Eqs. ( 5)-( 8) but with the substitution w ee (r 12 ) -→ w sr,µ ee (r 12 ). In Eq. ( 12), ĒCB,sr,µ c [n, j] is the complement short-range correlation functional including all interaction effects beyond Hartree and exchange. Finally, note that, for closed-shell systems, the current vanishes, j = 0, and thus we have simply short-range functionals of the density only. We only consider this case in this work.

In practice, the application of the relativistic RS-DFT approach of Eq. ( 11) requires to use approximations for the long-range wave-function part and for the short-range exchange and correlation functionals. For the long-range wavefunction part, one can use standard approximations: Hartree-Fock, MP2, MCSCF, etc. For the short-range functionals, no approximation including the relativistic effects has been proposed so far. The rest of paper is devoted to the development of a relativistic LDA for the short-range exchange functional

E CB,sr,µ x,LDA [n] = n(r) CB,sr,µ x n(r) dr, (15)
where CB,sr,µ x (n) is the exchange energy per particle of a RHEG with the short-range two-electron Coulomb-Breit interaction. The calculation of the latter quantity is quite complicated and was performed with the help of the software Wolfram Mathematica [START_REF]Mathematica, Version[END_REF].

III. EXCHANGE EFFECTS IN THE RELATIVISTIC HOMOGENEOUS ELECTRON GAS

We consider a homogeneous electron gas, i.e. a box of volume V containing N electrons with electronic density n = N/V characteristic of the system and a uniform background of positive charges with density n b equal to n in order to ensure the electroneutrality of the system. The electron gas is studied in the thermodynamic limit, i.e. N → ∞ and V → ∞ while n = N/V is kept constant, the positive background cancelling the diverging Hartree energy term. Such an electron gas is considered to be relativistic when the Fermi wave vector k F = (3π 2 n) 1/3 is non negligible compared to the speed of light c. This must really be understood as a comparison between the energy related to the Fermi wave vector k F c and the rest energy mc 2 where we take these quantities in atomic units.

In order to get an idea of the maximal value of k F that one encounters in an heavy element, one may consider a 1s hydrogen-like orbital, ψ 1s (r) = Z 3 /π 1/2 e -Zr where Z is the atomic number, and calculate the density of this doubly occupied 1s orbital at the nucleus: n = 2 |ψ 1s (0)| 2 = 2 Z 3 /π corresponding to a Fermi wave vector of k F = (6π) 1/3 Z = 2.66 Z. The heaviest elements having atomic numbers Z of about 100, this corresponds to a maximal value of k F of nearly 300 a.u..

A. One-electron part

The form of the non-interacting one-electron Dirac equation for this homogeneous electron gas is

c(α • p) + βmc 2 ψ k,σ (r) = E k ψ k,σ (r), ( 16
)
where ψ k,σ (r) is a one-electron wave function with wave vector k and "spin" index σ =↑, ↓ associated with the positive-energy eigenvalue

E k = √ k 2 c 2 + m 2 c 4 . (17)
The wave functions ψ k,σ (r) are four-component spinors of the form

ψ k,σ (r) = ϕ k,σ (r) χ k,σ (r) , (18)
where each component is itself a two-component spinor. The large component ϕ k,σ (r) has a plane-wave form

ϕ k,σ (r) = 1 √ V E k + mc 2 2E k e -ik•r ϕ σ , (19)
and the small component is obtained from the large component by

χ k,σ (r) = c(σ • k) E k + mc 2 ϕ k,σ (r), (20)
and ϕ σ is a two-component spin part with ϕ ↑ = 1 0 and

ϕ ↓ = 0 1
. What we call "spin" here really refers to the index σ identifying the two components of the large-component spinor.

As in the non-relativistic case, we will calculate the exchange energy of the relativistic homogeneous electron gas using these non-interacting one-electron wave functions ψ k,σ (r). We note however that, in the relativistic case, there is in principle an ambiguity on which one-electron wave functions to use to calculate the exchange energy. Indeed, adding for example the Hartree-Fock exchange potential in Eq. ( 16) will a priori lead to different one-electron wave functions since translational symmetry only imposes the realspace plane-wave form but not the spinor components. In the non-relativistic case, this ambiguity does not appear since in this case the one-electron plane-wave wave functions are completely determined by translational symmetry, independently of the one-electron potential used. We leave further discussions of this interacting point for future work.

B. Full-range Coulomb-Breit exchange energy

As the preparation for the case of the short-range Coulomb-Breit interaction, we recall the form of the full-range Coulomb-Breit exchange energy. The Coulomb exchange energy per particle can be expressed as [39,[START_REF] Jancovici | [END_REF]

(see Ap- pendix A) C x = - 3 k F 4π 5 6 + 1 3 c2 + 2 3 √ 1 + c2 arcsinh 1 c - 1 3 1 + c2 2 ln 1 + 1 c2 - 1 2 √ 1 + c2 -c2 arcsinh 1 c 2 , ( 21
)
where c = c/k F . The Breit exchange energy per particle has a similar form 46 (see Appendix A)

B x = 3 k F 4π 1 -2 1 + c2 1 -c2 -2 ln (c) + ln 1 + c2 + 2 √ 1 + c2 -c2 arcsinh 1 c 2 . ( 22
)
The Breit exchange energy per particle is an approximation to the exchange energy per particle obtained with the transverse component of the full QED photon propagator [38,39,[START_REF] Jancovici | [END_REF]. The exchange energy per particle obtained with the full QED photon propagator has in fact a simpler expression than the Coulomb-Breit one, thanks to the cancellation of many terms between the Coulomb and transverse components,

QED x = - 3 k F 4π 1 - 3 2 √ 1 + c2 -c2 arcsinh 1 c 2 . ( 23
)
It is interesting to consider the non-relativistic and ultrarelativistic limits of the Coulomb and Breit exchange energies. In the non-relativistic (NR) limit c → ∞ (or equivalently the low-density limit), only the Coulomb contribution survives

C,NR x = - 3 k F 4π and B,NR x = 0, ( 24
)
which is consistent with the fact that the Breit interaction is a purely relativistic phenomenon. In the opposite ultrarelativistic (UR) limit c → 0 (or equivalently the high-density limit), we have

C,UR x = - (1 + ln4) k F 4π and B,UR x = 3 k F 4π , (25)
i.e. the Breit contribution becomes the opposite of the nonrelativistic Coulomb exchange energy, and is larger in absolute value than the Coulomb contribution, implying that the total exchange energy becomes positive.

The different exchange energies per particle are plotted in Figure 1 as functions of the Fermi wave vector k F , up to k F = 300 a.u. which is about the largest value that could be encountered in an heavy element. We observe that, compared to the non-relativistic energy energy, relativity always reduces the exchange energy in absolute value, and this effect increases with the density. The relativistic effects for the Coulomb exchange contribution remain relatively small, even at high densities. By contrast, the Breit exchange contribution has a much larger effect at sufficiently high densities. The Coulomb-Breit exchange energy per particle is a good approximation to the exchange energy per particle obtained with the full QED photon propagator for k F c ≈ 137 a.u., which is consistent with the fact that the Breit interaction constitutes only the leading term in the expansion of the QED photon propagator in 1/c.

Finally, we note that the exchange energy per particle corresponding to the Gaunt interaction can also be similarly obtained, but it starts to deviate from the transverse exchange energy obtained from the QED photon propagator for substantially smaller values of k F than the Breit exchange energy does. For this reason, we prefer in this work to use the Breit

For the case of the Coulomb interaction, the associated exchange hole was derived by Ellis [47] and MacDonald and Vosko [39] (see Appendix A)

n C x (r 12 ) = - 9 4 n 1 (k F r 12 ) 2 j 1 (k F r 12 ) 2 + (1 -λ)A λ (k F r 12 ) 2 + λB λ (k F r 12 ) 2 , ( 27
)
with λ = 1/(1 + c2 ) and

A λ (k F r 12 ) = ∞ ν=0 (2ν + 1)!! (2ν + 1) j ν+1 (k F r 12 ) λ k F r 12 ν , B λ (k F r 12 ) = ∞ ν=0 (2ν + 1)!! (2ν + 1) j ν+2 (k F r 12 ) λ k F r 12 ν , (28)
where j ν are the spherical Bessel functions. We have extended this result to the case of the Breit interaction. The associated exchange hole is (see Appendix A)

n B x (r 12 ) = - 9 2 n 1 (k F r 12 ) 2 -j 1 (k F r 12 ) 2 + (1 -λ)A λ (k F r 12 ) 2 . ( 29
)
To the best of our knowledge this expression had not been derived before. This Breit exchange hole should not be interpreted as a modification of the pair density, but simply as giving after integration the Breit exchange energy.

Like the non-relativistic exchange hole, the relativistic Coulomb exchange hole is normalized to -1

n C x (r 12 )dr 12 = -1, (30)
but the Breit exchange hole is not. When the density increases from the low-density limit to the high-density limit, its integral varies from 0 to 1

0 ≤ n B x (r 12 )dr 12 ≤ 1, (31)
which confirms that it should be considered as an effective exchange hole with only purpose to give the Breit exchange energy.

In the non-relativistic limit, corresponding to λ = 0, the Coulomb exchange hole reduces to the well-known nonrelativistic form and the Breit exchange hole vanishes In the opposite ultra-relativistic limit, corresponding to λ = 1, we have 2 , and

n C,UR x (r 12 ) = - 9 4 n 1 (k F r 12 ) 2 j 1 (k F r 12 ) 2 + B 1 (k F r 12 )
n B,UR x (r 12 ) = 9 2 n 1 (k F r 12 ) 2 j 1 (k F r 12 ) 2 , (33)
where we observe again the ultra-relativistic limit of the Breit term is the exact opposite of the non-relativistic limit of the Coulomb contribution.

In Figure 2, we plot the effective Coulomb-Breit exchange pair-distribution function g CB x (r 12 ) = n CB x (r 12 )/n for different values of the density (or equivalently k F ). For a low value of the density, k F = 1 a.u., g CB x (r 12 ) has the usual shape of the exchange hole and is indistinguishable from the non-relativistic calculation, with the on-top value being g CB x (0) ≈ -1/2. For the higher values of the density, k F = 137 and 274 u.a., we observe that the Coulomb-Breit exchange hole becomes shallower. In particular, for k F = 274 a.u., the Coulomb-Breit exchange hole is almost flat, with a very small minimum not located at r 12 = 0 anymore. These dramatic modifications of the exchange effects due to relativity underline the importance of considering relativistic effects when calculating short-range exchange energies in systems containing high-density regions.

D. Short-range Coulomb-Breit exchange energy per particle

Since we have previously calculated the Coulomb-Breit exchange hole, the most straightforward way of calculating the short-range Coulomb-Breit exchange energy per particle is to integrate the exchange hole with the short-range electronelectron interaction CB,sr,µ x = 1 2 n CB x (r 12 ) w sr,µ ee (r 12 ) dr 12 .

Unfortunately, we did not manage to calculate this integral analytically. We thus instead follow the same route as for the calculation of the full-range exchange energy (see Appendix A), i.e. integrating over space variables first. Similarly to Eq. (A1), we obtain the short-range Coulomb exchange energy per particle as a Fourier-space integral where each wave vector spans the volume of the Fermi sphere

V k F C,sr,µ x = - 1 2n(2π) 6 V k F wsr,µ ee (k 12 ) E k 1 E k 2 + (k 1 • k 2 )c 2 + c 4 E k 1 E k 2 dk 1 dk 2 = 3k F 4π 1 0 1 0 d k1 d k2 k1 k2 1 c2 + k2 1 c2 + k2 2 k1 k2 + e -k1 + k2 2 μ 2 -e -k1 -k2 2 μ 2 μ2 + 2c 2 + k2 1 + k2 2 + 2 c2 + k2 1 c2 + k2 2 4 c2 + k2 1 c2 + k2 2        Ei        - k1 + k2 2 μ 2        -Ei        - k1 -k2 2 μ 2        + ln ( k1 -k2 ) 2 -ln ( k1 + k2 ) 2        , (35)
where wsr,µ ee (k 12 ) = (4π/k 2 12 )(1e -k 2 12 /4µ 2 ) is the Fourier transform of the short-range interaction in terms of the relative wave vector

k 12 = |k 12 | with k 12 = k 1 -k 2 , Ei(x) = - ∞ -x e -t /
t dt is the exponential integral function, and we have introduced the scaled variables μ = µ/k F , k1 = k 1 /k F , and k2 = k 2 /k F . Unfortunately we were unable to straightforwardly calculate the double integral in Eq. (35). To circumvent this difficulty we first make an asymptotic expansion of the integrand for c → ∞ and then integrate term by term to obtain the asymptotic series

C,sr,µ x ∼ k F ∞ i=0 α 2i ( μ) c2i , (36)
where the coefficients α 2i ( μ) can be obtained analytically.

Their derivation being quite lengthy, it is discussed in detail in the Supplementary Information [START_REF]See Supplementary Information for details on the calculation of the sums over spins in Eqs. (A1) and (A2), details on the calculation of the large-c expansions in Eqs. (36) and (41), and for a Mathematica notebook containing the complete expressions of the large-c expansions[END_REF]. We give here the first three coefficients.

α 0 ( μ) = 1 4π -3 + 4 √ π erf 1 μ μ + 2 2e -1 μ2 -3 μ2 -2 e -1 μ2 -1 μ4 , α 2 ( μ) = 1 12π 1 -6 √ π erf 1 μ μ3 -6 e -1 μ2 -3 μ4 + 12 e -1 μ2 -1 μ6 , α 4 ( μ) = 1 240π -13 + 12 √ π erf 1 μ 8 + 45 μ2 μ3 + 12 4e -1 μ2 -45 μ4 + 24 13e -1 μ2 -40 μ6 -648 e -1 μ2 -1 μ8 . ( 37
)
However, the asymptotic series of Eq. ( 36) diverges for c < 1, i.e. for k F 137 a.u.. To avoid this divergence, we use a diagonal Padé approximant [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory[END_REF] of even order M C,sr,µ

x,Padé = k F M/2 i=0 A 2i ( μ)/c 2i M/2 i=0 B 2i ( μ)/c 2i , ( 38
)
with the choice B 0 ( μ) = 1 without loss of generality, and the other coefficients A 2i ( μ) and B 2i ( μ) are uniquely determined, for a given M, from the coefficients α 2i ( μ) so that the asymptotic expansion of the Padé approximant matches the asymp-totic expansion of Eq. ( 36) up to order 2M (see Appendix B).

For example, we give here the coefficients for the diagonal Padé approximant of order 2

A 0 ( μ) = α 0 ( μ) and A 2 ( μ) = α 2 ( μ) - α 0 ( μ)α 4 ( μ) α 2 ( μ) , B 2 ( μ) = - α 4 ( μ) α 2 ( μ) . ( 39
)
We proceed similarly for the short-range Breit exchange energy per particle. Adapting Eq. (A2) to the case of the shortrange interaction gives the Fourier-space integral B,sr,µ x

= 1 2n(2π) 6 V k F wsr,µ ee (k 12 ) c 2 E k 1 E k 2 E k 2 + c 2 E k 1 + c 2 k 2 1 + E k 1 + c 2 E k 2 + c 2 k 2 2 dk 1 dk 2 = 3k F 4π 1 0 1 0 d k1 d k2 k1 k2 c2 -c2 + k2 1 c2 + k2 2 c2 + k2 1 c2 + k2 2        Ei        - k1 + k2 2 μ 2        -Ei        - k1 -k2 2 μ 2        + ln ( k1 -k2 ) 2 -ln ( k1 + k2 ) 2        ,( 40
)
which, after asymptotically expanding the integrand for c → ∞ and integrating term by term, turns into the asymptotic se-

∼ k F ∞ i=0 β 2i ( μ) c2i , ( 41
)
where the coefficients β 2i ( μ) can be obtained analytically. The first three of them are

β 0 ( μ) = 0, β 2 ( μ) = 1 60π 25 -12 √ π 3 + 5 μ2 erf 1 μ μ + 18 5 -2e -1 μ2 μ2 + 6 15 -7e -1 μ2 μ4 -48 1 -e -1 μ2 μ6 , β 4 ( μ) = 1 280π -77 + 12 √ π 10 + 42 μ2 + 105 μ4 erf 1 μ μ -60 7 -2e -1 μ2 μ2 -12 140 -37e -1 μ2 μ4
-24 70 -37e

-1 μ2 μ6 + 792 1 -e -1 μ2 μ8 . ( 42
)
The fact that β 0 ( μ) = 0 corresponds to the fact that the Breit exchange energy vanishes in the non-relativistic limit. Again, the asymptotic series of Eq. ( 41) diverges for c < 1, so we construct a diagonal Padé approximant of even order M B,sr,µ

x,Padé = k F M/2 i=0 C 2i ( μ)/c 2i M/2 i=0 D 2i ( μ)/c 2i , ( 43
)
with D 0 ( μ) = 1 and the other coefficients C 2i ( μ) and D 2i ( μ) are uniquely determined, for a given M, from the coefficients β 2i ( μ). We give here the coefficients of the diagonal Padé approximant of order 2

C 0 ( μ) = β 0 ( μ) and C 2 ( μ) = β 2 ( μ) - β 0 ( μ)β 4 ( μ) β 2 ( μ) , D 2 ( μ) = - β 4 ( μ) β 2 ( μ) . ( 44
)
The complete expressions of the large-c expansions in Eqs. ( 36) and ( 41) and their associated Padé approximants in Eqs. (38) and ( 43) are explicitly given up to an arbitrary order in the Mathematica notebook available in the Supplementary Information [START_REF]See Supplementary Information for details on the calculation of the sums over spins in Eqs. (A1) and (A2), details on the calculation of the large-c expansions in Eqs. (36) and (41), and for a Mathematica notebook containing the complete expressions of the large-c expansions[END_REF].

We want to check the accuracy of these Padé approximants to the short-range Coulomb-Breit exchange energy per particle and at which order M we can truncate the expansion. For this, we may check for the most demanding case of µ = 0 corresponding to the full-range interaction and for which we know the exact exchange energy per particle [Eqs. ( 21) and ( 22)]. We thus plot in Figure 3 the exact full-range Coulomb-Breit exchange energy per particle and its Padé approximants of orders 2, 4, and 6, as well as their relative errors ∆ CB x = ( CB x,Padé -CB x )/ CB x , as functions of the Fermi wave vector k F . We note in passing that the Padé approximants can be calculated either for the Coulomb and Breit exchange energies per particle separately, or directly for the total Coulomb-Breit exchange energy per particle, as done here for the plot. It turns out that these two ways of proceeding give very similar Padé approximants for the Coulomb-Breit exchange energy per particle, for example differing by at most about 10 -4 % for the Padé approximants of order 6. As seen in Figure 3, the error of all Padé approximants naturally increases with the density. The error of the Padé approximant of order 2 increases most rapidly with the density, the relative error exceeding 5% as soon as k F 200 a.u.. The relative error in the Padé approximant of order 4 is less than 0.5% until about k F = 200 a.u. and increases to about 4% at k F = 300 a.u.. The Padé approximant of order 6 has a maximal relative error of 0.5% at k F = 300 a.u.. We thus conclude that the Padé approximant of order 6 provides an excellent approximation for all the density range in which we are interested. We have also explicitly checked that this accuracy of the Padé approximant of order 6 is preserved for the short-range Coulomb-Breit exchange energy per particle for a non-zero value of µ by comparing to results obtained by numerical integration in Eq. ( 34). Hence, numerically, the Padé approximants appear to rapidly converge to the correct limit as M → ∞.

In Figure 4 we plot the relativistic short-range Coulomb-Breit exchange energy per particle, using the Padé approximant of order 6, as a function of the scaled range-separation parameter µ/k F for different values of k F , and compare it to the non-relativistic short-range exchange energy per particle (whose expression can be found in Refs. 1 and 24). For k F = 1 a.u., the relativistic exchange energy per particle is indistinguishable from its non-relativistic counterpart, for all values of µ. For higher values of the Fermi wave vector, k F = 137 and 274 a.u., the relativistic short-range Coulomb-Breit exchange energy per particle becomes much smaller, in absolute value, than the non-relativistic short-range exchange energy per particle. Also, it appears that, for large k F , the relativistic shortrange exchange energy per particle goes to zero when increasing µ/k F significantly faster than does its non-relativistic analogue.

E. Small-and large-µ expansions and simple approximation for the short-range Coulomb-Breit exchange energy per particle

We have determined the short-range Coulomb and Breit exchange energies per particle C,sr,µ x and B,sr,µ x in the form of Padé approximants with quite complicated coefficients. We now derive simple expressions valid for small and large values of the range-separation parameter µ and use them to construct simple approximations for C,sr,µ x and B,sr,µ x .

In order to obtain the expansions of C,sr,µ x and B,sr,µ x for µ → 0, we start from the asymptotic series in Eqs. (36) and (41), expand them with respect to µ, and extract the µ and µ 2 terms. For the short-range Coulomb exchange energy per particle, we obtain directly expressions independent of c for the linear and quadratic terms in µ C,sr,µ x

= C x + µ √ π - 3 µ 2 2π k F + O µ 3 , (45)
and for the short-range Breit exchange energy per particle we obtain expansions in c for the linear and quadratic terms in µ

B,sr,µ x = B x - µ √ π 3 5c 2 - 3 7c 4 + 3 9c 6 - 3 11c 8 + O 1 c10 + 3 µ 2 2π k F 1 c2 - 1 c4 + 1 c6 - 1 c8 + O 1 c10 + O µ 3 .( 46
)
In this expression, the series in 1/c can be exactly summed, and we can write Eq. ( 46) as B,sr,µ x

= B x - µ √ π (1 -f (c)) + 3 µ 2 2π k F (1 -g(c)) + O µ 3 , (47)
where

f (c) = 3 c2 -3 c3 arctan 1 c , (48)
and

g(c) = 1 1 + 1 c2 . ( 49
)
We note that the linear term in µ in Eq. ( 47) can also been found by inserting into Eq. ( 40) the Fourier transform of the linear term in µ of the expansion of the short-range interaction. However, such an approach cannot be used for higherorder terms in µ because these terms lead to divergences when inserting into Eq. ( 40). The expansion of the short-range Coulomb-Breit exchange energy per particle for µ → 0 is thus CB,sr,µ x

= CB x + µ √ π f (c) - 3 µ 2 2π k F g(c) + O µ 3 . (50)
Therefore, near µ = 0, the linear and quadratic terms in µ of CB,sr,µ x depend on c. In the non-relativistic limit, we have f (c → ∞) = 1 and g(c → ∞) = 1, and thus we correctly recover the non-relativistic expansion µ/ √ π -3 µ 2 /(2πk F ) [START_REF] Toulouse | [END_REF]24]. In the ultra-relativistic limit, we have f (c → 0) = 0 and g(c → 0) = 0, and thus the linear and quadratic terms in µ vanish.

In order to obtain the expansions of C,sr,µ x and B,sr,µ x for µ → ∞, we start from the asymptotic expansion of the Fourier transform of the short-range interaction wsr,µ ee (k 12 ) =

π µ 2 + O 1 µ 4 . (51)
Inserting then Eq. (51) into Eq. ( 35) leads to the expansion of the short-range Coulomb exchange energy per particle C,sr,µ x

= - k 3 F 24πµ 2 (1 + h(c)) + O 1 µ 4 , (52)
with

h(c) = 9 4 (c 2 + c4 ) - 9 4 c4 arcsinh 1 c 2 √ 1 + c2 -c2 arcsinh 1 c . (53)
Similarly, inserting Eq. (51) into Eq. ( 40) leads to the expansion of the short-range Breit exchange energy per particle

B,sr,µ x = k 3 F 12πµ 2 (1 -h(c)) + O 1 µ 4 . ( 54
)
The asymptotic expansion of the short-range Coulomb-Breit exchange energy per particle for µ → ∞ is thus CB,sr,µ x

= - k 3 F 24πµ 2 (3h(c) -1) + O 1 µ 4 . (55)
In the non-relativistic limit, we have h(c → ∞) = 1 and we thus recover the leading term of the asymptotic expansion of the non-relativistic short-range exchange energy per particle, -k 3 F /(12πµ 2 ) [START_REF] Toulouse | [END_REF]24]. In the ultra-relativistic limit, we have h(c → 0) = 0 and thus the leading term is k 3 F /(24πµ 2 ). We can see on Figure 5 that these small-and large-µ expansions indeed reproduce well, at the scale of the plot, the short-range Coulomb-Breit exchange energy per particle for µ/k F 0.25 and µ/k F 1, respectively.

We propose now simple approximations for C,sr,µ x and B,sr,µ x as rational interpolations between these small-and large-µ expansions. The short-range Coulomb exchange energy per particle can be approximated by C,sr,µ x

≈ C x + d C µ 1 + a C µ + b C µ 2 + c C µ 3 , (56)
where the coefficients are determined to recover the µ and µ 2 terms of the small-µ expansion of Eq. ( 45) and the 1/µ 2 term and the vanishing 1/µ 3 term of the large-µ expansion of Eq. ( 52)

a C = 3 2π 1/2 k F + 24π 3/2 k 3 F (1 + h(c)) ( C x ) 2 , b C = - 24π k 3 F (1 + h(c)) C x , c C = - 24π k 3 F (1 + h(c)) d C , d C = 1 π 1/2 + a C C x . (57)
In a similar way, we propose to approximate the short-range Breit exchange energy per particle by B,sr,µ x

≈ B x + d B µ 1 + a B µ + b B µ 2 + c B µ 3 , (58)
where the coefficients are again determined by imposing the small-µ expansion of Eq. ( 47) and the large-µ expansion of Eq. ( 54).

a B = 3(1 -g(c)) 2π 1/2 k F (1 -f (c)) + 12π 3/2 k 3 F (1 -h(c))(1 -f (c)) ( B x ) 2 , b B = 12π k 3 F (1 -h(c)) B x , c B = 12π k 3 F (1 -h(c)) d B , d B = - (1 -f (c)) π 1/2 + a B B x . (59)
Using the same form of approximation on directly the short-range Coulomb-Breit exchange energy per particle can lead to poles in µ, thus it is preferable to separately construct approximations for the Coulomb and Breit contributions and to sum them to get an approximation of the shortrange Coulomb-Breit exchange energy per particle. This simple approximation is reported on Figure 5 where it can be seen that it is a quite good approximation to the accurate shortrange Coulomb-Breit exchange energy per particle given by the Padé approximant of order 6. The relative error of this simple approximation is less than 5% for all values of k F and µ. For higher accuracy, one could construct rational interpolations with more coefficients determined from higher-order terms in the small-and large-µ expansions. However, while the large-µ expansion could easily be obtained at an arbitrary order, higher-order terms in the small-µ expansion would be obtained as large-c asymptotic series that might not be always easy to sum into a closed form.

IV. CONCLUSIONS

In this work, we have considered the extension of RS-DFT to a four-component relativistic framework using a Dirac-Coulomb-Breit Hamiltonian, and we have constructed a shortrange LDA exchange density functional based on calculations on the RHEG with a modified electron-electron interaction. More specifically, we have provided the relativistic shortrange Coulomb and Breit exchange energies per particle of the RHEG in the form of Padé approximants [Eqs. (38) and ( 43)], constructed from large-c asymptotic expansions (but without involving expansions with respect to µ) and which are systematically improvable to arbitrary accuracy. These quantities, as well as the associated effective Coulomb-Breit exchange hole of the RHEG, show the important impact of relativity on short-range exchange effects for high densities. We have also provided simpler approximations for the relativistic short-range Coulomb and Breit exchange energies per particle of the RHEG in the form of rational interpolations [Eqs. ( 56) and ( 58)] constructed from the exact small-and large-µ expansions (but without involving expansions with respect to c) of the short-range exchange energies per particle, which can also be used when a limited accuracy is sufficient (relative error less than 5%).

Possible continuation of this work includes the construction of a relativistic short-range LDA correlation density functional (even though relativistic effects are expected to be much smaller for correlation than for exchange), the construction of relativistic short-range GGA functionals, the inclusion of the dependence on the current for open-shell systems, and the implementation of these functionals in a four-component relativistic RS-DFT program for tests on atomic and molecular systems. In particular, for compounds with heavy elements, including relativistic effects in short-range functionals should have significant effects on quantities sensitive to atomic cores, such as total energies, NMR and EPR parameters, or X-ray spectra.

sphere

V k F C x = - 1 2N V 2 (2π) 6 V k F V σ 1 ,σ 2 =↓,↑ w ee (r 12 ) ψ † k 1 ,σ 1 (r 1 )ψ k 2 ,σ 2 (r 1 )ψ † k 2 ,σ 2 (r 2 )ψ k 1 ,σ 1 (r 2 ) dr 1 dr 2 dk 1 dk 2 = - 1 2N(2π) 6 V k F V w ee (r 12 ) e -ik 12 •r 12 E k 1 E k 2 + (k 1 • k 2 )c 2 + c 4 E k 1 E k 2 dr 1 dr 2 dk 1 dk 2 = - 1 2n(2π) 6 V k F wee (k 12 ) E k 1 E k 2 + (k 1 • k 2 )c 2 + c 4 E k 1 E k 2 dk 1 dk 2 = - 3 k F 4π 5 6 + 1 3 c2 + 2 3 √ 1 + c2 arcsinh 1 c - 1 3 1 + c2 2 ln 1 + 1 c2 - 1 2 √ 1 + c2 -c2 arcsinh 1 c 2 , (A1)
where we have introduced

k 12 = k 1 -k 2 , k 12 = |k 12 |, wee (k 12 ) = 4π/k 2
12 which is the Fourier transform of the Coulomb interaction potential, and c = c/k F . The final expression in Eq. (A1) to the expression given in Refs. 39 and 45. The Breit exchange energy per particle of the RHEG can be calculated in a similar way

B x = 1 4N V 2 (2π) 6 V k F V σ 1 ,σ 2 =↓,↑ w ee (r 12 ) ψ † k 1 ,σ 1 (r 1 )α 1 ψ k 2 ,σ 2 (r 1 ) • ψ † k 2 ,σ 2 (r 2 )α 2 ψ k 1 ,σ 1 (r 2 ) + ψ † k 1 ,σ 1 (r 1 )(α 1 • r 12 )ψ k 2 ,σ 2 (r 1 ) ψ † k 2 ,σ 2 (r 2 )(α 2 • r 12 )ψ k 1 ,σ 1 (r 2 ) r 2 12 dr 1 dr 2 dk 1 dk 2 = 1 2N(2π) 6 V k F V w ee (r 12 ) e -ik 12 •r 12 c 2 E k 1 E k 2 E k 2 + c 2 E k 1 + c 2 k 2 1 + E k 1 + c 2 E k 2 + c 2 k 2 2 dr 1 dr 2 dk 1 dk 2 = 1 2n(2π) 6 V k F wee (k 12 ) c 2 E k 1 E k 2 E k 2 + c 2 E k 1 + c 2 k 2 1 + E k 1 + c 2 E k 2 + c 2 k 2 2 dk 1 dk 2 = 3 k F 4π 1 -2 1 + c2 1 -c2 -2 ln (c) + ln 1 + c2 + 2 √ 1 + c2 -c2 arcsinh 1 c 2 , (A2)
which corresponds to the expression given in Ref. 46. The sums over spins in Eqs. (A1) and (A2) are explicitly calculated in the Supplementary Information [START_REF]See Supplementary Information for details on the calculation of the sums over spins in Eqs. (A1) and (A2), details on the calculation of the large-c expansions in Eqs. (36) and (41), and for a Mathematica notebook containing the complete expressions of the large-c expansions[END_REF].

Effective Coulomb-Breit exchange hole

The Coulomb exchange hole can be written as an integral over the wave vectors

n C x (r 12 ) = - 1 (2π) 6 n V k F e -ik 12 •r 12 E k 1 E k 2 + (k 1 • k 2 )c 2 + c 4 E k 1 E k 2 dk 1 dk 2 = - 2 16π 4 n k F 0 k 2 j 0 (kr 12 ) 1 + c √ k 2 + c 2 dk 2 + k F 0 k 2 j 0 (kr 12 ) 1 - c √ k 2 + c 2 dk 2 + 2 k F 0 k 3 √ k 2 + c 2 j 1 (kr 12 )dk 2 ,(A3)
where the last expression is obtained after integration over the angle coordinates, and j ν are the spherical Bessel functions. After repeated integrations by parts, using for a general function f , (A6) This is the expression given in Refs. 39 and 47. Note that there are some typos in the expression given in Ref. 47. We now extend the previous derivation to the case of the Breit interaction. The associated exchange hole can be expressed as

n B x (r 12 ) = 1 (2π) 6 n V k F e -ik 12 •r 12 c 2 E k 1 E k 2 E k 2 + c 2 E k 1 + c 2 k 2 1 + E k 1 + c 2 E k 2 + c 2 k 2 2 dk 1 dk 2 = 1 2π 4 n k F 0 k 2 1 j 0 (k 1 r 12 ) 1 - c k 2 1 + c 2 dk 1 × k F 0 k 2 2 j 0 (k 2 r 12 ) 1 + c k 2 2 + c 2 dk 2 , (A7)
which, after using the same integration by parts as before, gives n B x (r 12 ) = -

9 2 n 1 (k F r 12 ) 2 -j 1 (k F r 12 ) 2 + (1 -λ)A λ (k F r 12 ) 2 . (A8)       α N+1 α N+2 . . . α N+M                 ,(B3)
and then the coefficients A 0 , A 1 , ..., A N are simply given by

A i = i j=0 α i-j B j . (B4)
A lot about convergence of Padé approximants is known for the special case where F(z) is a so-called Stieltjes function. As prescribed in Ref. [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory[END_REF], it can be checked whether F(z) Obviously, if in addition the two limits are equal, then the two Padé sequences converge to F(z) as M → ∞ for any real positive z.

We have checked the convergence of the Padé approximants in the case of the full-range Coulomb and Breit exchange energies per particle for which we have explicit exact expressions, Eqs. ( 21) and (22). For the Coulomb term, defining F C

x (z) = -C x (c) with z = 1/c 2 , numerical investigations suggest that F C x (z) satisfies properties (1) to ( 4) and is thus a Stieltjes function. Moreover, the Padé approximants P M M+1 (z) and P M M (z) appear numerically to converge to the unique limit F C

x (z) as M → ∞ for real positive z. For the Breit term, defining F B x (z) = B x (c) with z = 1/c 2 , numerical investigations suggest that F B x (z) satisfies properties (1) to (3) but not (4), and thus is not a Stieltjes function. Nevertheless, the Padé approximants P M M+1 (z) and P M M (z) still appear numerically to converge to the unique limit F B x (z) as M → ∞ for real positive z. This does not come as a surprise since it is known that Padé approximants often converge for functions which are not Stieltjes. We expect the same convergence properties for the Padé approximants in the case of the short-range Coulomb and Breit interactions.

FIG. 1 .

 1 FIG. 1. Non-relativistic and relativistic exchange energies per particle for the Coulomb, Coulomb-Breit, and QED photon propagator electron-electron interactions as functions of the Fermi wave vector k F .

  k F r 12 ) 2 (k F r 12 ) 2 and n B,NR x (r 12 ) = 0. (32)

FIG. 2 .

 2 FIG. 2. Effective Coulomb-Breit exchange pair-distribution function at different values of k F as a function of k F r 12 .

FIG. 3 .

 3 FIG. 3. Left plot: Exact full-range Coulomb-Breit exchange energy per particle and its Padé approximants of orders 2, 4, and 6 as functions of the Fermi wave vector k F . Right plot: Relative errors, ∆ CB x = ( CB x,Padé -CB x )/ CB x , of the Padé approximants compared to the exact full-range Coulomb-Breit exchange energy per particle.

FIG. 4 .

 4 FIG.4. Non-relativistic and relativistic Coulomb-Breit short-range exchange energies per particle as functions of µ/k F . The relativistic Coulomb-Breit short-range exchange energy is obtained from the Padé approximant of order 6.

FIG. 5 .

 5 FIG.5. Relativistic Coulomb-Breit short-range exchange energy per particle, its small-and large-µ expansions of Eqs. (50) and (55), and its simple approximation of Eqs. (56) and (58) as functions of µ/k F . The reference relativistic Coulomb-Breit short-range exchange energy per particle is obtained from the Padé approximant of order 6.

k F 0 kr 12 j

 012 ν+2 j ν (kr 12 ) f (k)dk = k ν+2 r 12 j ν+1 (kr 12 ) f (k) ν+1 (kr 12 ) f (k)dk, (A4)we obtain the expression of the Coulomb exchange hole asn C x (r 12 ) = -9 4 n 1 (k F r 12 ) 2 j 1 (k F r 12 ) 2 + (1 -λ)A λ (k F r 12 ) 2 + λB λ (k F r 12 ) 2 , (A5)where λ = 1/(1 + c2 ) andA λ (k F r 12 ) = ∞ ν=0 (2ν + 1)!! (2ν + 1) j ν+1 (k F r 12 ) λ k F r 12 ν , B λ (k F r 12 ) = ∞ ν=0 (2ν + 1)!! (2ν + 1) j ν+2 (k F r 12 ) λ k F r 12 ν .

  is a Stieltjes function by verifying the following four properties: (1) F(z) is analytic in the cut complex plane |arg z| < π; (2) lim z→∞ F(z) = C where C is a nonnegative real constant; (3) F(z) has an asymptotic series representation of the signalternating form ∞ i=0 a i (-z) i where a i ≥ 0; (4) -F(z) is Herglotz, i.e. sgn[Im(-F(z))] = sgn[Im z]. If these four properties are satisfied, then it can shown that, for any real positive z, the Padé sequences P M M+1 (z) and P M M (z) both converge as M → ∞ and sandwich F(z)

Appendix A: Coulomb-Breit exchange energy per particle and exchange hole of the relativistic homogeneous electron gas

In this appendix, we review the calculation of the Coulomb-Breit exchange energy per particle and the exchange hole of the RHEG.

Coulomb-Breit exchange energy per particle

The Coulomb exchange energy per particle of the RHEG can be straightforwardly calculated by summing over spins, integrating over each space coordinates in the volume V, and finally integrating over each wave vector in the volume of the Fermi

Appendix B: Padé approximants

We recall here how to calculate Padé approximants and check their convergence [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory[END_REF]. Consider a function F which is asymptotic to the following (divergent) power series, as z → 0,

Its Padé approximant of order (N, M) is

where B 0 = 1 without loss of generality, and the other N +