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

TL;DR

This paper extends range-separated density-functional theory to a four-component relativistic framework, combining wave-function and density functional methods for accurate short-range exchange calculations in high-density relativistic systems.

Contribution

It introduces a relativistic short-range exchange local-density approximation based on a relativistic homogeneous electron gas, incorporating relativistic effects into short-range exchange energies.

Findings

Relativistic effects significantly impact short-range exchange energies at high densities.

The developed approximation is systematically improvable using Padé approximants.

The formalism effectively combines wave-function and density functional approaches for relativistic systems.

Abstract

We lay out the extension of range-separated density-functional theory to a four-component relativistic frame-work using a Dirac-Coulomb-Breit Hamiltonian in the no-pair approximation. This formalism combines a wave-function method for the long-range part of the electron-electron interaction with a density(-current) functionalfor the short-range part of the interaction. We construct for this formalism a short-range exchange local-densityapproximation based on calculations on a relativistic homogeneous electron gas with a modified Coulomb-Breitelectron-electron interaction. More specifically, we provide the relativistic short-range Coulomb and Breit ex-change energies per particle of the relativistic homogeneous electron gas in the form of Pad e approximantswhich are systematically improvable to arbitrary accuracy. These quantities, as well as the associated effectiveCoulomb-Breit exchange hole, show the important impact of relativity on short-range exchange effects for highdensities.

Figures

Equations

$${\hat{H}}_{+}^{\text{DCB}}$$

$${\hat{H}}_{+}^{\text{D}}$$

$$\displaystyle\bm{\alpha}=\left(\begin{array}[]{cc}\mathbf{0}_{2}&\bm{\sigma}\\ \bm{\sigma}&\mathbf{0}_{2}\\ \end{array}\right)~{}\text{and}~{}~{}\bm{\beta}=\left(\begin{array}[]{cc}\mathbf{I}_{2}&\mathbf{0}_{2}\\ \mathbf{0}_{2}&-\mathbf{I}_{2}\\ \end{array}\right),$$

$${\hat{\psi }}_{+}(\mathbf{r})=\sum _p^{+}{\psi }_p(\mathbf{r}){\hat{a}}_p\text{and}{\hat{\psi }}_{+}^{†}(\mathbf{r})=\sum _p^{+}{\psi }_p^{†}(\mathbf{r}){\hat{a}}_p^{†},$$

$${\hat{W}}_{+}^{\text{CB}}$$

$${\mathbf{w}}^{\text{CB}}({\mathbf{r}}_{12})={\mathbf{w}}^{\text{C}}(r_{12})+{\mathbf{w}}^{\text{B}}({\mathbf{r}}_{12}),$$

$${\mathbf{w}}^{\text{C}}(r_{12})=w_{\text{ee}}(r_{12})({\mathbf{I}}_4{)}_1({\mathbf{I}}_4{)}_2,$$

$${\mathbf{w}}^{\text{B}}({\mathbf{r}}_{12})=-\frac{1}{2}{w}_{\text{ee}}(r_{12})\left({\mathbf{\alpha }}_1\cdot {\mathbf{\alpha }}_2+\frac{({\mathbf{\alpha }}_1\cdot {\mathbf{r}}_{12})({\mathbf{\alpha }}_2\cdot {\mathbf{r}}_{12})}{r_{12}^2}\right),$$

$$E$$

$$w_{\text{ee}}(r_{12})=w_{\text{ee}}^{\text{lr},\mu }(r_{12})+w_{\text{ee}}^{\text{sr},\mu }(r_{12}),$$

$$E$$

$${\bar{E}}_{\text{Hxc}}^{\text{CB,sr},\mu }[n,\mathbf{j}]$$

$$E_{\text{H}}^{\text{CB,sr},{\mu}}[n,\mathbf{j}]=\frac{1}{2}\iint w_{\text{ee}}^{\text{sr},{\mu}}(r_{12})~{}n(\mathbf{r}_{1})n(\mathbf{r}_{2})~{}\text{d}\mathbf{r}_{1}\text{d}\mathbf{r}_{2}~{}~{}~{}~{}~{}~{}~{}~{}\\ -\frac{1}{4c^{2}}\Bigg{[}\iint w_{\text{ee}}^{\text{sr},{\mu}}(r_{12})~{}{\mathbf{j}(\mathbf{r}_{1})\cdot\mathbf{j}(\mathbf{r}_{2})}\;\text{d}\mathbf{r}_{1}\text{d}\mathbf{r}_{2}~{}~{}~{}~{}~{}~{}\\ +\iint w_{\text{ee}}^{\text{sr},{\mu}}(r_{12})~{}\frac{\mathbf{j}(\mathbf{r}_{1})\cdot\mathbf{r}_{12}~{}~{}~{}\mathbf{j}(\mathbf{r}_{2})\cdot\mathbf{r}_{12}}{r_{12}^{2}}\text{d}\mathbf{r}_{1}\text{d}\mathbf{r}_{2}\Bigg{]}.$$

$$E_{\text{x}}^{\text{CB,sr},\mu }[n,\mathbf{j}]=⟨\Phi [n,\mathbf{j}]|{\hat{W}}_{+}^{\text{CB,sr},\mu }|\Phi [n,\mathbf{j}]⟩-E_{\text{H}}^{\text{CB,sr},\mu }[n,\mathbf{j}],$$

$$\displaystyle E_{\text{x,LDA}}^{\text{CB,sr},{\mu}}[n]=\int n(\mathbf{r})~{}{\epsilon}_{\text{x}}^{\text{CB,sr},{\mu}}\Big{(}n(\mathbf{r})\Big{)}~{}\text{d}\mathbf{r},$$

$$\left(c(\mathbf{\alpha }\cdot \mathbf{p})+\mathbf{\beta }m{c}^2\right){\psi }_{\mathbf{k},\sigma }(\mathbf{r})=E_k{\psi }_{\mathbf{k},\sigma }(\mathbf{r}),$$

$$E_k=\sqrt{k^2{c}^2+m^2{c}^4}.$$

$${\psi }_{\mathbf{k},\sigma }(\mathbf{r})$$

$${\phi }_{\mathbf{k},\sigma }(\mathbf{r})=\frac{1}{\sqrt{V}}\sqrt{\frac{E_k+m{c}^2}{2E_k}}{e}^{-i\mathbf{k}\cdot \mathbf{r}}{\phi }_{\sigma },$$

$${\chi }_{\mathbf{k},\sigma }(\mathbf{r})=\frac{c(\mathbf{\sigma }\cdot \mathbf{k})}{E_k+m{c}^2}{\phi }_{\mathbf{k},\sigma }(\mathbf{r}),$$

$${ϵ}_{\text{x}}^{\text{C}}$$

$${\epsilon}_{\text{x}}^{\text{B}}=\frac{3~{}k_{\text{F}}}{4{\pi}}\Bigg{(}1-2\Big{(}1+{\tilde{c}}^{2}\Big{)}\Bigg{(}1-{\tilde{c}}^{2}\Bigg{(}-2~{}\text{ln}\left({\tilde{c}}\right)+\text{ln}\left(1+{\tilde{c}}^{2}\right)\Bigg{)}\Bigg{)}+2\left(\sqrt{1+{\tilde{c}}^{2}}-{\tilde{c}}^{2}\text{arcsinh}\left(\frac{1}{\tilde{c}}\right)\right)^{2}\Bigg{)}.$$

$${ϵ}_{\text{x}}^{\text{QED}}$$

$${ϵ}_{\text{x}}^{\text{C,NR}}$$

$${ϵ}_{\text{x}}^{\text{C,UR}}=-\frac{(1+\text{ln}4)k_F}{4\pi }\text{and}{ϵ}_{\text{x}}^{\text{B,UR}}=\frac{3k_{\text{F}}}{4\pi },$$

$${ϵ}_{\text{x}}^{\text{CB}}$$

$$n_{\text{x}}^{\text{C}}(r_{12})=-\frac{9}{4}n\frac{1}{(k_{\text{F}}r_{12})^{2}}\Bigg{[}j_{1}(k_{\text{F}}r_{12})^{2}+(1-{\lambda})A_{\lambda}(k_{\text{F}}r_{12})^{2}\\ +{\lambda}B_{\lambda}(k_{\text{F}}r_{12})^{2}\Bigg{]},$$

$$A_{\lambda }(k_{\text{F}}{r}_{12})=\sum _{\nu =0}^{\infty }\frac{(2\nu +1)!!}{(2\nu +1)}{j}_{\nu +1}(k_{\text{F}}{r}_{12}){\left(\frac{\lambda }{k_{\text{F}}{r}_{12}}\right)}^{\nu },$$

$$B_{\lambda }(k_{\text{F}}{r}_{12})=\sum _{\nu =0}^{\infty }\frac{(2\nu +1)!!}{(2\nu +1)}{j}_{\nu +2}(k_{\text{F}}{r}_{12}){\left(\frac{\lambda }{k_{\text{F}}{r}_{12}}\right)}^{\nu },$$

$$n_{\text{x}}^{\text{B}}(r_{12})=-\frac{9}{2}n\frac{1}{(k_{\text{F}}r_{12})^{2}}\Bigg{[}-j_{1}(k_{\text{F}}r_{12})^{2}+(1-{\lambda})A_{\lambda}(k_{\text{F}}r_{12})^{2}\Bigg{]}.$$

$$\int n_{\text{x}}^{\text{C}}(r_{12})\text{d}{\mathbf{r}}_{12}=-1,$$

Full text

I Introduction

Range-separated density-functional theory (RS-DFT) (see, e.g., Refs. Savin, 1996; Toulouse et al., 2004a) is an alternative to Kohn-Sham density-functional theory (DFT) Kohn and Sham (1965) for electronic-structure calculations of atoms, molecules, and solids. It permits to rigorously combine an explicit wave-function calculation for the long-range 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 Franck et al. (2015) and can provide improvement over usual Kohn-Sham DFT approximations for the description of strong-correlation effects (see, e.g., Refs. Pollet et al., 2002; Fromager et al., 2007), weak intermolecular interactions (see, e.g., Refs. Ángyán et al., 2005; Goll et al., 2005; Toulouse et al., 2009), fractionally charged subsystems (see, e.g., Refs. Vydrov et al., 2007; Mussard and Toulouse, 2017), and electronic excitation energies (see, e.g., Refs. Tawada et al., 2004; Rebolini et al., 2013; Fromager et al., 2013).

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 Fromager et al. (2009); Alam and Fromager (2012). A more sophisticated approach was developed by Kullie and Saue Kullie and Saue (2012) 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. Engel, 2002; C. van Wüllen, 2010), 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 E. Engel, S. Keller, A. Facco Bonetti, H. Müller, and R. M. Dreizler (1995); E. Engel, S. Keller, and R. M. Dreizler (1996), 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 Toulouse et al. (2004a); Toulouse et al. (2005). 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 four-component relativistic Dirac-Coulomb-Breit Hamiltonian. In Sections III.1, III.2, and III.3, 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.4, 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 Savin (1996); Gill et al. (1996); Toulouse et al. (2004b), 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 express 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.5, we derive expansions of the relativistic short-range 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 Paq .

II Relativistic range-separated density-functional theory

II.1 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 Sucher (1980)

[TABLE]

where $\hat{H}^{\text{D}}_{+}$ is the one-electron Dirac Hamiltonian and $\hat{W}_{+}^{\text{CB}}$ is the two-electron Coulomb-Breit (CB) interaction. In second-quantized position representation, $\hat{H}^{\text{D}}_{+}$ can be written as, in atomic units,

[TABLE]

where $\mathbf{p}=-i\bm{\nabla}_{\mathbf{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_{\text{ne}}(\mathbf{r})$ is the electron-nuclei scalar potential energy operator, $\mathbf{I}_{4}$ is the $4\times 4$ identity matrix, and $\bm{\alpha}$ and $\bm{\beta}$ are the $4\times 4$ Dirac matrices

[TABLE]

where $\bm{\sigma}=(\bm{\sigma}_{x},\bm{\sigma}_{y},\bm{\sigma}_{z})$ is the 3-dimensional vector of the $2\times 2$ Pauli matrices, and $\mathbf{0}_{2}$ and $\mathbf{I}_{2}$ are the $2\times 2$ zero and identity matrices, respectively. In Eq. (2), $\hat{\psi}_{+}(\mathbf{r})$ and $\hat{\psi}_{+}^{\dagger}(\mathbf{r})$ are the (projected) annihilation and creation field operators

[TABLE]

where the sum is over a set of orthonormal 4-component-spinor orbitals $\{\psi_{p}(\mathbf{r})\}$ which are positive-energy eigenfunctions of a Dirac Hamiltonian with some potential Mittleman (1981), and $\hat{a}_{p}$ and $\hat{a}_{p}^{\dagger}$ 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

[TABLE]

where $\mathbf{w}^{\text{CB}}(\mathbf{r}_{12})$ is the sum of the Coulomb and Breit contributions

[TABLE]

where $\mathbf{r}_{12}=\mathbf{r}_{1}-\mathbf{r}_{2}$ and $r_{12}=|\mathbf{r}_{12}|$. The Coulomb interaction is

[TABLE]

where $w_{\text{ee}}(r_{12})=1/r_{12}$, and $(\mathbf{I}_{4})_{1}$ and $(\mathbf{I}_{4})_{2}$ are the $4\times 4$ identity matrices acting on electron 1 and 2, respectively. The Breit interaction is

[TABLE]

where the Dirac matrices $\bm{\alpha}_{1}$ and $\bm{\alpha}_{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_{\text{ee}}(r_{12})\bm{\alpha}_{1}\cdot\bm{\alpha}_{2}$, and the remaining lowest-order retardation correction (see, e.g., Ref. Gorceix and Indelicato, 1988).

The no-pair Hamiltonian is not unique since it depends on the choice of the set of orbitals $\{\psi_{p}\}$. It has been proposed Saue and Visscher (2003); A. Almoukhalalati, S. Knecht, H. J. Aa. Jensen, K. G. Dyall, and T. Saue (2016) to define the no-pair relativistic ground-state energy of a $N$-electron system using a minmax principle, that we will formally write as

[TABLE]

where the maximization is done with respect to the set of positive-energy orbitals $\{\psi_{p}\}$ [on which the Hamiltonian is projected via Eq. (8)] by rotations with its complement set of negative-energy orbitals, and the minimization is done with respect to normalized multideterminant wave functions $\Psi$ within the $N$-electron space generated by the set of positive-energy orbitals $\{\psi_{p}\}$. In practice, Eq. (13) can be realized with a multiconfiguration self-consistent-field (MCSCF) algorithm which targets a saddle point in the parameter space Jensen et al. (1996); Thyssen et al. (2008); A. Almoukhalalati, S. Knecht, H. J. Aa. Jensen, K. G. Dyall, and T. Saue (2016). For one-electron Dirac Hamiltonians, this type of minmax principle appears well founded (see Refs. Talman, 1989; Datta and Devaiah, 1988; Grisemer and Siedentop, 1999; Dolbeault et al., 2000). For the no-pair 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. A. Almoukhalalati, S. Knecht, H. J. Aa. Jensen, K. G. Dyall, and T. Saue, 2016 on two-electron systems suggests that it is indeed a reasonable definition of the no-pair ground-state energy.

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

The Dirac-Coulomb RS-DFT introduced by Kullie and Saue Kullie and Saue (2012) can be readily extended to a Dirac-Coulomb-Breit Hamiltonian. We will do so using the general formalism of relativistic current-density-functional theory Rajagopal and Callaway (1973); Rajagopal (1978); MacDonald and Vosko (1978), 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_{\text{ee}}(r_{12})=1/r_{12}$ appearing in Eqs. (11) and (12)

[TABLE]

where $w^{\text{lr},\mu}_{\text{ee}}(r_{12})=\operatorname{erf}(\mu r_{12})/r_{12}$ and $w^{\text{sr},\mu}_{\text{ee}}(r_{12})=\operatorname{erfc}(\mu r_{12})/r_{12}$ are long-range and short-range potentials, respectively. Here, $\operatorname{erf}(x)=(2/\sqrt{\pi})\int_{0}^{x}e^{-t^{2}}\text{d}t$ is the error function, $\operatorname{erfc}(x)=1-\operatorname{erf}(x)$ is the complementary error function, and $\mu$ is the range-separation parameter acting like an inverse smooth cut-off radius. We then assume that the no-pair relativistic ground-state energy can be expressed as

[TABLE]

where, as before, the maximization is over the set of projecting positive-energy orbitals $\{\psi_{p}\}$ by rotations with its complement set of negative-energy orbitals and the minimization is over normalized multideterminant wave functions $\Psi$ within the $N$-electron space generated by the set of positive-energy orbitals. In Eq. (15), $\hat{W}_{+}^{\text{CB,lr},{\mu}}$ is the long-range version of the two-electron Coulomb-Breit interaction as defined by Eqs. (9)-(12) but with the substitution $w_{\text{ee}}(r_{12})\longrightarrow w^{\text{lr},\mu}_{\text{ee}}(r_{12})$, and $\bar{E}_{\text{Hxc}}^{\text{CB,sr},{\mu}}[n_{\Psi},\mathbf{j}_{\Psi}]$ is a complement short-range Coulomb-Breit Hartree-exchange-correlation functional of the electron density $n_{\Psi}(\mathbf{r})=\langle\Psi|\hat{n}(\mathbf{r})|\Psi\rangle$ and current electron density $\mathbf{j}_{\Psi}(\mathbf{r})=\langle\Psi|\hat{\mathbf{j}}(\mathbf{r})|\Psi\rangle$, where $\hat{n}(\mathbf{r})=\hat{\psi}^{\dagger}_{+}(\mathbf{r})\hat{\psi}_{+}(\mathbf{r})$ and $\hat{\mathbf{j}}(\mathbf{r})=\hat{\psi}^{\dagger}_{+}(\mathbf{r})c\bm{\alpha}\hat{\psi}_{+}(\mathbf{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 Toulouse and Savin (2006). We note that, even though Eq. (15) seems a natural extension of non-relativistic RS-DFT Savin (1996); Toulouse et al. (2004a), the existence of a universal density-current functional $\bar{E}_{\text{Hxc}}^{\text{CB,sr},{\mu}}[n,\mathbf{j}]$ giving the no-pair relativistic ground-state energy via the minmax procedure of Eq. (15) is not established. Indeed, similarly to the problem of defining static density functionals for excited-state energies (see, e.g., Refs. Görling, 1999; Ayers and Levy, 2009), 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 $\{\psi_{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 Engel (2002). These aspects are beyond the scope of the present work.

The limiting cases of the relativistic RS-DFT approach of Eq. (15) should be mentioned. For $\mu\to\infty$, the long-range interaction reduces to the Coulomb form, $w^{\text{lr},\mu\to\infty}_{\text{ee}}(r_{12})=1/r_{12}$, and the complement short-range functional vanishes, $\bar{E}_{\text{Hxc}}^{\text{CB,sr},{\mu\to\infty}}[n]=0$, so Eq. (15) reduces to the wave-function theory of Eq. (13). For $\mu=0$, the long-range interaction vanishes, $w^{\text{lr},\mu=0}_{\text{ee}}(r_{12})=0$, and the complement short-range functional reduces to a full-range density functional, $\bar{E}_{\text{Hxc}}^{\text{CB,sr},{\mu=0}}[n]=E_{\text{Hxc}}^{\text{CB}}[n]$, so it reduces to a no-pair relativistic Kohn-Sham DFT method (see, e.g., Ref. Saue and Helgaker, 2002), and the corresponding minimizing wave function in Eq. (15) is a single Slater determinant $\Psi^{\mu=0}=\Phi$.

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

[TABLE]

The short-range Hartree contribution, containing Coulomb and Breit contributions, is an explicit functional of the density $n$ and the current $\mathbf{j}$

[TABLE]

Note that, in relativistic Kohn-Sham DFT, the Hartree energy is usually defined with the full QED photon propagator in the Feynman gauge Engel (2002). 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

[TABLE]

where ${\Phi[n,\mathbf{j}]}$ is the Kohn-Sham determinant associated with density $n$ and current $\mathbf{j}$, and $\hat{W}_{+}^{\text{CB,sr},{\mu}}$ is the short-range version of the two-electron Coulomb-Breit interaction obtained as in Eqs. (9)-(12) but with the substitution $w_{\text{ee}}(r_{12})\longrightarrow w^{\text{sr},\mu}_{\text{ee}}(r_{12})$. In Eq. (LABEL:EHxcsrdecomp), $\bar{E}_{\text{c}}^{\text{CB,sr},{\mu}}[n,\mathbf{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, $\mathbf{j}=\mathbf{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. (15) requires to use approximations for the long-range wave-function part and for the short-range exchange and correlation functionals. For the long-range wave-function 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

[TABLE]

where ${\epsilon}_{\text{x}}^{\text{CB,sr},{\mu}}(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 Mat .

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_{\text{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\rightarrow{\infty}$ and $V\rightarrow{\infty}$ 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_{\text{F}}=(3{\pi}^{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 ${\hbar}k_{\text{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_{\text{F}}$ that one encounters in an heavy element, one may consider a 1s hydrogen-like orbital, {\psi}_{1\text{s}}(r)=\big{(}Z^{3}/{\pi}\big{)}^{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~{}|{\psi}_{1\text{s}}(0)|^{2}=2~{}Z^{3}/{\pi}$ corresponding to a Fermi wave vector of $k_{\text{F}}=(6{\pi})^{1/3}~{}Z=2.66~{}Z$. The heaviest elements having atomic numbers $Z$ of about 100, this corresponds to a maximal value of $k_{\text{F}}$ of nearly 300 a.u..

III.1 One-electron part

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

[TABLE]

where ${\psi}_{\mathbf{k},{\sigma}}(\mathbf{r})$ is a one-electron wave function with wave vector $\mathbf{k}$ and “spin” index $\sigma=\uparrow,\downarrow$ associated with the positive-energy eigenvalue

[TABLE]

The wave functions ${\psi}_{\mathbf{k},{\sigma}}(\mathbf{r})$ are four-component spinors of the form

[TABLE]

where each component is itself a two-component spinor. The large component ${\varphi}_{\mathbf{k},{\sigma}}(\mathbf{r})$ has a plane-wave form

[TABLE]

and the small component is obtained from the large component by

[TABLE]

and ${\varphi}_{\sigma}$ is a two-component spin part with {\small{\varphi}_{\uparrow}=\Big{(}\begin{array}[]{c}1\\ 0\\ \end{array}\Big{)}} and {\small{\varphi}_{\downarrow}=\Big{(}\begin{array}[]{c}0\\ 1\\ \end{array}\Big{)}}. What we call “spin” here really refers to the index $\sigma$ 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 ${\psi}_{\mathbf{k},{\sigma}}(\mathbf{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. (20) will a priori lead to different one-electron wave functions since translational symmetry only imposes the real-space 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.

III.2 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 Jancovici (1962); MacDonald and Vosko (1978) (see Appendix A)

[TABLE]

where $\tilde{c}={c}/{k_{\text{F}}}$. The Breit exchange energy per particle has a similar form Ramana and Rajagopal, 1982 (see Appendix A)

[TABLE]

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 Jancovici (1962); Rajagopal (1978); MacDonald and Vosko (1978). 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,

[TABLE]

It is interesting to consider the non-relativistic and ultra-relativistic limits of the Coulomb and Breit exchange energies. In the non-relativistic (NR) limit ${\tilde{c}}\to\infty$ (or equivalently the low-density limit), only the Coulomb contribution survives

[TABLE]

which is consistent with the fact that the Breit interaction is a purely relativistic phenomenon. In the opposite ultra-relativistic (UR) limit ${\tilde{c}}\to 0$ (or equivalently the high-density limit), we have

[TABLE]

i.e. the Breit contribution becomes the opposite of the non-relativistic 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_{\text{F}}$, up to $k_{\text{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_{\text{F}}\lesssim c\approx 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_{\text{F}}$ than the Breit exchange energy does. For this reason, we prefer in this work to use the Breit interaction.

III.3 Effective Coulomb-Breit exchange hole

A convenient way of analyzing the exchange effects in terms of the interelectronic distance $r_{12}$ is to introduce an effective Coulomb-Breit exchange hole $n_{\text{x}}^{\text{CB}}(r_{12})$ such that

[TABLE]

For the case of the Coulomb interaction, the associated exchange hole was derived by Ellis Ellis (1977) and MacDonald and Vosko MacDonald and Vosko (1978) (see Appendix A)

[TABLE]

with $\lambda=1/(1+\tilde{c}^{2})$ and

[TABLE]

where $j_{\nu}$ 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)

[TABLE]

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$

[TABLE]

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 [math] to $1$

[TABLE]

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 ${\lambda}=0$, the Coulomb exchange hole reduces to the well-known non-relativistic form and the Breit exchange hole vanishes

[TABLE]

In the opposite ultra-relativistic limit, corresponding to ${\lambda}=1$, we have

[TABLE]

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_{\text{x}}^{\text{CB}}(r_{12})=n_{\text{x}}^{\text{CB}}(r_{12})/n$ for different values of the density (or equivalently $k_{\text{F}}$). For a low value of the density, $k_{\text{F}}=1$ a.u., $g_{\text{x}}^{\text{CB}}(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_{\text{x}}^{\text{CB}}(0)\approx-1/2$. For the higher values of the density, $k_{\text{F}}=137$ and $274$ u.a., we observe that the Coulomb-Breit exchange hole becomes shallower. In particular, for $k_{\text{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.

III.4 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 electron-electron interaction

[TABLE]

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. (66), 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_{\text{F}}}$

[TABLE]

where $\tilde{w}_{\text{ee}}^{\text{sr},{\mu}}(k_{12})=(4\pi/k_{12}^{2})(1-e^{-k_{12}^{2}/4{\mu}^{2}})$ is the Fourier transform of the short-range interaction in terms of the relative wave vector $k_{12}=|\mathbf{k}_{12}|$ with $\mathbf{k}_{12}=\mathbf{k}_{1}-\mathbf{k}_{2}$, $\text{Ei}(x)=-\int_{-x}^{\infty}e^{-t}/t~{}\text{d}t$ is the exponential integral function, and we have introduced the scaled variables $\tilde{\mu}={\mu}/k_{\text{F}}$, $\tilde{k}_{1}=k_{1}/k_{\text{F}}$, and $\tilde{k}_{2}=k_{2}/k_{\text{F}}$. Unfortunately we were unable to straightforwardly calculate the double integral in Eq. (41). To circumvent this difficulty we first make an asymptotic expansion of the integrand for $\tilde{c}\rightarrow{\infty}$ and then integrate term by term to obtain the asymptotic series

[TABLE]

where the coefficients ${\alpha}_{2i}(\tilde{\mu})$ can be obtained analytically. Their derivation being quite lengthy, it is discussed in detail in the Supplementary Information Paq . We give here the first three coefficients.

[TABLE]

However, the asymptotic series of Eq. (42) diverges for $\tilde{c}<1$, i.e. for $k_{\text{F}}\gtrsim 137$ a.u.. To avoid this divergence, we use a diagonal Padé approximant Bender and Orszag (1999) of even order $M$

[TABLE]

with the choice $B_{0}(\tilde{\mu})=1$ without loss of generality, and the other coefficients $A_{2i}(\tilde{\mu})$ and $B_{2i}(\tilde{\mu})$ are uniquely determined, for a given $M$, from the coefficients ${\alpha}_{2i}(\tilde{\mu})$ so that the asymptotic expansion of the Padé approximant matches the asymptotic expansion of Eq. (42) up to order $2M$ (see Appendix B). For example, we give here the coefficients for the diagonal Padé approximant of order 2

[TABLE]

We proceed similarly for the short-range Breit exchange energy per particle. Adapting Eq. (67) to the case of the short-range interaction gives the Fourier-space integral

[TABLE]

which, after asymptotically expanding the integrand for $\tilde{c}\rightarrow{\infty}$ and integrating term by term, turns into the asymptotic series

[TABLE]

where the coefficients ${\beta}_{2i}(\tilde{\mu})$ can be obtained analytically. The first three of them are

[TABLE]

The fact that ${\beta}_{0}(\tilde{{\mu}})=0$ corresponds to the fact that the Breit exchange energy vanishes in the non-relativistic limit. Again, the asymptotic series of Eq. (47) diverges for $\tilde{c}<1$, so we construct a diagonal Padé approximant of even order $M$

[TABLE]

with $D_{0}(\tilde{\mu})=1$ and the other coefficients $C_{2i}(\tilde{\mu})$ and $D_{2i}(\tilde{\mu})$ are uniquely determined, for a given $M$, from the coefficients ${\beta}_{2i}(\tilde{\mu})$. We give here the coefficients of the diagonal Padé approximant of order 2

[TABLE]

The complete expressions of the large-$\tilde{c}$ expansions in Eqs. (42) and (47) and their associated Padé approximants in Eqs. (44) and (49) are explicitly given up to an arbitrary order in the Mathematica notebook available in the Supplementary Information Paq .

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 $\mu=0$ corresponding to the full-range interaction and for which we know the exact exchange energy per particle [Eqs. (27) and (28)]. 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 $\Delta{\epsilon}_{\text{x}}^{\text{CB}}=\left|({\epsilon}_{\text{x},\text{Pad\'{e}}}^{\text{CB}}-{\epsilon}_{\text{x}}^{\text{CB}})/{\epsilon}_{\text{x}}^{\text{CB}}\right|$, as functions of the Fermi wave vector $k_{\text{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_{\text{F}}\gtrsim 200$ a.u.. The relative error in the Padé approximant of order 4 is less than 0.5% until about $k_{\text{F}}=200$ a.u. and increases to about 4% at $k_{\text{F}}=300$ a.u.. The Padé approximant of order 6 has a maximal relative error of 0.5% at $k_{\text{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 $\mu$ by comparing to results obtained by numerical integration in Eq. (40). Hence, numerically, the Padé approximants appear to rapidly converge to the correct limit as $M\to\infty$.

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 $\mu/k_{\text{F}}$ for different values of $k_{\text{F}}$, and compare it to the non-relativistic short-range exchange energy per particle (whose expression can be found in Refs. Savin, 1996; Toulouse et al., 2004b). For $k_{\text{F}}=1$ a.u., the relativistic exchange energy per particle is indistinguishable from its non-relativistic counterpart, for all values of $\mu$. For higher values of the Fermi wave vector, $k_{\text{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_{\text{F}}$, the relativistic short-range exchange energy per particle goes to zero when increasing ${\mu}/k_{\text{F}}$ significantly faster than does its non-relativistic analogue.

III.5 Small- and large-$\mu$ 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 ${\epsilon}_{\text{x}}^{\text{C,sr},{\mu}}$ and ${\epsilon}_{\text{x}}^{\text{B,sr},{\mu}}$ 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 ${\mu}$ and use them to construct simple approximations for ${\epsilon}_{\text{x}}^{\text{C,sr},{\mu}}$ and ${\epsilon}_{\text{x}}^{\text{B,sr},{\mu}}$.

In order to obtain the expansions of ${\epsilon}_{\text{x}}^{\text{C,sr},{\mu}}$ and ${\epsilon}_{\text{x}}^{\text{B,sr},{\mu}}$ for ${\mu}\to 0$, we start from the asymptotic series in Eqs. (42) and (47), expand them with respect to $\mu$, and extract the ${\mu}$ and ${\mu}^{2}$ terms. For the short-range Coulomb exchange energy per particle, we obtain directly expressions independent of ${\tilde{c}}$ for the linear and quadratic terms in ${\mu}$

[TABLE]

and for the short-range Breit exchange energy per particle we obtain expansions in ${\tilde{c}}$ for the linear and quadratic terms in ${\mu}$

[TABLE]

In this expression, the series in $1/\tilde{c}$ can be exactly summed, and we can write Eq. (52) as

[TABLE]

where

[TABLE]

and

[TABLE]

We note that the linear term in $\mu$ in Eq. (LABEL:Breitsmallmu) can also been found by inserting into Eq. (46) the Fourier transform of the linear term in $\mu$ of the expansion of the short-range interaction. However, such an approach cannot be used for higher-order terms in $\mu$ because these terms lead to divergences when inserting into Eq. (46). The expansion of the short-range Coulomb-Breit exchange energy per particle for ${\mu}\to 0$ is thus

[TABLE]

Therefore, near $\mu=0$, the linear and quadratic terms in $\mu$ of ${\epsilon}_{\text{x}}^{\text{CB,sr},\mu}$ depend on $\tilde{c}$. In the non-relativistic limit, we have $f(\tilde{c}\to\infty)=1$ and $g(\tilde{c}\to\infty)=1$, and thus we correctly recover the non-relativistic expansion $\mu/\sqrt{\pi}-3~{}{\mu}^{2}/(2{\pi}k_{\text{F}})$ Toulouse et al. (2004b, a). In the ultra-relativistic limit, we have $f(\tilde{c}\to 0)=0$ and $g(\tilde{c}\to 0)=0$, and thus the linear and quadratic terms in $\mu$ vanish.

In order to obtain the expansions of ${\epsilon}_{\text{x}}^{\text{C,sr},{\mu}}$ and ${\epsilon}_{\text{x}}^{\text{B,sr},{\mu}}$ for ${\mu}\to\infty$, we start from the asymptotic expansion of the Fourier transform of the short-range interaction

[TABLE]

Inserting then Eq. (57) into Eq. (41) leads to the expansion of the short-range Coulomb exchange energy per particle

[TABLE]

with

[TABLE]

Similarly, inserting Eq. (57) into Eq. (46) leads to the expansion of the short-range Breit exchange energy per particle

[TABLE]

The asymptotic expansion of the short-range Coulomb-Breit exchange energy per particle for ${\mu}\to\infty$ is thus

[TABLE]

In the non-relativistic limit, we have $h(\tilde{c}\to\infty)=1$ and we thus recover the leading term of the asymptotic expansion of the non-relativistic short-range exchange energy per particle, $-k_{\text{F}}^{3}/(12\pi\mu^{2})$ Toulouse et al. (2004b, a). In the ultra-relativistic limit, we have $h(\tilde{c}\to 0)=0$ and thus the leading term is $k_{\text{F}}^{3}/(24\pi\mu^{2})$.

We can see on Figure 5 that these small- and large-${\mu}$ expansions indeed reproduce well, at the scale of the plot, the short-range Coulomb-Breit exchange energy per particle for $\mu/k_{\text{F}}\lesssim 0.25$ and $\mu/k_{\text{F}}\gtrsim 1$, respectively.

We propose now simple approximations for ${\epsilon}_{\text{x}}^{\text{C,sr},{\mu}}$ and ${\epsilon}_{\text{x}}^{\text{B,sr},{\mu}}$ as rational interpolations between these small- and large-${\mu}$ expansions. The short-range Coulomb exchange energy per particle can be approximated by

[TABLE]

where the coefficients are determined to recover the $\mu$ and $\mu^{2}$ terms of the small-$\mu$ expansion of Eq. (51) and the $1/\mu^{2}$ term and the vanishing $1/\mu^{3}$ term of the large-$\mu$ expansion of Eq. (58)

[TABLE]

In a similar way, we propose to approximate the short-range Breit exchange energy per particle by

[TABLE]

where the coefficients are again determined by imposing the small-$\mu$ expansion of Eq. (LABEL:Breitsmallmu) and the large-$\mu$ expansion of Eq. (60).

[TABLE]

Using the same form of approximation on directly the short-range Coulomb-Breit exchange energy per particle can lead to poles in ${\mu}$, thus it is preferable to separately construct approximations for the Coulomb and Breit contributions and to sum them to get an approximation of the short-range 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 short-range 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_{\text{F}}$ and ${\mu}$. For higher accuracy, one could construct rational interpolations with more coefficients determined from higher-order terms in the small- and large-$\mu$ expansions. However, while the large-$\mu$ expansion could easily be obtained at an arbitrary order, higher-order terms in the small-$\mu$ 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 short-range LDA exchange density functional based on calculations on the RHEG with a modified electron-electron interaction. More specifically, we have provided the relativistic short-range Coulomb and Breit exchange energies per particle of the RHEG in the form of Padé approximants [Eqs. (44) and (49)], constructed from large-$c$ asymptotic expansions (but without involving expansions with respect to $\mu$) 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. (62) and (64)] constructed from the exact small- and large-$\mu$ 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.

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.

A.1 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 sphere $V_{k_{\text{F}}}$

[TABLE]

where we have introduced $\mathbf{k}_{12}=\mathbf{k}_{1}-\mathbf{k}_{2}$, $k_{12}=|\mathbf{k}_{12}|$, $\tilde{w}_{\text{ee}}(k_{12})=4{\pi}/k_{12}^{2}$ which is the Fourier transform of the Coulomb interaction potential, and $\tilde{c}={c}/{k_{\text{F}}}$. The final expression in Eq. (66) corresponds to the expression given in Refs. Jancovici, 1962; MacDonald and Vosko, 1978. The Breit exchange energy per particle of the RHEG can be calculated in a similar way

[TABLE]

which corresponds to the expression given in Ref. Ramana and Rajagopal, 1982. The sums over spins in Eqs. (66) and (67) are explicitly calculated in the Supplementary Information Paq .

A.2 Effective Coulomb-Breit exchange hole

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

[TABLE]

where the last expression is obtained after integration over the angle coordinates, and $j_{\nu}$ are the spherical Bessel functions. After repeated integrations by parts, using for a general function $f$,

[TABLE]

we obtain the expression of the Coulomb exchange hole as

[TABLE]

where $\lambda=1/(1+\tilde{c}^{2})$ and

[TABLE]

This is the expression given in Refs. Ellis, 1977; MacDonald and Vosko, 1978. Note that there are some typos in the expression given in Ref. Ellis, 1977. We now extend the previous derivation to the case of the Breit interaction. The associated exchange hole can be expressed as

[TABLE]

which, after using the same integration by parts as before, gives

[TABLE]

Appendix B Padé approximants

We recall here how to calculate Padé approximants and check their convergence Bender and Orszag (1999). Consider a function $F$ which is asymptotic to the following (divergent) power series, as $z\to 0$,

[TABLE]

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

[TABLE]

where $B_{0}=1$ without loss of generality, and the other $N+M+1$ coefficients $A_{i}$ and $B_{i}$ are determined so that the power expansion of Eq. (75) matches the power series of Eq. (74) up to order $N+M$. This gives the following matrix equation for the determining the coefficients $B_{1},B_{2},...,B_{M}$

[TABLE]

and then the coefficients $A_{0},A_{1},...,A_{N}$ are simply given by

[TABLE]

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. Bender and Orszag, 1999, it can be checked whether $F(z)$ is a Stieltjes function by verifying the following four properties: (1) $F(z)$ is analytic in the cut complex plane $|\text{arg}\;z|<\pi$; (2) $\lim_{z\to\infty}F(z)=C$ where $C$ is a nonnegative real constant; (3) $F(z)$ has an asymptotic series representation of the sign-alternating form $\sum_{i=0}^{\infty}{a}_{i}(-z)^{i}$ where $a_{i}\geq 0$; (4) $-F(z)$ is Herglotz, i.e. $\text{sgn}[\text{Im}(-F(z))]=\text{sgn}[\text{Im}\;z]$. If these four properties are satisfied, then it can shown that, for any real positive $z$, the Padé sequences $P_{M+1}^{M}(z)$ and $P_{M}^{M}(z)$ both converge as $M\to\infty$ and sandwich $F(z)$

[TABLE]

Obviously, if in addition the two limits are equal, then the two Padé sequences converge to $F(z)$ as $M\to\infty$ 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. (27) and (28). For the Coulomb term, defining $F_{\text{x}}^{\text{C}}(z)=-{\epsilon}_{\text{x}}^{\text{C}}(\tilde{c})$ with $z=1/\tilde{c}^{2}$, numerical investigations suggest that $F_{\text{x}}^{\text{C}}(z)$ satisfies properties (1) to (4) and is thus a Stieltjes function. Moreover, the Padé approximants $P_{M+1}^{M}(z)$ and $P_{M}^{M}(z)$ appear numerically to converge to the unique limit $F_{\text{x}}^{\text{C}}(z)$ as $M\to\infty$ for real positive $z$. For the Breit term, defining $F_{\text{x}}^{\text{B}}(z)={\epsilon}_{\text{x}}^{\text{B}}(\tilde{c})$ with $z=1/\tilde{c}^{2}$, numerical investigations suggest that $F_{\text{x}}^{\text{B}}(z)$ satisfies properties (1) to (3) but not (4), and thus is not a Stieltjes function. Nevertheless, the Padé approximants $P_{M+1}^{M}(z)$ and $P_{M}^{M}(z)$ still appear numerically to converge to the unique limit $F_{\text{x}}^{\text{B}}(z)$ as $M\to\infty$ 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.