Relativistic reduced density matrix functional theory.
M. Rodrı́guez-Mayorga,∗ K.J.H. Giesbertz, and L. Visscher
Department of Theoretical Chemistry, Vrije Universiteit Amsterdam,
De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands.
(Dated: April 25, 2022)
As a new approach to efficiently describe correlation effects in the relativistic quantum world we
propose to consider reduced density matrix functional theory, where the key quantity is the firstorder reduced density matrix (1-RDM). In this work, we first introduce the theoretical foundations
to extend the applicability of this theory to the relativistic domain. Then, using the so-called no-pair(np) approximation, we arrive at an approximate treatment of the relativistic effects by focusing
on electronic wavefunctions and neglecting explicit contributions from positrons. Within the np
approximation the theory becomes similar to the nonrelativistic case, with as unknown only thefunctional that describes the electron-electron interactions in terms of the 1-RDM. This requires
the construction of functional approximations, and we therefore also present the relativistic versionsof some common RDMFT approximations that are used in the nonrelativistic context and discuss
their properties.

CONTENTS

I. Introduction

2

II. ReRDMFT including electron-positron pair creation and annihilation processes.
A. The free particle Dirac equation and quantization of the Dirac field
B. The Hamiltonian operator including an external nonlocal potential, charge sectors in Fock space, andBogoliubov transformations
C. Vacuum polarization
D. The interacting particle Hamiltonian and the configuration-interaction vacuum
E. The 1-RDM and its natural orbital representation
F. Introducing ReRDMFT through the constrained-search formalism

3
3
4
6
7
8
9

III. The no-pair approximation and ReRDMFT
A. Vacuum polarization effects at fermion-fermion interacting case
B. The no-pair approximation including vacuum polarization ReRDFMT approach
C. The no-pair approximation ReRDFMT approachcf |Ψ+ ⟩ term as an explicit functional of the second-order reduced density matrix.
D. The ⟨Ψ+ |W
E. np-ReNOFT functional approximations by imposing Kramers’ symmetry
1. The Dirac–Hartree–Fock approximation
2. f (nI , nJ )-functional approximations
3. Approximations built imposing N -representability conditions
F. Properties of rel-GNOF/PNOFx functionals.

10
10
11
12

IV. Closing remarks

20

Acknowledgements

13
13
15
15
16
19

20

A. Commutation relation between the Hamiltonian and the opposite charge operator

21

B. Extending Gilbert’s theorem to the relativistic domain.
The no-pair vacuum polarization framework

22
22

C. The npvp energy functional23∗ marm3.14@gmail.com

2
D. The N -representability conditions within the np approximation.
The D-, Q-, and G- N -representability conditions

23
24

E. The relativistic Fixed-Phases functional.

27

References

29

I.

INTRODUCTION

When relativistic effects play a role in the quantum world, the Schrödinger equation must be replaced by the
Dirac equation. This change of paradigm is required to describe the electronic structure of heavy elements in the
periodic table [1–5] as electrons reach high velocities due to their strongly attractive nuclear potentials. When suchelements are present, often (near) degeneracies in the electronic energies (e.g. due to spin-orbit coupling effects) occur
and demand a good treatment of the so-called nondynamic (strong/static) correlation effects [6–14]. While methods
like complete active space including second-order perturbation theory (CASPT2) [15–18], DMRG [19, 20] and multireference coupled cluster theory [21–26] can be employed they all exhibit a high computational scaling with systemsize that limits their general applicability. An interesting efficient alternative is offered by reduced density matrix
functional theory (RDMFT) as it has an intrinsic low-order scaling with system size [27, 28].
RDMFT is emerging as a strong competitor to the widely used density functional theory (DFT) due to the possibility
to use fractional occupation numbers, which facilitates the study of electronic systems where the so-called nondynamic
correlation effects are enhanced[29–35]. Indeed, the popular working horse of physicists and chemists (i.e. the useof Kohn–Sham DFT approach with the usual density-functional approximations), in general, fails to account for
nondynamic correlation effects [36] (except for a few cases like the approximation based on the fractional spin localized
orbital scaling correction [37], some functional approximations developed using strictly correlated electrons [38],Becke’s B13 functional [39, 40], among others.). Actually, the capability of RDMFT to render nondynamic correlation
effects has lead to a recent burst of this theory into new domains like the study of superconductivity [41] and of Bose–
Einstein condensates [42]. In an attempt to extend it to the relativistic context, we validate its applicability by setting
the theoretical foundations of relativistic RDMFT (ReRDMFT).
Already in 2002 Ohsaku et al. introduced a local (relativistic) quantum electro-dynamical (QED) RDMFT [43]theory in which the nonlocality properties of the first-order reduced density matrix (1-RDM) were not exploitedand the so-called noninteracting kinetic energy needed to be evaluated with an auxiliary noninteracting system.
This complication may be avoided by using the full 1-RDM as all one-body interactions can then be evaluated as
explicitly known functionals of the 1-RDM. In our work, we do therefore exploit also the nonlocality of the 1-RDM
and introduce relativistic RDMFT (ReRDMFT). In this theory, we consider an external nonlocal potential (as it
was employed by Gilbert in the nonrelativistic context [44]) and define the energy as a functional of the 1-RDM. In
this way, the functional expression for all the one-body interactions is fully known in terms of this matrix and only
the energy functional for electron-electron (as well as the positron-positron and the electron-positron) interactions
remains unknown.
Recently, Toulouse proposed a relativistic density-functional theory based on a Fock-space effective QED Hamiltonian using the Coulomb or Coulomb–Breit two-particle interaction [45]. This theory, based on the works of Chaix etal. (see Refs. 46, 47), includes vacuum polarization effects through the creation of electron-positron pairs. Following
his work we propose a similar approach called npvp-ReRDMFT that is capable to account for the effects of vacuum
polarization within the no-pair vacuum-polarization approximation. Finally, neglecting the effect of vacuum polarization and assuming that a floating vacuum state is taken as reference each time spinor rotations are applied, we arrive
to our last formulation called np-ReRDMFT. The np-ReRDMFT corresponds to the usual application of the no-pairapproximation [48, 49].
This work is organized as follows: 1) we introduce the most general ReRDMFT approach, where creation andannihilation of electron-positron pairs is allowed. To that end, we first discuss the single-particle problem subject to anonlocal external potential; then, we present its many-particle generalization and introduce the Fock space. Next, weanalyze some properties of this Fock space and show how wavefunctions can be split according to the different charge
sectors. Then, we focus on the effects of changing the basis representation that leads to the vacuum polarization
effects. Finally, we present the correlated wavefunction (that includes creation and annihilation of electron-positronpairs processes) and use it in combination with the constrained search formalism to introduce ReRDMFT. 2) As an
approximation to the full ReRDMFT we include the so-called no-pair approximation at two levels. An initial approach
where vacuum polarization effects are taken into account (leading to npvp-ReRDMFT) and a second approach wherethese effects are neglected (np-ReRDMFT). It is within the np-ReRDMFT framework that we evoke the Kramers’
symmetry [50] and use it to adapt the nonrelativistic RDMFT functional approximations to the relativistic context.

3
Finally, we discuss some of the properties of the functional approximations. Before proceeding, let us present Table I
where we have collected the index conventions that will be employed in this work.
TABLE I. Indices used throughout this work.
Indices
I,J,K,L
R,S
A,B

labeling
positive energy spinors (PS)negative energy spinors (NS)all spinors (NS ∪ PS)
half of the PS
i,j,k,l(not related by Kramers’ symmetry)half of the PSī,j̄,k̄,¯l(Kramers partners of the unbarred spinors)µ, ν, τ , η scalar orbitals (components of a 4-spinor)p,qelectron pairs

II.

ReRDMFT INCLUDING ELECTRON-POSITRON PAIR CREATION AND ANNIHILATIONPROCESSES.
A.

The free particle Dirac equation and quantization of the Dirac fieldLet us start by defining the time-independent free particle Dirac equationb D (r)ψA (r) = −ic(αr · ∇r ) + c2 mβ ψA (r)T= EA ψA (r),(1)

b D (r) is the usual first-quantized 4 × 4 Dirac kinetic + rest mass operator, i =where T
speed of light, m = 1 a.u. is the electron mass,√−1, c = 137.036 a.u. is theα = (αx , αy , αz )
!
!
!!
02 σx
02 σy
02 σz
=
,
,
,
σx 02
σy 02
σz 02

(2)

02 is the 2×2 null matrix,
σx =

!
0 1
,
1 0

σy =

!
0 −i
,
i 0

β=

I2 02
02 −I2

σz =

!
1 0
;
0 −1

(3)

!
,

I2 is the 2×2 unit matrix. Solutions of the free Dirac equation (ψA ) are 4-component-spinor orbitals
ϕA,1 (r)ϕ (r)
 A,2 ψA (r) = ,ϕA,3 (r)ϕA,4 (r)(4)(5)

whose conjugate-transpose form reads
†ψA(r) = ϕ∗A,1 (r) ϕ∗A,2 (r) ϕ∗A,3 (r) ϕ∗A,4 (r) .

(6)

4
Let us remark that the scalar functions (ϕ) represent spatial orbitals because the spin is accounted by their arrangementin 4-component spinors. These solutions can be partitioned into a set of positive-energy 4-component spinors (EA > 0)
and a set of negative-energy ones (EA < 0), i.e. {ψA } = {ψR } ∪ {ψI }. From now on, we will designate the “positiveenergy spinors” (“negative energy spinors”) as PS (NS). Hence, the Dirac field is quantized as
X
X
X †
b
bb ψI (r) +ψ(r)=baA ψA (r) =dbR ψR (r),
(7)
I
A

I

R

where the sum has been decomposed into contributions involving electron(positron) annihilation(creation) operators.
These operators obey the usual anticommutation relations{baA , baB† } = δAB and {baA† , baB† } = {baA , baB } = 0.
B.

(8)

The Hamiltonian operator including an external nonlocal potential, charge sectors in Fock space, andBogoliubov transformations

Next, let us introduce a Hermitian nonlocal external potential that can be expressed in 4-component form as



nl
nl
nl
nl
vext,1,1(r, r′ ) vext,1,2(r, r′ ) vext,1,3(r, r′ ) vext,1,4(r, r′ )
v nl
nl
nl
nl
(r, r′ ) vext,2,2(r, r′ ) vext,2,3(r, r′ ) vext,2,4(r, r′ )


nl
vext
(r, r′ ) =  ext,2,1
,
nl
′
nl
′
nl
′
nl
vext,3,1 (r, r ) vext,3,2 (r, r ) vext,3,3 (r, r ) vext,3,4 (r, r′ )
nl
nl
nl
nl
(r, r′ ) vext,4,2(r, r′ ) vext,4,3(r, r′ ) vext,4,4(r, r′ )vext,4,1

(9)

∗
 nl
nl
(r′ , r) to ensure the hermiticity of the(r, r′ ) = vext,ν,µwith the matrix elements obeying the constraint vext,µ,ν
operator.
The bound particle Dirac equation that describes the states of a single particle (electron or positron) subject to anonlocal external potential as in (9) readsZb D (r)ψA (r) + dr′ vnl (r, r′ )ψA (r′ ) = EA ψA (r).
T
(10)
ext
R
†
The states ψ depend on the choice of nonlocal external potential and form an orthonormal set ( drψA(r)ψB (r) =δAB ).
This single-particle equation can be readily generalized to noninteracting many-particle systems. To that end, let
us recast the Hamiltonian to operate in Fock space asnl
b 0v = TbD + Vbext
H
Z
i Z
h
 nl ′
′′′bn1 (r, r ) + drdr′ Tr vext
(r , r)bn1 (r, r′ )= drdr δ(r − r )Tr TD (r)b
XZ


nl
=
drdr′ δ(r − r′ )TD,µ,τ (r) + vext,µ,τ
(r′ , r) nb1,τ,µ (r, r′ ),

(11)

µ,τ

where we have introduced the one-particle density matrix operator whose elements read ashinb1,τ,µ (r, r′ ) = N ψbµ† (r′ )ψbτ (r)X †XX † †
bb bb ϕ∗ (r′ )ϕJ,τ (r) +bb db ϕ∗ (r′ )ϕS,τ (r)=I J I,µI S I,µ
I,J

+

I

XX
R

J

dbRbbJ ϕ∗R,µ (r′ )ϕJ,τ (r) −

(12)

S

X

dbS† dbR ϕ∗R,µ (r′ )ϕS,τ (r)R,Sthat is defined using creation and annihilation Dirac field operators with normal ordering [51, 52] (using Wick’s
theorem) N [· · · ] of the elementary operators [53] bbI† , dbR† , bbI , and dbR w.r.t. the effective vacuum state (|0v ⟩). The
formal definition of the effective vacuum state requires a detailed discussion that we briefly summarize in the following

5lines. First we note that any Hamiltonian leads to a set of PS and NS. One could then define a bare vacuum state inwhich all spinors are unoccupied. This definition of the vacuum is problematic, however. Since NS are much lower inenergy than PS, electrons occupying PS would be able to decay into empty NS by emitting photons resulting in an
endless cascade of emission processes. To overcome this issue, Dirac proposed to fill all NS with electrons, introducing
the so-called Dirac sea[54]. This is also problematic as the energy of such a vacuum state in which all NS are occupied
is minus infinity. To avoid dealing with infinities, a reinterpretation was thus suggested in QED where all filled (byelectrons) NS are redefined as empty positronic spinors. Thus, the effective vacuum state, which sets the zero of theenergy scale, corresponds to the state where all PS do not contain electrons and all NS do not contain positrons. In
nl
bv = Hb 0 ≡ TbD ), we refer to its effective vacuum state as |00 ⟩ and to its spinor
the particular case when vext
= 04×4 (H000basis as ψA (r).
nl0b v (with vext
Let us assume that a given H
̸= 04×4 ) is built initially in the free particle ψA(r) basis, with the normal0vbordering taken w.r.t. free particle vacuum |00 ⟩. In this basis H0 is not diagonal and it does not commute with the
be = P bb †bb (Nbp = P db † db ). Moreover, Hb v does not commute with theelectron (positron) number operators N
0
I I I
R R R
b =N
be + Nbp . A useful operator that does commute with Hb v is, however, thetotal number of particles operator N
0
opposite charge operator [49, 55] that is defined asX †X †bb bb −
b=N
be − Nbp =Q
dbR dbR
(13)
I I
I

R

v
that provide the diagonal rep(see Appendix A for more details). We then consider using the eigensolutions ψAcfv with the tilde symbol indicating that the Hamiltonian is transformed with respect to its original
resentation H
0representation and has normal ordering w.r.t. its own vacuum |0v ⟩. For consistency we will also indicate this vacuumas |e0v ⟩ with the tilde symbol indicating that these quantities are transformed with respect to their original representations. In the current case |e0v ⟩ = |0v ⟩, but this will change when two-particle interactions are considered below. In
eA (r) = ψ v (r) basis the Hamiltonian commutes with the number of electrons (positrons) operators, the numberthe ψAof particles operator, and the opposite charge operator. Due to the diagonal representation, the pair creation andannihilation processes vanish; thence, the Hamiltonian conserves the particle number. Since both basis are orthonormal, they are related by a unitary transformation V = eκ [56], where the antihermitian matrix κ (formed by theparameters κAB ∈ C) can be chosen to contain all its diagonal entries equal to zero to avoid redundant phase-shiftsof the spinors. Therefore, the relationship between the two bases reads as
X
0eA (r) =ψψB(r)VBA .
(14)
B
†
†
c
be
eb , be b
e
fv operator is written in terms of the transformed operators (bThe H
I dR , bI , and dR ) that are given by
0
X
∗
b
e
aA = eκb b
aA e−bκ =b
aB VBA,(15)

B

where the spinor rotation operator (eκb ) is expressed as the exponential of an antihermitian operator κb with
X
κ
b=
κAB b
aA† b
aB
A,B

=

X

κIJ bbI†bbJ +

XX

I,J

+

I

XX
R

κISbbI† dbS†

(16)

S

κRJ dbRbbJ +JXκRS dbR dbS† .

R,SActually, the spinor rotation operator eκb corresponds to a Bogoliubov transformation mixing electron annihilationand positron creation operators [57, 58], producing modified electron and positron creation and annihilation operatorsthat are consistent with a transformed vacuum.
Recalling that the Hamiltonian commutes with the opposite charge operator in any basis, it is convenient to classify
its eigenstates within a Fock space that gathers together different particle-number sectors
(∞,∞)

F=

M

H(Ne ,Np ) ,(Ne ,Np )=(0,0)

(17)

6
L
where
designates the direct sum. Alternatively, the Fock space can also be decomposed into charge (Q = Ne − Np )sectors
F=

∞
M

HQ .

(18)

Q=−∞

b v , the wavefunction of a system containing a fixed number of electron (Ne ) and positrons
For a given Hamiltonian H
0(Np ) can be written as a single determinant (SD) [59]. The SD wavefunction reads as|Φ⟩ = bbI†1 · · · bbI†Ne dbR†1 · · · dbR†Np |00 ⟩,(19)

and is is antisymmetric under the exchange of particles (as required for systems formed by fermions). Furthermore,
the wavefunction in the basis where the Hamiltonian is diagonal can be written as a unitary rotation of the original
wave function by operator Vb ase = Vb |Φ⟩ = eκb |Φ⟩ = eκb bb † · · · bb † db † · · · db † |00 ⟩|Φ⟩I1
IN e R 1
RN p= eκb bbI†1 e−bκ · · · eκbbbI†Ne e−bκ eκb dbR†1 e−bκ · · · eκb dbR†Np e−bκ eκb |00 ⟩
b† b†
b†
b†
0v ⟩,
= ebI1 · · · ebINe deR1 · · · deRNp |e(20)

eA (r) basis aswith |e0v ⟩ = eκb |00 ⟩. Let us also rewrite the annihilation Dirac field operators in the ψ
b
ψ(r)
=

X

b
eA (r) =
e
aA ψ

Xb †
Xb
eb ψ
e
eR (r),deR ψ
I I (r) +
I

A

(21)

R

be†
a similar expression can be written for the creation Dirac field operator (ψ(r)). In summary, as in the nonrelativistic
context, the spinor rotation operator allows us to change the basis representation of the creation and annihilation
operators, which may be used to obtain a diagonal representation of a particular Hamiltonian. As the only difference
between Hamiltonians is in the nonlocal external potential, this diagonal representation always depends on the giveneA (r) in
nonlocal external potential. Finally, taking the Taylor expansion of eκb , it is easy to recognize that basis ψwhich the Hamiltonian is diagonal is given by the set of {κAB } parameters that lead to a stationary point [60] (i.e. aminimum, a maximum, or a saddle point) of the energy w.r.t. variations of these parameters.
C.

Vacuum polarization

From now on, we will assume that an initially given orthonormal basis ψA (r) is employed to express the Hamiltonian0
operator, this could be the ψA(r) basis, but any orthonormal basis can in principle be used as long as it is possible todistinguish between initial positive and negative energy states, e.g. by diagonalizing an initial Hamiltonian. To stress
this freedom of initial reference, we drop the symbol 0 from this definition as the choice of a free particle vacuum is only
one possibility. To avoid artefacts in the calculation of the kinetic energy due to basis set incompleteness, it is only
necessary to build this initial basis with adequate kinetic balance [61–63] between the bases for the large and small
components of the spinors. The initial effective vacuum state is now defined as |0⟩ with the initial normal ordering
taken w.r.t. this vacuum. The absence of tilde symbols on the basis functions and the vacuum hereby indicates thatthe Hamiltonian is non-diagonal in this basis. In the previous section we introduced the rotation operator (eκb ) anddefined the transformed effective vacuum|e0v ⟩ = eκb |0⟩(22)that provides a new reference for normal ordering the transformed operators. As we previously mentioned, the initiallycbb v operator can be rewritten in another basis (Hfv ) employing the operator n
defined H
e1,τ,µ (r, r′ ) (i.e. expressing the00b† b† bb
one-body density matrix operator in terms of the elementary operators ebI , deR , ebI , and deR ). Recalling that the normal
ordering is introduced to guarantee that the energy of the effective vacuum is zero, i.e.,cb v |0⟩ = ⟨e
fv |e
⟨0|H
0v |H
0
0 0v ⟩ = 0,

(23)

7
when a Hamiltonian is not evaluated with its reference effective vacuum state it leads to nonzero contributionscb v |efv |0⟩. Actually, the relationship between the noninteracting Hamiltonian operators can simply⟨e
0 v |H̸ 0 ̸= ⟨0|H
0 0v ⟩ =0
be written ascbv = H
fv + E
e0 ,
H
0
0
0

(24)

where

e00 = ⟨eb 0v |eE0v |H0v ⟩ =

Z

drdr′ Trhi′b D (r) + vnl (r′ , r) n
e vp
(r,
r
δ(r − r′ )T
)
,
ext
1

(25)

is known as the vacuum polarization [46] (vp) energy (also present in the interacting particle picture, see below). And
with the vacuum polarization one-particle density matrix defined asX
X†′eR (r)ψ
e† (r′ ) −e vpnψψR (r)ψR(r′ )(26)
1 (r, r ) =R
R

R

whose components read as′e n1,η,ν (r, r′ )|e
n
evp
0v ⟩
1,η,ν (r, r ) = ⟨0v |bX
X=
ϕeR,η (r)ϕe∗R,ν (r′ ) −ϕR,η (r)ϕ∗R,ν (r′ ).
R

(27)

R

Before proceeding, let us remark that vp effects are present whenever spinor rotations mixing the positive andnegative spinors are involved and the reference effective vacuum changes. When only one-body operators are present
in the Hamiltonian, changes on it can be accounted for by orbital relaxation. This is no longer possible when also
two-body operators are considered as will be discussed in the next section.

D.

The interacting particle Hamiltonian and the configuration-interaction vacuum

In the more complete description of relativistic electronic dynamics provided by QED also the electromagnetic field
is quantized such that the electron-electron (and electron-positron) interaction is mediated by the exchange of photons.
This allows for a proper account of retardation effects due to the photons finite velocity. In the present work we takea simpler approach and treat the electromagnetic field semi-classically with the approximate interaction provided by
the Coulomb–Breit interaction. This operator corresponds to the single-photon exchange electron-electron scattering
amplitude in QED evaluated with the zero-frequency limit of the photon propagator in the Coulomb electromagnetic
gauge [55]. Unfortunately, this approximation makes the present theory not Lorentz invariant, but it is the most widelyused interaction at present; it improves the description of the interaction w.r.t. nonrelativistic Coulomb interaction(i.e. |r1 − r2 |−1 ). Thus, the approximate interacting Hamiltonian reads as
b =H
b 0v + WcHZb v + 1 dr1 dr2 Tr [W(r1 , r2 )bn2 (r1 , r2 )]
=H
0
2
X Z
b 0v + 1=H
dr1 dr2 Wµ,ν,τ,η (r1 , r2 )bn2,τ,η,µ,ν (r1 , r2 ),2 µ,ν,τ,η

(28)

where
Wµ,ν,τ,η (r1 , r2 ) =11(αµ,τ · r12 ) (αν,η · r12 )δµ,τ δν,η − (αµ,τ · αν,η ) −
2
r12
2
2r12

(29)

with r12 = r1 − r2 and r12 = |r12 |; and the pair-density matrix operator is defined using creation and annihilation
Dirac field operatorshinb2,τ,η,µ,ν (r1 , r2 ) = N ψbν† (r2 )ψbµ† (r1 )ψbτ (r1 )ψbη (r2 ) ,(30)

8where the normal ordering is taken w.r.t. the effective vacuum state |0⟩. Due to the interaction, a SD wavefunctionb and representations of H
b are in general non-diagonal.
(Eqs. (19) and (20)) is no-longer an eigenstate of Hb
Let us start the search for the wavefunctions of H by noticing that the Fock space can be split into charge sectors,
b can be written as [64]such that the ground state energy of H
bv + Wc |Ψ⟩.
E = min ⟨Ψ|H
0(31)

|Ψ⟩∈HQbb
Thence, the wavefunctions must also be eigenstates
R of Q with eigenvalue Q (i.e. Q|Ψ⟩ = Q|Ψ⟩). Consequently, |Ψ⟩can be constrained to have a given charge Q (i.e. dr⟨Ψ|Tr [bn1 (r, r)] |Ψ⟩ = Q). Assuming that Ψ is normalized to 1,
a state that belongs to HQ for a particular charge sector Q ≥ 0 can be written (parameterized) as|Ψ⟩ = X

cI1 ...IQ bbI†1 · · · bbI†Q +I1 ,...,IQ

X

X

I1 ,...,IQ ,IQ+1

R1

cI1 ...IQ IQ+1 R1 bbI†1 · · · bbI†Q bbI†Q+1 dbR†1


X

+

X

cI1 ...IQ IQ+1 IQ+2 R1 R2 bbI†1 · · · bbI†Q bbI†Q+1 bbI†Q+2 dbR†1 dbR†2 + · · · |0⟩.

(32)

I1 ,...,IQ ,IQ+1 ,IQ+2 R1 ,R2Since we account for all contributions to all orders there is no need for spinor rotations, i.e. spinor rotations become
redundant, c.f. full configuration interaction (FCI) in the nonrelativistic context.
The charge sectors Q < 0 present the same structure as the Q ≥ 0 ones. Thus, without any loss of generality andfor practical reasons (i.e. as we are usually interested in electronic wavefunctions more than in positronic ones), fromnow on will assume that our wavefunctions belong to the latter one.
The Q = 0 case is especially interesting as it could be used for a redefinition of the vacuum state. Up to now,
we have considered an effective vacuum (|0⟩) that could be obtained for some effective one-body potential v. This is
not possible when we include an explicit two-electron interaction. Nevertheless, from the partition of the Fock space
into charge sectors, it is possible to define instead a CI vacuum wavefunction (|0CI ⟩) as a linear combination of stateswith arbitrary number of electron-positron pairs that belongs to the charge sector Q = 0. Thus, the CI vacuum
wavefunction reads as|0CI ⟩ = c0 +

X X
I1

R1

cI1 R1 bbI†1 dbR†1 +X X

cI1 I2 R1 R2 bbI†1 bbI†2 dbR†1 dbR†2 + · · · |0⟩,(33)

I1 ,I2 R1 ,R2

where spinor rotations are redundant because this wavefunction accounts for all vacuum contributions to all orders.
Let us remark that although the two-electron operator does not change, these CI coefficients vary with the nonlocal
external potential (e.g. when the molecular structure is changed). Only if this nonlocal external potential as well as
the charge (Q) is kept fixed, then normal ordering [51, 52] can in principle be defined w.r.t. this |0CI ⟩ vacuum state.
b CI ⟩ = 0) can then be omitted. In the usual situation, when the external potential is
The vp effects (i.e. ⟨0CI |H|0
allowed to change, vp effects cannot be neglected even when a CI vacuum is employed as the coefficients of CI vacuum
depend on the external potential and this then affects the normal ordering.
In the nonrelativistic limit, when only electrons are considered in the wavefunction, the effective vacuum and the
CI one coincide because in the reinterpretation of the NS the electrons are decoupled from the positrons and the latterdo not contribute to the wavefunction. Moreover, in that limit the external potential does not affect the definitionof the vacuum state; thence, the effective vacuum can be used independently of the interactions considered for theelectrons and for any different external potential employed.

E.

The 1-RDM and its natural orbital representation

Let us comment on some properties of the 1-RDM obtained from a wavefunction |Ψ⟩ given by Eq. (32). Using theone-particle density matrix operator Eq. (12) we define the matrix elements of the 1-RDM in the spinor basis asn1,ν,µ (r, r′ ) = ⟨Ψ|bn1,ν,µ (r, r′ )|Ψ⟩

9
=

X

+

XX

1

DIJ ϕ∗I,µ (r′ )ϕJ,ν (r) +

XX

I,J

R

I
1

I ∗
DR
ϕR,µ (r′ )ϕI,ν (r) +1

DIR ϕ∗I,µ (r′ )ϕR,ν (r)1

DSR ϕ∗R,µ (r′ )ϕS,ν (r),

R

X

(34)

R,S

I

I
where we have introduced the 1-RDM coefficients 1 DIJ = ⟨Ψ|bbI†bbJ |Ψ⟩, 1 DIR = ⟨Ψ|bbI† dbR† |Ψ⟩, 1 DR
= ⟨Ψ|dbRbbI |Ψ⟩, and†b
1 Rb
DS = −⟨Ψ|dS dR |Ψ⟩. In a nonrelativistic context, the coefficients arising from the expectation values ⟨Ψ|bbI† dbR† |Ψ⟩and ⟨Ψ|dbRbbI |Ψ⟩ do not contribute since the wavefunction is formed by a fixed number of particles and due to theabsence of positrons in this limit, the last contribution ⟨Ψ|dbS† dbR |Ψ⟩ then trivially vanishes.
Collecting the matrix elements we form the corresponding 4×4 matrix, and after trivial algebra we obtain the full
1-RDM as
ϕ∗I,1 (r′ )ϕJ,1 (r) ϕ∗I,2 (r′ )ϕJ,1 (r) ϕ∗I,3 (r′ )ϕJ,1 (r) ϕ∗I,4 (r′ )ϕJ,1 (r)
 ∗

′
∗
′
∗
′
∗
′
X
1 J ϕI,1 (r )ϕJ,2 (r) ϕI,2 (r )ϕJ,3 (r) ϕI,3 (r )ϕJ,2 (r) ϕI,4 (r )ϕJ,2 (r)n1 (r, r′ ) =DI  ∗
ϕI,1 (r′ )ϕJ,3 (r) ϕ∗I,2 (r′ )ϕJ,3 (r) ϕ∗I,3 (r′ )ϕJ,3 (r) ϕ∗I,4 (r′ )ϕJ,3 (r)I,J
ϕ∗I,1 (r′ )ϕJ,4 (r) ϕ∗I,2 (r′ )ϕJ,4 (r) ϕ∗I,3 (r′ )ϕJ,4 (r) ϕ∗I,4 (r′ )ϕJ,4 (r)
XX
XX
X
1 R
1 I
1 R
DI (· · ·)4×4 +DR (· · ·)4×4 +DS (· · ·)4×4
+
=

I

R

R

X

1

DIJ ψJ (r)ψI† (r′ ) +

XX

I,J

I

I

R,S
1

DIR ψR (r)ψI† (r′ ) +

R

XX
R

1

†
I
DR
ψI (r)ψR(r′ ) +IX

1†
DSR ψS (r)ψR(r′ ). (35)

R,S

P
B
A
form an Hermitian
The trace of the 1-RDM is equal to the charge Q, i.e. A 1 DA
= Q. Moreover, the coefficients 1 DA
matrix that can be diagonalized, which leads to the natural orbital [65] (NO) representation of the 1-RDM
n1 (r, r′ ) =X

nA χA (r)χ†A (r′ ),

(36)

A

where {nA } are known as the occupation numbers
P (ONs) that take values in the interval between [−1, 1] and the NOs
are written as 4-component spinors χA (r) = B ψB (r)UBA (the columns of the U matrix contain the eigenvectors
obtained from the diagonalization of the 1 D matrix).
In the particular case of a SD wavefunction (19) the (canonical) ψ spinors coincide with the NOs with ONs thatare then either 0 or ±1. Consequently, the 1-RDM for |Φ⟩ reads as′nSD
1 (r, r ) =NeX

ψI (r)ψI† (r′ ) −

I

F.

Np
X

†
ψR (r)ψR(r′ ).

(37)

R

Introducing ReRDMFT through the constrained-search formalism Considering the constrained-search formalism [66–68] we define the universal functional of the 1-RDM W n1 for
N -representable 1-RDMs (n1 ∈ DQ ) that come from a state |Ψ⟩ ∈ HQ [69],W [n1 ] =

min

c |Ψ⟩
⟨Ψ|W

|Ψ⟩∈HQ (n1 )c |Ψ [n1 ]⟩,= ⟨Ψ [n1 ] |W(38)

where HQ (n1 ) is the set of states |Ψ⟩ ∈ HQ that yield a constrained 1-RDM (n1 ), and |Ψ [n1 ]⟩ designates the statethat minimizes this energy contribution [70]. Thence, we may write the ground state energy functional of the 1-RDM
as

Z
h

i
′
′ b
nl
′
′
EQ = min W [n1 ] + drdr Tr δ(r − r )TD (r) + vext (r , r) n1 (r, r ) .
n1 ∈DQ(39)

10
At this point we would like to stress the similarity in structure of the expressions above with relativistic DFT, which
is based on local potentials and the opposite charge densityZn(r) = dr′ δ(r − r′ )Tr [n1 (r, r′ )] .
(40)
The disadvantage of relativistic DFT is, however, that the free Dirac operator must also be included into the constrained search expression in this case, i.e.,
F [n] =minc |Ψ⟩⟨Ψ|TbD + W|Ψ⟩∈HQ (n)c |Ψ [n]⟩,= ⟨Ψ [n] |TbD + W(41)

where HQ (n) is the set of states |Ψ⟩ ∈ HQ that yield a constrained opposite charge density (n), and |Ψ [n]⟩ designates
the state that minimizes this energy contribution. Consequently, the energy functional within the constrained search
densformalism for N -representable opposite charge densities (n ∈ DQ
) reads asZ
EQ = min F [n] + drvext (r)n(r) .
(42)
dens
n∈DQ

Returning to ReRDMFT, clearly, the advantage of the functional of the 1-RDM over the one based on the opposite
charge density is that all one-body interactions are an explicit functional of n1 . For instance, adding an energy
contribution arising from an external vector field (Aext ) to the above expression is straight forward: only the termnlδ(r − r′ )cα · Aext (r) needs to be added and fits the general form of vext
(r′ , r) in (9). Thus, the functional of the 1RDM contains more information than the functional of the current (i.e. j(r) = cψ † (r)αψ(r) = Tr [cαn1 (r, r)]), whichclearly shows the generality of the present theory and makes it interesting also for the description of magnetism inmolecular systems. The functional presented in (39) establishes relativistic reduced density matrix functional theory(ReRDMFT). In particular, when the 1-RDM is given in the NO representation we will refer to it as relativistic
natural orbital functional theory (ReNOFT). For practical purposes and to resemble nonrelativistic applications of
RDMFT, we will use ReNOFT in this work when building functional approximations. An alternative approach to theconstrained-search formalism is presented in the Appendix B, where we discuss the extension of Gilbert’s theorem [44]to the relativistic domain.
III.
A.

THE NO-PAIR APPROXIMATION AND ReRDMFT

Vacuum polarization effects at fermion-fermion interacting case

The Hamiltonian given in Eq. (28) is written in normal-ordering w.r.t. to the initial effective vacuum state |0⟩[71],bwhich leads to ⟨0|H|0⟩= 0. The relationship between the interacting Hamiltonian and its transformed counterpart(i.e. written in terms of transformed operators and basis) is given by the equation
vp
be b
b =H
e0 ,
H
+ Ve + E

(43)

where
b vp
b vp b vp
Ve = Ve H + Ve x
"
#
XZ
XZ
vp
b
=
dr1
dr2 Wµ,ν,τ,η (r1 , r2 )en1,η,ν (r2 , r2 ) ne1,τ,µ (r1 , r1 )
µ,τ

−

X Z Z

(44)

ν,η

b
dr1 dr2 Wµ,ν,τ,η (r1 , r2 )envp
e1,τ,ν (r1 , r2 ),
1,η,µ (r2 , r1 )nµ,ν,τ,η

e0 = ⟨ebe
e0 + 1E
0v |H|
0v ⟩ = E02Z

dr1 dr2 Tr [W(r1 , r2 )envp
2 (r1 , r2 )] ,(45)

where we have introduced the vp pair-density matrix (written in terms of its components) as
n
evp
evp
nvp
evp
nvp
2,τ,η,µ,ν (r1 , r2 ) = n
1,η,ν (r2 , r2 )e1,τ,µ (r1 , r1 ) − n
1,τ,ν (r1 , r2 )e1,η,µ (r2 , r1 ).

(46)

11
Let us highlight that the definition of the vacuum as an effective vacuum (|0⟩ or |e0v ⟩, i.e. at the mean-field level), leads
vp
vp
b
b
to only Hartree (Ve H ) and exchange (Ve x ) like terms for the interaction between the fermions (electrons or positrons)
e vp
present in the |Ψ⟩ and the polarization of the vacuum. Indeed, the n
1 can be regarded as the difference between thee states, i.e.,
1-RDM obtained from a SD wavefunction built from all positronic ψ
†
b†
e =b
0v ⟩,
|Φ⟩
deR1 · · · deRM |e(47)

and the same kind of function build with all positronic ψ states, i.e.,|Φ⟩ = dbR†1 · · · dbR†M |0⟩,(48)

where M is the total number of positronic states (i.e. the positronic states are fully occupied). Thus, the energy terms
that account for the interaction between particles-vacuum and the ‘vacuum-vacuum’ effects, resemble the contribution
arising from the frozen-core electrons in a multi-configuration self-consistent field calculations in the nonrelativistic
context (e.g. complete active space self-consistent field). Consequently, only Hartree and exchange like terms arepresent.
The effective vacuum state |e0v ⟩ that minimizes Eq. 45 can actually produce an energy that diverges to −∞ due
to infrared and ultraviolet divergences, which should be accounted for in a practical implementation of the present
theory. Nevertheless, these issues go beyond the scope of the present work; therefore, we assume that a proper
renormalization scheme is applied in order to keep everything finite in the rest of this work. Further details about
infrared and ultraviolet divergences and how renormalization should be applied can be found in Refs. 45, 72–79.
Finally, let us recall that when the normal ordering is taken w.r.t. the CI vacuum state (for a fixed nonlocal
external potential), vp effects do not need to be considered anymore even when spinor rotations are applied. In otherwords, only the CI vacuum leads to reference state that is invariant under spinor rotations. Such a definition of
vacuum can, however, not be used in practice as the number of CI coefficients involved is formally infinite, whichmakes the CI expansion unmanageable from the practical perspective. Furthermore, working with finite basis sets
still introduces a large number of CI coefficients (much larger than in non-relativistic Full CI calculations); thus, theuse of the CI vacuum becomes rapidly intractable also for finite basis. In the following, we introduce the so-called
no-pair approximation, where spinor rotations are required; therefore, they are always be accompanied by vp energy
contributions.

B.

The no-pair approximation including vacuum polarization ReRDFMT approach

The no-pair (np) approximation [48, 80] is a widely used simplification when dealing with relativistic calculations.
It relies on the fact that most physical and certainly chemical processes correspond to much smaller energy exchanges(e.g. X-rays that interact with the electronic density contain ∼ 150keV) than required to create electron-positron pairs(> 1MeV) [81]. States that contain one or more positron pairs do therefore hardly contribute to the wave functions
changes induced by such processes. Thus, in the description of common processes like chemical reactions or various
kinds of spectroscopic phenomena, pure electronic wavefunctions are generally employed (i.e. Ne = N and Np = 0).
Within the np approximation only the PS are populated while the NS remain empty. The np vp ground state energy
reads as
npvp
EN
=

min
e (N,0)|Ψ+ ⟩∈Hb + ⟩,⟨Ψ+ |H|Ψ

(49)

e (N,0) = eκb H(N,0) that is the set of states generated bywhere the minimization is over normalized states in the set He (N,0) can be written as
all spinor rotations of N -electron states. A state |Ψ+ ⟩ ∈ H
X|Ψ+ ⟩ = eκb
cI1 ...IN bbI†1 · · · bbI†N |0⟩I1 ,...,IN

=

X

b†
b†
cI1 ...IN ebI1 · · · ebIN |e
0v ⟩.

(50)

I1 ,...,INb
Since the state |Ψ+ ⟩ can also be written as |Ψ+ ⟩ = Pe+ |Ψ⟩, where |Ψ⟩ ∈ HN is an arbitrary state constrained to haveb†b
N negative charges, and Pe+ projects onto the N -electron Hilbert space constructed by the {ebI } operators. Thence,

12
the energy minimization procedure given by Eq. (49) can also be regarded as an optimization w.r.t. the projector,
where the energy minimization procedure implies finding an optimal vacuum state that depends on the Ne throughEqs. (44) and (45). Lastly, using the constrained search formalism [66–68] and Eq. (49) we can readily introduce
npvp-ReRDMFT, where vp effects are accounted at the effective QED level [82].
For making npvp-ReRDMFT more explicit, let us denote the np 1-RDM as
X′
1 eJ ee† (r′ ),n+
DI ψJ (r)ψ(51)
1 (r, r ) =
I
I,J
†
+
e
eb b
e J = ⟨Ψ+ |b
e vp
where 1 D
I bJ |Ψ+ ⟩. Making use of n1 and n
1 all one-body interactions of Eq. (49) (i.e. the electronic andIthe vp effects) can be readily evaluated asnl⟨Ψ+ |TbD + Vbext|Ψ+ ⟩ =

Z

drdr′ Trhb D (r) + vnl (r′ , r)δ(r − r′ )T
ext



vp
′
′
n+
1 (r, r ) + n1 (r, r )

i

.

(52)

vp
c
Also, the explicit vp contributions in terms of n+
1 and n1 to ⟨Ψ+ |W |Ψ+ ⟩ is known through Eqs. (44) and (45). Indeed,cf |Ψ+ ⟩ in terms of n+ is unknown and it needs to be approximated (seeonly the pure electronic contribution to ⟨Ψ+ |W
1
Appendix C for more details).

C.

The no-pair approximation ReRDFMT approachvp
e0 and ⟨Ψ+ |Vb
e |Ψ+ ⟩) is usually neglected in practical applications
The vp contribution to the total energy (i.e. Eof the np approximation. Redefining the normal ordering w.r.t. an effective vacuum state used whenever a spinor
rotation is applied, the np energy reads
np
EN
=

min
e (N,0)|Ψ+ ⟩∈Hbe⟨Ψ+ |H|Ψ
+ ⟩,

(53)

where the Hamiltonian is written in normal ordering w.r.t. the floating vacuum state |e0v ⟩. From Eq. (53) we may+
write, within the constrained-search formalism, the functional expression for N -representable np 1-RDMs (n+1 ∈ DN )e (N,0)that come from a state |Ψ+ ⟩ ∈ H
 
W np n+1 =mincf |Ψ+ ⟩
⟨Ψ+ |We (N,0) (n+ )|Ψ+ ⟩∈H
1

 +
  c
f
= ⟨Ψ+ n+1 |W |Ψ+ n1 ⟩,

(54)

 
e (N,0) (n+ ) is the set of states |Ψ+ ⟩ ∈ H
e (N,0) that yield a constrained 1-RDM (n+ ), and |Ψ+ n+ ⟩ designates
where H1
1
1
the state that minimizes this energy contribution. Let us remark that we have explicitly written the Hamiltonian
operator as a functional of the np 1-RDM to highlight that it changes during the minimization procedure. Actually,
during the optimization procedure the gap between the floating vacuum state and the N -electron ground state energy
is increased.
Consequently, the total energy functional that designates np-ReRDMFT can be defined as

np
EN
=


Z
h

i
 +
+
np
′
′ b
nl
′
′
min+ Wn1 + drdr Tr δ(r − r )TD (r) + vext (r , r) n1 (r, r ) .
+(55)

n1 ∈DN+The np approximationh i resembles the nonrelativistic result because the trace of the n1 matrix is equal to the number
e = Ne . The n+ matrix can also be expressed in the NO representation,of electrons, i.e. Tr 1 D
1′n+
1 (r, r ) =XeI (r)χ
e†I (r′ ),
nI χ

I

with

P

I nI = Ne which allows us to introduce np-ReNOFT.

(56)

13
To conclude this section, as it is mentioned in Ref. 45, the np approximation regains the concept of an N -electron
wavefunction,XeI (r1 ) ∧ · · · ∧ ψeI (rN ),cI1 ...IN ψ(57)
Ψ+ (r1 , r2 , . . . , rN ) =1N
I1 <...<IN

where ∧ denotes the normalized antisymmetrized tensor product.

D.

cf |Ψ+ ⟩ term as an explicit functional of the second-order reduced density matrix.
The ⟨Ψ+ |W

c
f operator. In the np
To conclude this section, let us comment on some properties of the expectation value of the W
approximation this term reads asXcf |Ψ+ ⟩ =⟨Ψ+ |W

2

KL
DIJ

Z

h
i
eL (r1 ) ⊗ ψeK (r2 ))(ψe† (r2 ) ⊗ ψe† (r1 ))dr1 dr2 Tr W(r1 , r2 )(ψ
I
J

(58)

I,J,K,Lb †b †b b
KL
where ⊗ denotes the tensor product, and 2 DIJ= 21 ⟨Ψ+ |ebI ebJ ebLebK |Ψ+ ⟩ is the np second-order reduced densitycf |Ψ+ ⟩ contribution is an explicit functional of the 2-RDM. In the np
matrix (2-RDM) element. Clearly, the ⟨Ψ+ |W
approximation, the properties of the np 2-RDM are the same as in the nonrelativistic case. In particular, the trace of
this matrix reads  X 2 IJ
Ne (Ne − 1)
Tr 2 D =
DIJ =
.
(59)
2
I,J

Similarly, this matrix must fulfill the same N -representability conditions [83, 84] as in the nonrelativistic case, which we
exploit in the following section to adapt/re-build some functional approximations for the relativistic problem. Finally,cf.
practical applications of relativistic quantum chemistry/physics normally retain only part of the contributions to W
Indeed, the Coulomb–Breit interaction (29) consists of two physically distinct contributions: one representing themagnetic interaction between the electrons, and another one representing the retardation due to the finite velocity of
the interaction [85]. The magnetic one was also discovered by Gaunt [86] and can be combined with the Coulombinteraction to yield
Wµ,ν,τ,η (r1 , r2 ) =1[δµ,τ δν,η − (αµ,τ · αν,η )] ,
r12

(60)

which will be used in the rest of this work (i.e. inserting it in Eq. (58)).
npvp
np
Before proceeding, let us stress that both energy expressions EN
and ENcan be minimized since the Hamiltonianis written using normal ordering. Hence, they are bounded from below and the state |Ψ+ ⟩ can actually be considered
as a projected state from the |Ψ⟩ ∈ HN .
E.

np-ReNOFT functional approximations by imposing Kramers’ symmetry

Since the np approximation is the most common initial starting point for relativistic calculations, we propose theconstruction of functional approximations
  in this framework.
The explicit dependence of the W n+functional in terms of the n+
1
1 is unknown (from now on we will drop the+super-index when referring to the np 1-RDM). In nonrelativistic NOFT the most accurate functionals [87–90] arebuilt by imposing the so-called N -representability conditions [83, 84, 91] and using only up to two different indices todefine 2-RDM elements, because the Hartree and the exchange contributions can be fully accounted for by using only
two indices in the NO representation. These approximations usually employ a restricted approach, whose equivalent
in the relativistic contexts corresponds to using Kramers’ pairing symmetry when external magnetic fields are notincluded [85], which leads to a restricted formulation of spin magnetic moments. Then, a pair of NOs forms a Kramers’bχbχei = χ
eī and K
eī = −χei , where we have defined the time-reversal operator aspair (i, ī) if they transform as K!
σy 02 bb
K = −i
K0 ,
(61)
02 σy

14
b 0 is the complex conjugation operator. From now on we form two subsets of PS by splitting all spinors
and where K
into Kramers’ pairs [92]. Actually, the proper definition of Kramers’ pairs is not unique because the pair (i, ī) is
degenerate; thus, any unitary transformation applied to this pair leads to an equivalent pair of spinors that is still a
Kramers’ pair. In the non-relativistic limit, this corresponds exactly to the arbitrariness in the orientation of the spinquantization axis of the restricted formalism. These subsets are labeled with lowercase barred and unbarred indexes.
Indeed, in practical calculations containing M PS, we have M/2 Kramers’ pairs and thus, M/2 (un)barred PS [93].
beImposing the Kramers’ symmetry we may rewrite the H
operator (including only the Coulomb and the Gauntinteractions) in the NOs representation as [55]
X
X h
be KR
b+ + 1H=
hii X⟨ij|kl⟩b
x++
x−−
x++
− ⟨īj|α ⊙ α|kl⟩bx−−ii
ik,jl − ⟨ij|α ⊙ α|kl⟩bik,jl + ⟨īj|kl⟩bīk,jlīk,jl2
ii,j,k,lii 1 X h
++
−−
−−
¯l⟩b
¯l⟩b
⟨
īj|k
x
−
⟨
īj|α
⊙
α|k
x
+ ⟨ij|k̄l⟩b
x++
−
⟨ij|α
⊙
α|
k̄l⟩b
x
+
īk,j l̄īk,j l̄ik̄,jlik̄,jl4
i,j,k,l
i
1 X h
−−
++
−−
¯l⟩b
¯l⟩b
⟨īj̄|kl⟩b
x++
−
⟨
ī
j̄|α
⊙
α|kl⟩b
x
+
⟨ij|
k̄
x
+
−
⟨ij|α
⊙
α|
k̄
x
,
īk,j̄līk,j̄lik̄,j l̄ik̄,j l̄8

(62)

i,j,k,lwhere α ⊙ α = αx ⊗ αx + αy ⊗ αy + αz ⊗ αz and the integrals are defined asZb D (r) + vnl (r, r′ ) χe†I (r′ ) δ(r − r′ )TeI (r),
hII = drdr′ χ
ext


Z
e†I (r) ⊗ χe†J (r′ ) χeK (r) ⊗ χeL (r′ )′ χ,⟨IJ|KL⟩ = drdr|r′ − r|Ze†I (r) ⊗ χe†J (r′ ) α ⊙ α χeK (r) ⊗ χeL (r′ )′ χ,⟨IJ|α ⊙ α|KL⟩ = drdr|r′ − r|

(63a)
(63b)
(63c)

†
† †
b †b
b †b
b †b
b †b
2 s1
eb b
e bs1 s2 = (1+s1 Tbik )(1+s2 Tbjl )b
eb b
e b
e b
e
b s = (1+sTbij )bX
bsJL,IK
with Tbijebi ebj = ebj̄ ebī and Tbijebī ebj = −ebj̄ ebi [94].
i bj , x
I bJ bL bK = x
ij
IK,JL
Let us remark that the + and − signs placed in Eq. (62) do not refer to positive or negative energy spinors; they only
b
np
1 s2
e KR |Ψ+ ⟩.
b s and x. Using the above Hamiltonian, the np energy reads as EN= ⟨Ψ+ |H
bsIJ,KL
enter in the definition of X
IJ
For keeping the notation as concise as possible, we will drop the KR superscript on the Hamiltonian from now on.
Within nonrelativistic NOFT, most of the approximations rely on the usage of only up to two different indices inthe electron-repulsion integrals to account for Hartree, exchange, and (some) correlation effects. Actually, retaining
up to two indices most of the (nonrelativistic) NOFT approximations are able to retrieve the so-called nondynamic
electron correlation energy, but they fail to account for the dynamic one (that can be accounted using different
strategies [90, 95, 96]). Following the same strategy for np-ReNOFT and retaining only up to two different indices in
Eqs. (63b) and (63c) the np energy in the NO representation, Eq. (55), reads as

np
EN
≈

X

hii (ni + nī ) +

i

−

X

X

2


ij
īj
Dij
+ 2 Dīj+ 2 Diij̄j̄ + 2 Dīīj̄j̄ Jij

i,j
2

i,j

ij
īj
Dij
− 2 Dīj− 2 Diij̄j̄ + 2 Dīīj̄j̄



G
Jij
+

X h

2

ji
Dij
+ 2 Dj̄īj̄īG
Kij − Kij

ii,ji X h2 j j̄ 2 īi 
i
1 X h2 j ī 2 ij̄ 2 īj 2 j̄i GDīj + Dj̄i + Dj ī + Dij̄ Lij − LG
+
D
+
D
K
+
K
+
ij
ij
ij
j̄j
iī
2 i,j
i̸=j
i
1 X h2 iī 2 j j̄ 2 īi 2 j̄j −
Dj̄j + Dīi + Dj j̄ + Diī Lij + LG
,
ij
2

(64)

i̸=j

G
G
where Jij = ⟨ij|ij⟩, Jij= ⟨ij|α ⊙ α|ij⟩, Kij = ⟨ij|ji⟩, Kij= ⟨ij|α ⊙ α|ji⟩, Lij = ⟨īj|j ī⟩ (notice that Lii = 0 [85]),
and LG
=
⟨
īj|α
⊙
α|j
ī⟩.
Interestingly,
some
NOFT
approximations(i.e. GNOF/PNOFx approximations) employedij
in nonrelativistic calculations [88, 97] already use L integrals to account for correlation effects among opposite-spin
electrons.
The properties of the matrix elements of the np 2-RDM are the same as of the nonrelativistic 2-RDM ones (i.e.
the matrix elements present the same symmetry and antisymmetry properties upon exchange of indices, they enter

15
similarly in the evaluation of the N -representability conditions, etc.); hence, we have adopted the same strategy thatis used in the nonrelativistic context to develop np(vp)-ReRDMFT functional approximations. Indeed, keeping only
up to two different indices and using the natural orbital basis permits us to account for all Hartree and exchange interactions (now also including the exchange among ‘opposite spin’ terms in np(vp)-ReNOFT[np(vp)-ReRDMFT]). But,
it also facilitates the search for analytic constraints on the approximated 2-RDM matrix elements upon evaluationof the D-, Q-, and G-N -representability conditions. Indeed, using the same strategy as in the nonrelativistic context
leads to relativistic functional approximations that correspond to a generalization of their nonrelativistic counterparts(like the Dirac–Hartree–Fock method represents a generalization of the nonrelativistic Hartree–Fock method). Consequently, with this strategy we can also ensure that the correct nonrelativistic limit is retrieved with the functionals
proposed in the following sections.

1.

The Dirac–Hartree–Fock approximation

The np Dirac–Hartree–Fock (DHF) approximations uses a SD wavefunction to approximate |Ψ+ ⟩ (occupying only
PS). Actually, the DHF energy can be readily written from Eq. (64) for ‘spin-compensated’ systems where ni = nīand defining the matrix elements of the 2-RDM with the SD approximation (2 DSD )KL SD(2 DIJ) =nI nJ(δIK δJL − δIL δJK ).
2(65)This approximation produces the correct symmetry and antisymmetry properties of the 2-RDM [98]. Thence, the np
DHF energy reads as
E np,DHF = 2=2

Xhii ni +

X

i

i,j

Ne /2

Ne /2

X

hii +

i




G
ni nj 2Jij − Kij − Kij− Lij − LG
ij

X


G
2Jij − Kij − Kij− Lij − LG
.
ij

(66)

i,j

In the last expression we have use ni = nī = 1 for the PS that form the SD wavefunction, and 0 otherwise. An efficientalgorithm to optimize this functional has been recently proposed by Sun et al. [99]. Moreover, this functional reducesδδ
to the one provided by Hafner for a purely Coulomb electron-electron interaction (i.e. Wµ,ν,τ,η (r1 , r2 ) = µ,τr12ν,η ) from
Ref. 100.
In the following, we introduce two paths to build ReNOFT approximations that are also used in the nonrelativistic
context. The functionals presented in this work are proposed for ‘spin-compensated’ systems, i.e. the NOs forming
Kramers’ pairs always show the same ONs (ni = nī ). The first family of approximations is based on mimicking theexchange-correlation hole, while the second one on imposing N -representability conditions. Nevertheless, both paths
approximate the 2-RDM elements as functions of the ONs, i.e.
2

KL
KL
DIJ
= 2 DIJ(nI , nJ , nK , nL )(67)

and does not take into account the actual shape of the spinors.

2.

f (nI , nJ )-functional approximations

The family of nonrelativistic f (nI , nJ )-functional approximations is introduced as an attempt to separate theelectron-electron interactions into Hartree and exchange-correlation effects [98]. These approximations account for
exchange-correlation effects by approximating the so-called exchange-correlation hole. The exchange-correlation contribution is accounted replacing the product nI nJ in the second term on the r.h.s. of the DHF energy by a f (nI , nJ )function that attenuates the exchange; thus, accounts not only for exchange but also for exchange-correlation effects.
In this context, the relativistic extension of the f (nI , nJ )-functional approximations (i.e. the 2RDM elements) reads
as
2

KL
DIJ

approx.

=

nI nJ
f (nI , nJ )
δIK δJL −δIL δJK ,
2
2

(68)

In Table II we have collected some f (nI , nJ ) functions that lead to some representative functional approximations.

16

The functionals given by Eq. (68) retrieve their nonrelativistic counterparts. Although, in nonrelativistic NOFT
ji
and X 2 Dīj̄j̄ī ); therefore, theapproximations the f (nI , nJ ) functions only affect the ‘same-spin’ terms (i.e. X 2 Dij
n n

j ī
) could be defined as − i2 j (or using any other function), which may lead
‘opposite-spin’ ones (i.e. X 2 Dij̄ij̄ and X 2 Dīj
to an extended definition of f (nI , nJ )-functional approximations in the relativistic scenario.
Let us briefly introduce the origin of the functionals presented in Table II that have their origin in the nonrelativistic
context. The Müller, Buijse and Baerends (MBB) functional was introduced independently by Müller and by Buijse and Baerends [101–103], it produces an approximated functional which fulfills the sum rule, (i.e. Tr 2 D = Ne (N2e −1) ).
This functional was derived from the requirement of minimal violation of the Pauli principle and from the analysis of
Fermi and Coulomb holes.
In the case of the Power functional, the α parameter was first proven to have to be α ≥ 0.5565 to produce admissible
densities (i.e. solutions that are: stable with respect to the corresponding Euler equations, N -representable and whose
electron-electron interaction energy satisfies the Lieb–Oxford bound [109]). This functional is mostly used in solidstate physics and, according to Sharma et al. [105], the MBB functional overcorrelates the electrons and the Powerfunctional mediates with the overcorrelation of the electrons by using a parameter α > 1/2. In some applications, thebound α ≥ 0.5565 is not fixed and values like 0.53 are employed [110].
The Goedecker and Umrigar (GU) and the Marques and Lathiotakis (MLSIC) functionals are examples of approxii
imations developed to remove the so-call self-interaction [107, 108] terms (i.e. they remove the nonphysical 2 Dii
and
2 īī
Dīī elements). The former corrects these terms on the MBB functional, while the latter corrects a previous version
of a similar functional also based on a Padé approximant expression. The MLSIC parameters were optimized toreproduce the G2 test.
The energy expression for the f (nI , nJ ) functionals reads as

E rel-x = 2

X

hii ni +iX

2ni nj Jij −i,jX

f (ni , nj )





G
Kij − Kij
+ Lij − LG
ij

(69)

i,j

with x =MBB, Power, GU, and MLSIC.

3.

Approximations built imposing N -representability conditions

A second family of functionals approximations is based on the reconstruction of the 2RDM elements imposing theD, Q, G N -representability conditions [83, 88, 98, 111].
From Eq. (64), we recognize the 2-RDM matrix elements that must be approximated. Indeed, the 2-RDM matrix
elements can be organized in blocks of the form

TABLE II. Representative f (nI , nJ )-functional approximations.
Functional
MBB
Power
GU

f (nI , nJ )√nI nJ(nI nJ )α
√
nI nJ
nI nJ
0 +a1 nI nJ
nI nJ a1+b
1 nI nJ

MLSIC

case

Ref.

all
101–103all
104–106
I ̸= J
107
otherwiseI ̸= JnI nJ
otherwisea0 = 1298.78,
a1 = 35114.4, andb1 = 36412.2

108

17







2
D=






2

2 ij
Dij
ij , Dji



0

0

0
!

!
2

0

2

0

, 2 Dīj
Dīi
j̄j
īj
Dij̄jī , 2 Djījī

02i̸=j !
2i̸=j

Dīi, 2 Dījj j̄
j īDijīj̄ , 2 Djj īī00i̸=j !

0
i̸=j



2

0

Dīīj̄j̄ , 2 Dīj̄j̄ī







.






(70)





2 ij
2 īj̄ 2 īj̄
We may identify three decoupled blocks containing different spinor combinations. The 2 Dij
Dīj̄ , Dj̄ ī
ij , Dji andblocks resemble the αα and ββ blocks of the nonrelativistic case, while the large middle block corresponds to theαβαβ, αββα, etc. terms. In the nonrelativistic limit, the off-diagonal terms of the middle block do not contribute to
the energy as they cancel upon integration over spin (e.g. the block 2 Dīi= 0). Proposing functional approximationsj j̄
is facilitated by the introduction of the so-called cumulant matrix, which can be defined as
Λ = 2 D − 2 DSD .

(71)

Then, the auxiliar matrices ∆ and Π are defined to approximate the cumulant matrix (as in the nonrelativistic
context). These auxiliar matrices facilitate the evaluation of the N -representability conditions. Upon evaluation of
the D-, Q-, and G- N -representability conditions (see the Appendix D for more details), and dividing the spinor
space Ω into subspaces mutually disjoint Ωp (see Fig. 1 for more details) we arrive to the relativistic version of the
GNOF/PNOFx (x = 5, 7) functional approximations [90, 95, 112] that can be written as
Ne /2

E rel-GNOF/PNOFx =

X
p=1

Ne /2

Ep +

X

Eqp .

(72)

q̸=p

The first sum accounts for all intra-subspace contributions and reads asXEp =ni (2hii + Jii + JiiG + LG
ii )
i∈Ωp

+

X

G
G
Πintrai,j (Kij + Kij + Lij + Lij ),(73)

i,j∈Ωp
i̸=j

where
( √
− ni nj ,Πintra=√i,j+ ni nj ,i or j ≤ Ne /2
i, j > Ne /2,(74)

the second sum accounts for inter-subspace contributions (Eba ) that can be defined as
X X


G
Eqp =
ni nj 2Jij − (Kij − Kij) − (Lij − LGij )i∈Ωq j∈Ωp

+

X X

G
G
Πinteri,j;p,q (Kij + Kij + Lij + Lij ),(75)

i∈Ωq j∈Ωp

p
inter
inter
where Πinteri,j;p,q = 0 in rel-PNOF5, Πi,j;p,q = − ni hi nj hj in rel-PNOF7 [113], Πi,j;p,q = −4ni hi nj hj in relPNOF7s [114] [115], and
q

p
d d

n
n
−
n
h
n
h
−
ndi ndj , i j
i i j jqp
Πinteri,j;p,qndi ndj − ni hi nj hj + ndi ndj ,



0,

i > q, j = p or i = q, j > p
i > q, j > p
i = q, j = p

(76)

18

FIG. 1. Partition of the spinor space Ω into subspaces mutually disjoint Ωp . While arrows refer to α or β spin in thenonrelativistic context, in this work they correspond to electrons occupying spinors that are Kramers’ pairs. Let us remarkthat some NOs may remain uncoupled (a.k.a. virtual spinors) like the orange ones in this scheme.

√


ni hd
with ndi = hp p , hp = 1 − np , and hdp = hp exp −(hp /(0.02 2))2 in GNOF [116]. Each electron pair is describedby an even number of PS forming Kramers’ pairs. Moreover, the rel-GNOF/PNOFx functionals presented in thiswork produce 2-RDM elements that are real, but complex elements could be introduced by taking other functions toapproximate the Πintra and Πinter sub-matrices.

As we did for the f (nI , nJ )-functionals, we briefly introduce the origin of the nonrelativistic PNOFx approximations
whose relativistic adaptations are presented in this work for completeness. PNOF5 was developed to produce thecorrect electron distribution upon bond cleavage processes[112]. This functional was initially introduced in the perfectpairing approach, where each sub-space Ωp contained only two NOs. Although, it was later extended beyond theperfect-pairing approach to allow more orbitals in each sub-space. The lack of inter-subspace interaction in PNOF5
lead to the later versions of PNOFx functionals (like PNOF7). The PNOF7 functional accounts for the interaction
among pairs and has been proven account for the so-called nondynamic correlation energy, but it misses the dynamic
one [34, 35, 113]. In an attempt to account for both types of electron correlation the GNOF functional was recently
proposed [116]; it has proven to be capable of accounting for an important part of the so-called dynamic correlationenergy. Though, it still misses important contributions that lead to weak interactions (like ones present in van der
Waals interactions and Hydrogen bonds). Finally, let us mention that in the limit when the ONs tend to 0 and 1,
all these functional approximations retrieve the HF energy expression. The same holds for their relativistic versions,that lead to the np DHF expression for the total energy (66) in this limit.

19
F.

Properties of rel-GNOF/PNOFx functionals.

For ‘spin-compensated’ two-electron systems (see Appendix E) the rel-GNOF/PNOFx functionals introduced inthis work reduce toX
X√
rel-GNOF/PNOFx
G
n1 ni (K1i + K1i+ L1i + LGE2e−=ni (2hii + Jii + JiiG + LG
1i )ii ) − 2
i

+

i
i̸=1

X√

G
ni nj (Kij + Kij+ Lij + LG
ij ),

(77)

i,j
i̸=j,i̸=1,j̸=1where the label 2e− indicates that it is for a two-electron system (also used in the Appendix E) and the PS are ordered
in descending order w.r.t. their ON (and we have labeled as 1 the one with the largest ON). This functional can beproven to be the relativistic extension of the Fixed-Phases functional [117], which is based on the NO representationof the (singlet) wavefunction for two-electron systems [87] (see the Appendix E for more details).
The nonrelativistic PNOF5 energy functional expression is known to be equivalent to a constrained version of theenergy expression of an antisymmetrized product of strongly-orthogonal geminals (APSG) wavefunction [51, 118]. In
this section, we prove that the rel-PNOF5 functional is also equivalent to a (relativistic) APSG wavefunction within
the np approximation. To that end, let us introduce the N -electron wavefunction
Ne /2

Y

ΨAPSG
(r1 , r2 , ..., rN ) = Ab+Ξp (r2p−1 , r2p ),

(78)

p=1

where
1 X
ei (r) ⊗ χeī (r′ ) − χei (r′ ) ⊗ χeī (r)]Ξp (r, r′ ) = √ci [χ2 i∈Ωp(79)√
is a geminal wavefunction written in the NO basis including phase factors (i.e. ci = ni eiζi , see Appendix E for
more details), and the Ab operator
produces the antisymmetrized product of geminal wavefunctions. The geminalR
wavefunctions are normalized drdr′ Ξ†p (r, r′ )Ξp (r, r′ ) = 1 for each electron pair, wherei
1 X ∗h †ei (r) ⊗ χe†ī (r′ ) − χe†i (r′ ) ⊗ χe†ī (r) .
Ξ†p (r, r′ ) = √
ci χ
(80)
2 i∈Ωp
Imposing the strong-orthogonality conditionZ∀pq
dr1 dr2 Ξ†p (r1 , r2 )Ξq (r1 , r2 ) = δpq ,(81)

we arrive to the following expression for the 1-RDM
Ne /2

n1+,APSG (r, r′ ) =

X X

h
i
ei (r)χ
e†i (r′ ) + χeī (r)χ
e†ī (r′ )
ni χ

(82)

p=1 i∈ΩpThus, the relativistic APSG wavefunction permits us to write the energy functionalNe /2ei }] = 2
E APSG [{ζi }, {ni }, {χXni hii +

i
Ne /2

+

X


p=1

X

p̸=q

i,j
i∈Ωp ,j∈ΩqGni nj 2Jij − (Kij − Kij) − (Lij − LG
ij )

(83)





X

X

√
G

+ Lij + LGei(ζj −ζi ) ni nj (Kij + Kijij )i̸=j∈ΩpP
which can be optimized subject to the conditions: a) ∀i ni ∈ [0, 1], and b) Ne = 2 i ni . Clearly, the energy functional
given by Eq. (83) is a lower bound w.r.t. the rel-PNOF5 functional. Noticing that in rel-PNOF5 functional the phaseseiζi = ±1 are fixed (see Eq. (74)) while in E APSG they are parameters, we can readily conclude that in E APSG there
is more flexibility during the optimization procedure that can lead to a lower energy.

20
IV.

CLOSING REMARKSIn this work we have introduced ReRDMFT at three different levels of theory. At the first level we have presented
ReRDMFT including electron-positron pair creation/annihilation processes. At this level, the trace of the 1-RDM(n1 (r, r′ )) is the charge Q; the minimum of the energy (EQ ) is attained for the CI ground state wavefunction (|ΨQ ⟩).
This minimum is guaranteed due to the normal ordering procedure applied and it permits us to introduce ReRDMFT
through the constrained search formalism. In principle, when the nonlocal external potential is fixed, the normal
ordering is taken w.r.t. the CI vacuum, and the total energy of the the system of interest is referenced w.r.t. thisvacuum, its energy is invariant under spinor transformations (rotations). The energy as functional of the 1-RDM for all
one-body interactions is presented, but functional approximations are required for the particle-particle interactions.
The development of approximations for the particle-particle interactions will remain an open task, but a proper
definition of these interactions would lead to a more precise description of the relativistic problem as it would also bevalid for systems including positrons (e.g. electron-positron pairs).
Secondly, we have discussed the effect of taking the normal ordering w.r.t. the effective vacuum state |e0v ⟩ andhow this leads to the so-called vp when a spinor rotation is applied. These effects become unavoidable when the np
approximation is adopted. Therefore, at the second level we have introduced npvp-ReRDMFT. We have discussed
the importance of the NS spinors when spinor transformations are applied in the energy minimization procedure.
vp′′
At this level, two 1-RDMs play a crucial role (i.e. n+
1 (r, r ) and n1 (r, r )). Within npvp-ReRDMFT, QED effects
are taken at the mean-field level and the advantage of using a formulation based on the 1-RDMs becomes evident
vp′′
because all vp effects as well as electron-vacuum interactions are explicit functionals of n+
1 (r, r ) and n1 (r, r ). Only+′
the functional of the electron-electron interaction in terms of n1 (r, r ) remains unknown. Fortunately, the functional
approximations presented in this work (as part of the next level) are also valid for npvp-ReRDMFT.
Thirdly, the np-ReRDMFT is the last level of theoretical background introduced in this work. In np-ReRDMFT
vp effects are neglected and a floating effective vacuum state is employed as reference. Within this framework theconcept of an N -electron wavefunction is recovered. When the Hamiltonian preserves time-reversal symmetry it ispossible to exploit Kramers’ pairing symmetry; with Kramers’ pairs we have written the energy expression using only
up to two different indices. This expression allowed us to recognize the 2-RDM elements that need to be approximated as functions of the occupation numbers. In the end, the approximate 2-RDM elements proposed in this work
were properly adapted from the two major families of functional approximations used in the nonrelativistic context.
Subsequently, some properties of the functionals based on the inclusion of N -representability conditions were also
discussed. The performance of np-ReRDMFT approximations is an open question that will be addressed in future
works. Nevertheless, the computational cost of using np(vp)-ReRDFMT will increase w.r.t. its nonrelativistic counterpart, simply because we are forced to always use complex algebra and a spinor representation. The computational
cost for the optimization should be similar to the cost required by a multiconfigurational-self-consistent field approach(specially for the optimization of the spinors). With the advantage that in np(vp)-ReRDMFT the optimization of the
CI vector is replaced by the (usually cheaper) optimization over occupation numbers.
Finally, let us note that the analysis of the so-called generalized Pauli constraints [119–126] (GPC) of the ONs
obtained when using some nonrelativistic functional approximations has attracted some attention [127] in the last few
years. These conditions serve to approach to the possibility of obtaining pure N -representable 1-RDMs. They areparticularly important in systems that are not ‘spin-compensated’ (i.e. that in our case do not preserve time-reversal
symmetry [128–130]); thus, they are not contemplated in this work but they might be explored in future studies
for treating ‘spin-uncompensated’ systems. Furthermore, pure N -representability conditions can also be imposed to
the 2-RDM [130], and together with the D-, Q-, and G-conditions could also lead to improve results especially instrongly-correlated systems as Mazziotti suggested [130]; we may expect the same improvement for the relativistic
functionals presented in this work.

ACKNOWLEDGEMENTS

M.R.-M. acknowledges the European Commission for an Horizon 2020 Marie Sklodowska-Curie Individual Fellowship (891647-ReReDMFT). M.R.-M. would also like to acknowledge Prof. M. Piris for the fruitful discussions.

21
Appendix A: Commutation relation between the Hamiltonian and the opposite charge operator

b v and the opposite charge operator Qb let us first rewrite theIn order to prove the commutation relation between H
0
Hamiltonian operator in the 4-component spinor basis as
X †XX † †
bb bb HIJ +bb db HIR
bv =
H
(A1)
0
I J
I R
=

I,J

I

R

X

X

dbS† dbR HRS

dbRbbI HRI −

R,I

R,S

where
Z

h
i
†
b D + vnl (r′ , r) ψB (r).
drdr′ ψA(r′ ) δ(r − r′ )T
ext

HAB =

(A2)

The evaluation of the commutator between opposite charge operator and first term of Hamiltonian lead us to
X

h
i X
i
Xh †
bb bb , db † db
HIJ bbI†bbJ , bbK† bbK −
HIJ
I J
R R

I,J,K

I,J

=

X

(A3)

R



HIJ bbI†bbJ bbK† bbK − bbK† bbK bbI†bbJ

I,J,K

=

X



HIJ bbI†bbJ bbK† bbK − bbI†bbJ bbK† bbK + bbI†bbK δJK − bbK† bbJ δIK

I,J,K

=0
h
i
the commutator bbI†bbJ , dbR† dbR = 0 because positronic and electronic creation and annihilation operators do commute.
P
The proof that the fourth term of the Hamiltonian (i.e. − R,S dbS† dbR HRS ) also commutes with the opposite charge
operator follows closely the one presented for the first term; we therefore omit its detailed description. Next, thecommutation relation between the second term of the Hamiltonian and the opposite charge operator reads as
XX
I,J

R

XX
I,J

I

I,J


 XX


HIR bbI† dbR†bbJ†bbJ − bbJ†bbJ bbI† dbR† −HIR bbI† dbR† dbS† dbS − dbS† dbSbbI† dbR†
I

R,S

 XX



HIR bbI† dbR† dbS† dbS − bbI† dbR† dbS† dbS − bbI† dbS† δSR
HIR bbI† dbR†bbJ†bbJ − bbI† dbR†bbJ†bbJ − bbJ† dbR† δIJ −
I

R

XX
I

(A4)

R,S

R

XX

=

h
i XX
h
i
HIR bbI† dbR† , bbJ†bbJ −HIR bbI† dbR† , dbS† dbS

HIR

−bbI† dbR† + bbI† dbR†



R,S

= 0,

R

which proves that the second term also commutes with the Hamiltonian. Finally, the commutator between the third
term of the Hamiltonian and the opposite charge operator leave us with
XX
I,J

R

XX
I,J

=

I


 XX


HRI dbRbbI bbJ†bbJ − bbJ†bbJ dbRbbI −HRI dbRbbI dbS† dbS − dbS† dbS dbRbbI
I

R,S


 XX


HRI dbRbbI bbJ†bbJ − dbRbbI bbJ†bbJ + dbRbbJ δIJ −HRI dbRbbI dbS† dbS − dbRbbI dbS† dbS + dbSbbI δRS

R

I

XX
I

(A5)

R,S

R

XX
I,J

i XX
i
h
h
HRI dbRbbI , bbJ†bbJ −HRI dbRbbI , dbS† dbS

HRI



R,S


dbRbbI − dbRbbI = 0,Rthat completes the demonstration of the commutation relation between the Hamiltonian and the opposite charge
operator.

22
Appendix B: Extending Gilbert’s theorem to the relativistic domain.

In 1975, Gilbert extended the Hohenberg and Kohn theorems [131–133] for external nonlocal potentials. In thiswork we extend this theorem to the relativistic domain. Based on the existence of a lower bound of the energy (ensured
by writing the creation an annihilation operators in normal ordering) we may define the ground state wavefunction|Ψ⟩ ∈ HQ (see Eq. (32)) as the one that minimizes the energy for a particular Hamiltonian (Eq. (28)) and charge
sector Q of the Fock space.
We generalize Gilbert’s theorem for relativistic nondegenerated ground states in the next theorem.
b = TbD + Vb nl + W
c , whose nondegenerated ground state wavefunction
Theorem 1 Given a relativistic Hamiltonian H
extreads as ΨQ (for a particular charge sector Q). There exists a one-to-one correspondence between ΨQ and its 1-RDM(n1 (r, r′ )).
1,nl
2,nl
Proof. Assume that two-different nonlocal external potentials Vbext
and Vbext
that differ in more than a constant lead
1,nl
12
b 1 Ψ1 = (TbD + Vbext
to the nondegenerated ground states ΨQ and ΨQ (for a charge sector Q). That is to say, H
+
Q
1
1
2
2
1
2
1
′
2
′
c
b
W )ΨQ = E1 ΨQ and H2 ΨQ = E2 ΨQ . The 1-RDMs of ΨQ and ΨQ are n1 (r, r ) and n1 (r, r ), respectively. If the two
2,nl
1,nl
b2 − Hb 1 = Vbext
Hamiltonians only differ in the nonlocal external potential (i.e. H− Vbext), by Rayleigh-Ritz variational
principle we can writeb 1 |Ψ1Q ⟩ < ⟨Ψ2Q |H
b 1 |Ψ2Q ⟩ = E12
E1 = ⟨Ψ1Q |H(B1)

b 2 |Ψ2Q ⟩ < ⟨Ψ1Q |H
b 2 |Ψ1Q ⟩ = E21 ,E2 = ⟨Ψ2Q |H(B2)∆E = (E12 − E1 ) + (E21 − E2 ) > 0(B3)∆E = (E12 − E1 ) + (E21 − E2 )b2 − Hb 1 |Ψ1Q ⟩ + ⟨Ψ2Q |H
b1 − Hb 2 |Ψ2Q ⟩= ⟨Ψ1Q |H
Z
 nl

nl
= − drdr′ Tr v2,ext(r′ , r) − v1,ext(r′ , r)× n21 (r, r′ ) − n11 (r, r′ ) > 0

(B4)

and

so that

and

1,nl
2,nl
Since we imposed that Vbext
and Vbext
are different, the corresponding 1-RDMs (n11 (r, r′ ) and n21 (r, r′ )) must also be
different. Therefore, there exists a one-to-one correspondence between the 1-RDM and the wavefunction.

ΨQ ↔ n1 (r, r′ ).

□

(B5)

Consequently, there exists a functional of the energy
EQ = E [n1 ]for

n1 (r, r′ ) ∈ DQ .

(B6)This functional retrieves its nonrelativistic counterpart for a fixed number of electrons (i.e. Q = Ne , see below).
Furthermore, this functional leads to the existence of ReRDMFT and it complements the functional introduced
through the constrained-search formalism.
The no-pair vacuum polarization framework′
It is possible to establish a one-to-one correspondence between |Ψ+ ⟩ and n+
1 (r, r ) at the no-pair vacuum polarization
approximation level of theory. Using the same interacting Hamiltonian operator but restricting the wavefunctions
to be given by Eq. (50) (or Eq. (57)), similar arguments to the above presented lead to the desired one-to-one′
correspondence (i.e. |Ψ+ ⟩ ↔ n+
1 (r, r )). Finally, neglecting vacuum polarization effects, we may obtain the energyfunctional
 
+
(B7)
ENe = E n+forn1 ∈ DN1ethat allows us to introduce np-ReRDMFT and clearly retrieves the nonrelativistic limit.

23
Appendix C: The npvp energy functional

The E npvp functional of the np 1-RDM and the vp 1-RDM reads as
Z

h

i
′
b D (r) + vnl (r′ , r) n+ (r, r′ ) + n
e vp
drdr′ Tr δ(r − r′ )Text
1
1 (r, r )
#
"
XZ
XZ
vp
dr1
dr2 Wµ,ν,τ,η (r1 , r2 )en1,η,ν (r2 , r2 ) ne++
1,τ,µ (r1 , r1 )



e vp
E npvp n+
=
1 ,n
1

µ,τ

−

(C1)

ν,η

X Z Z

dr1 dr2 Wµ,ν,τ,η (r1 , r2 )envp
n+1,η,µ (r2 , r1 )e1,τ,ν (r1 , r2 )µ,ν,τ,η

1
+
2

Z

 
f +
dr1 dr2 Tr [W(r1 , r2 )envp
2 (r1 , r2 )] + W n1 ,

 
  c
 
f n+ = ⟨Ψ+ n+ |W
f |Ψ+ n+ ⟩ is unknown and needs to be approximated. Finally, thewhere the explicit form of W
1
1
1



npvp
usual np approximation corresponds to E np n+=En+1
1 ,0 .
Appendix D: The N -representability conditions within the np approximation.

The so-called N -representability conditions of the 2-RDM aim to ensure that this matrix is associated to a wavefunction (i.e. 2 D ↔ |Ψ+ ⟩). They also serve to propose a systematic way to build approximations for the 2-RDM
matrix elements in terms of the 1-RDM ones. Following the reconstruction procedure proposed in Ref. 88, let us
define the matrix elements of the cumulant matrix as [134]λij,kl = −∆ij∆ij
δik δjl +δil δjk ,
2
2

(D1a)

and for the middle blockΠīi,k̄k∆ījδik δjl +δij δkl ,22Πiī,kk̄
∆ij̄
λij̄,kl̄ = −δik δjl +δij δkl ,22Πiī,k̄k∆ij̄
λij̄,k̄l =δil δjk −δij δkl ,22λīj,k̄l = −

(D1b)
(D1c)
(D1d)

and
λīj,kl̄ =Πīi,kk̄∆ījδil δjk −δij δkl .
2
2

(D1e)

In the i ̸= j case, the ∆ sub-matrices with a fixed i and j values must be Hermitian; thence, ∆∗ij̄ = ∆j̄i . From theantisymmetry properties of the 2-RDM, we obtain that ∆j̄i = ∆ij̄ , which makes these matrix elements be real. Thus,the ∆ sub-matrices must be symmetric. The hermiticity of the 2-RDM makes us to impose Πīi,kk̄ = Π∗kk̄,īi . And, theantisymmetry of the 2-RDM imposes Πīi,kk̄ = Πīi,k̄k = Πiī,k̄k = Πiī,kk̄ .
Using the above definitions for the SD contribution and the cumulant matrix elements, the 2-RDM elements canbe approximated asni nj − ∆ij
ni nj − ∆ijδik δjl −δil δjk ,22nī nj̄ − ∆īj̄nī nj̄ − ∆īj̄
2 k̄l̄
Dīj̄ =
δik δjl −δil δjk ,22Πiī,kk̄ni nj̄ − ∆ij̄
2 kl̄
Dij̄ =
δik δjl +δij δkl ,22Πīi,k̄kni nj̄ − ∆j̄i
2 l̄k
Dj̄i =
δik δjl +δij δkl ,
2
2
2

kl
Dij
=

(D2a)
(D2b)
(D2c)
(D2d)

24
2

kl̄
=−
Dj̄i

Πīi,kk̄ni nj̄ − ∆j̄i
δik δjl −δij δkl ,
2
2

(D2e)

2

Dil̄kj̄ = −Πiī,k̄kni nj̄ − ∆ij̄
δik δjl −δij δkl .
2
2

(D2f)

and

The antisymmetry properties of the 2-RDM elements implies that Eqs. (D2a), (D2b), and (D2c) provide enoughinformation to build the rest of 2-RDM matrix elements. In contrast to the nonrelativistic approach, the 2-RDM
matrix elements given by Eqs. (D2e) and (D2f) contribute to the energy; we must approximate them in terms of the∆ and Π matrices.
The D-, Q-, and G- N -representability conditions

The so-called N -representability conditions evaluated at the np approximation level are the same as in the nonrelativistic approach; they are associated with the positive semidefinite character of the following Hermitian matrices(see Refs. 84 and 98 for more details)
KL
=
DIJ

1
b †b †b b⟨Ψ+ |ebI ebI ebLebK |Ψ+ ⟩ ,
2

(D3)

QKL
IJ =

1
b b b †b †⟨Ψ+ |ebI ebJ ebL ebK |Ψ+ ⟩ ,
2

(D4)

GKL
IJ =

1
b †b b †b⟨Ψ+ |ebI ebJ ebL ebK |Ψ+ ⟩ .
2(D5)

In the nonrelativistic approach, these matrices can be split in terms of different spin-blocks for the evaluation of thepositive semi-definite character.
• The D-Condition:
Eq. (D3) tests the positive semidefinite character of the 2-RDM. From the block structure of the 2-RDM
2 ij(see Eq. (70)), let us first focus on the block (2 Dij
ij , Dji ), which using Eq. (D2a) gives us the condition thatfor i ̸= j there is a set of [(M/2) × (M/2 − 1)]/2 two-by-two sub-blocks of the form[135]!
(ni nj −∆ij )(n n −∆ )− i j2 ij
2,(D6)(n n −∆ )(ni nj −∆ij )− i j2 ij
2which upon diagonalization produce the eigenvalues 0 and 2(ni nj − ∆ij ). Since the D matrix must be positive
semidefinite, we can conclude that∆ij ≤ ni nj .

(D7)

Similarly, the third block produces the same type of eigenvalues and conditions (i.e. ∆īj̄ ≤ nī nj̄ ). Finally, forthe middle block and i ̸= j we have a set of [(M/2) × (M/2)] two-by-two blocks of the form!
(ni nj̄ −∆ij̄ )(ni nj̄ −∆ij̄ )
−
2
2
,
(D8)
(n n −∆ )(ni nj̄ −∆ij̄ )− i j̄2 ij̄
2whose eigenvalues are 0 and 2(ni nj̄ − ∆ij̄ ), and to fulfill the D-condition the following constraint must hold∆ij̄ ≤ ni nj̄ .

(D9)Also, for the middle block we have one large [(M/2) × (M/2)] block whose off-diagonal elements are given by
all the Π matrix off-diagonal elements. The diagonal part of this block contains the i = j case, which includes
the ∆ contribution. This large block is basis set dependent and does not lead to analytic expressions of any of
the auxiliar matrices; thence, we can not extract any information (constraints) from it.

25• The Q-Condition:
Using the anticommutation rules for the creation and annihilation operators, we may rewrite Eq. (D4) as

1
I
J
I
J
δJL (δIK − 1 DK
) − δIL (δJK − 1 DK
) + δJK 1 DL
− δIK 1 DL
2
12 IJ
= DKL + [δJL δIK hK − δIL δJK hK + δJK δLI hL − δIK δJL hL ] ,
2

2 IJ
QKL
IJ = DKL +

(D10)
(D11)

where we have defined hI = 1 − nI , and taken the 1-RDM in the NO representation in the last expression. It
is straightforward to recognize that the Q matrix shows the same block structure as the D matrix. Hence, wehave one large block that depends on the Π and ∆ matrices (and on the basis set), and two-by-two blocks fori ̸= j of the form!
(1−ni −nj +ni nj −∆ij )(1−ni −nj +ni nj −∆ij )
−
2
2
,
(D12)
(1−ni −nj +ni nj −∆ij )(1−ni −nj +ni nj −∆ij )−22(1−ni −nj̄ +ni nj̄ −∆ij̄ )2(1−ni −nj̄ +ni nj̄ −∆ij̄ )−2(1−ni −nj̄ +ni nj̄ −∆ij̄ )2(1−ni −nj̄ +ni nj̄ −∆ij̄ )
2

−

!
,

(D13)

that upon diagonalization lead us to the following conditions∆ij ≤ hi hj

and∆ij̄ ≤ hi hj̄ .

(D14)• The G-Condition:
As we did for the Q-condition, we may rewrite Eq. (D5), using the anticommutation relations of the creation and annihilation operators, as
1
2 KJ
1 K
GKL
IJ = − DIL + δJL DI
2
1
KJ= −2 DIL+ δJL δIK nI ,
2

(D15)
(D16)

where in the last expression we have used the 1-RDM in the NO representation. For i ̸= j the first block2 ij 2 ij
Dij , Dji (and also for the third block) the G matrix is formed by a large basis set dependent block that
does not provide any extra analytical constraint. But it also contains 1 × 1 blocks with elements of the formij
−2 Dij+ 12 ni that lead to a condition of the form
∆ij ≥ −hi nj ,(D17)which is easy to satisfy with ONs between 0 and 1 and ∆ij ≥ 0. Using the matrix elements of the central blockof the 2-RDM let us introduce the Gcentral matrix that reads as

!
!
īj
īj
 Gj̄i , Gīj
i̸=j !
Gcentral = 


Gij
, Gjj
j̄ ī
īī
i̸=j

Gījij̄ , Gījjī



,


i̸=j ! Gijj̄ī , Gjj īī

(D18)

i̸=j

whose elements can be built using! Eq. (D15). It!is easy to recognize that this block can be split into three disjoint
sub-blocks (i.e.

Gīj, Gīj
j̄i
īj

,
i̸=j

Gjij̄ī , Gjj īī, and the rest). The latter sub-block is basis set dependent andi̸=jdoes not lead to new information about the auxiliar matrices. On the contrary, the former sub-block provides
new information. This block leads to M/2 eigenvalues of the form−(ni nī − ∆iī ) − Πīi,īi + nī
,
2

(D19)

26
which fulfill the N -representability condition if we let for example Πīi,īi = nī and ∆iī = ni nī . Nevertheless, thissub-block also produces for i ̸= j two-by-two sub-blocks of the form
!
!
Gīj
Gīj
1 ∆īj + nī hj −Πj̄j,īi
j̄i
īj
=
,
(D20)
2
−Πīi,j̄j ∆j̄i + nj̄ hi
Gj̄i
Gj̄i
j̄i
īj
whose eigenvalues obtained upon diagonalization read as1nī hj + ∆īj + ∆j̄i + nj̄ hi
4
1q 2
4Πīi,j̄j + (nī hj − nj̄ hi )2 ;
±
4

(D21)

where we have used the symmetry properties of the Π matrix. The eigenvalue with the positive square rootis greater or equal to zero and it satisfies the requirement. On the other hand, the eigenvalue with a negativesquare root leads to some bounds for the elements of the Π matrix. These bounds depend on the ∆ matrix andthey can be written asnī hj + ∆īj + ∆j̄i + nj̄ hi

2≥ 4Π2īi,j̄j + (nī hj − nj̄ hi )2 ,nī hj (2∆īj + 2∆j̄i + 4nj̄ hi ) + ∆2īj + ∆2j̄i + 2∆j̄i nj̄ hi + ∆īj (2∆j̄i + 2nj̄ hi ) ≥ 4Π2īi,j̄j ,which in the particular case when ∆īj = ∆j̄i = 0 lead toqnī hi nj̄ hj ≥ |Πīi,j̄j |.

(D22)
(D23)

(D24)

TABLE III. Definition of the Π matrix elements, where ndi = ni hdp /hp with i ∈ Ωp , hp = 1 − np , and√
hdp = hp exp −(hp /(0.02 2))2 .
functional Πintra; i, j ∈ Ωp , i ̸= j Πinter; i ∈ Ωp , j ∈ Ωq , p ̸= q
Ref.
ij
ij
√
− ni nj i = p or j = p
PNOF5 √
0
all
136
ni nj
otherwise√p− ni nj i = p or j = p
all95,113
PNOF7 √
− ni nj hi hj
ni nj
otherwise√− ni nj i = p or j = p
PNOF7s √−4ni nj hi hj
all
90
ni nj
otherwiseqp√− ni nj i = p or j = p ndi ndj − ni nj hi hj − ndi ndj i > p and j = q or i = p and j > q
q
p
√
116
GNOF
nn
otherwise nd nd − n n h h + nd ndi > p, j > q
i

j

i

j

i

j

0

i j

i

j

otherwise

From this analysis, we conclude that at the np-ReNOFT level the ansatz proposed in Eqs. (D1a)-(D1c) provides
equivalent eigenvalues for the D-, Q-, and G-conditions as the nonrelativistic context. Indeed, only the definitionof the auxiliar matrix Π has slightly changed w.r.t. its nonrelativistic counterpart. In principle, in the relativistic
approach four indices are required to define this matrix (instead of the two indices used in nonrelativistic NOFT)
ij̄
because now the elements of the form 2 Dījcontribute to the energy. Nevertheless, when the Π matrix is real and theONs for barred and unbarred states are the same (i.e. ni = nī ), only two indices are needed; thus, the Eq. (D24) canbe rewritten as in the nonrelativistic approachpni hi nj hj ≥ |Πi,j |,(D25)with Πi,j = Πīi,j̄j .
To ensure that the rel-GNOF/PNOFx functionals used in this work recover their nonrelativistic counterparts inthe nonrelativistic limit, we impose the following constraints:

27
1. Partition the PS space into subspaces {Ωp } (see Fig. 1 for more details).
2. Let ∆ij = ∆ij̄ = ∆īj = ∆īj̄ .
3. Define the diagonal terms Πi,i = ni and ∆ii = n2i to also satisfy Eq. (D19).
4. Set ∆ij = ni nj if i, j ∈ Ωp and ∆ij = 0 otherwise (this selection satisfies the above inequalities).
And we arrive to the rel-GNOF/PNOFx functionals presented in this work. The only missing term is the definitionof the Π matrix elements. In Table III we have collected different definitions, which are also employed in thenonrelativistic context; that (as we have seen) are also valid in the relativistic np approximation scenario. Let us
remark that the Π matrix is separated into two contributions. One contribution formed by indices belonging to thesame Ωp subspace (named as Πintra ) that accounts for intra-subspace interactions; the second contribution where
the indices belong to different Ω subspaces (denoted as Πinter ). Using these terms in the definitions of the 2-RDM
elements (Eqs. D2a-D2f), we arrive to the approximated (reconstructed) 2-RDM matrix. Inserting this matrix elements
in Eq. (64) we arrive to the functionals presented in Eq. (72).

Appendix E: The relativistic Fixed-Phases functional.

In 1956 Löwdin and Shull proved that the configuration interaction coefficients become simple functions of the√ONs (i.e. cI = ± nI ) when NOs are employed to express the two-electron wavefunction [87]. Thence, in this work,
we extend this result to the relativistic case. Let us start our discussion introducing the most general relativistic
two-electron wavefunction [137]X
Ψ2e−
CIJ ψI (r1 ) ⊗ ψJ (r2 ),(E1)+ (r1 , r2 ) =I,J

where the coefficients CIJ ∈ C and form an antisymmetricmatrix C (i.e. CIJ = −CJI and CII = 0). From thePnormalization of the wavefunction we have that I,J |CIJ |2 = 1. Applying the Carlson–Keller expansion [138] (a.k.a.
the Schmidt decomposition [139]) to the wavefunction (E1) we arrive to
1 X 
Ψ2e−
c̃ĩ ϑĩ (r1 ) ⊗ ϑ˜ĩ (r2 )+ (r1 , r2 ) = √2ĩ− ϑĩ (r2 ) ⊗ ϑ˜ĩ (r1 ) ,(E2)where the index ĩ runs only over half of the PS, and we have introduced the ϑ spinors that form Schmidt pairs (ĩ, ˜ĩ).
Note that these Schmidt orbitals are only defined up to a unitary transformation within the degenerate subspace.
The 1-RDM of the wavefunction given by E2 reads asihXn+,2e−(r, r′ ) =|c̃i |2ĩ ϑĩ (r)ϑ†ĩ (r′ ) + ϑ˜ĩ (r)ϑ˜† (r′ ) ,
(E3)
1
ĩ

ĩ

which is already given in its diagonal representation and from which we find |c̃i |2 = nĩWhen the Hamiltonian operator preserves time-reversal symmetry (i.e. it is time-independent and does not contain
external magnetic fields) it commutes with the Kramers’ operator [140]. This implies that we can choose all oureigenstates to satisfy [141]
b 2e−
Ψ2e−
+ (r1 , r2 ) = KΨ+ (r1 , r2 ).

(E4)which has the right symmetry for real CI coefficients. This means that the degenerate configurations in (E2) transform
into each other. Since the Schmidt orbitals are only defined up to a unitary transformation with the degenerate
subspace, we can use this freedom to transform the Schmidt pairs to Kramers pairs. So the Schmidt pairs can be
replaced with Kramers pairs
1 X ei (r1 ) ⊗ χeī (r2 )ci χ
Ψ2e−+ (r1 , r2 ) = √2 i
ei (r2 ) ⊗ χeī (r1 ) ,
−χ

(E5)

28
where
ci =

√

ni e−iζi .

(E6)

The phases ζi = kπ with k ∈ Z to preserve the Kramers’ symmetry.
The 1-RDM is given byn1+,2e− (r, r′ ) =Xhiei (r)χ
e†i (r′ ) + χeī (r)χ
e†ī (r′ ) .
ni χ

(E7)

i

Evaluating the electron-electron interaction using Eq. (E5) we arrive to
c
f 2e−
⟨Ψ2e−
+ |W |Ψ+ ⟩ =X √n i n j
i,j

2

e

i(ζj −ζi )Z
dr1 dr2 Tr

1(I16×16 − α ⊙ α)r12e†j̄ (r1 ))ei (r1 ) ⊗ χeī (r2 ))(χe†j (r2 ) ⊗ χe†j̄ (r2 )) − (χei (r1 ) ⊗ χeī (r2 )(χ
e†j (r1 ) ⊗ χ× (χe†j̄ (r1 ))ei (r2 ) ⊗ χeī (r1 ))(χe†j (r2 ) ⊗ χe†j̄ (r2 )) + (χei (r2 ) ⊗ χeī (r1 ))(χe†j (r1 ) ⊗ χ
− (χ

i

,

(E8)

where I16×16 is the 16×16 identity matrix. In this form, we notice that the energy minimization procedure also
implies the optimization w.r.t. the phases ({ζi }). Consequently, fixing the phases as in the nonrelativistic case andexploding Kramers’ symmetry on the integrals we arrive to the energy expression given by Eq. (77). Therefore,
the rel-GNOF/PNOFx functionals are equivalent to the relativistic version of the Fixed-Phases functional for ‘spincompensated’ two-electron systems.
Finally, let us show that the relativistic Fixed-Phases functional contains its nonrelativistic counterpart. To that
end, let us first focus on the NOs forming the Kramers’ pairs, which are built as four component spinors, i.e.
ϕei,1 (r)ϕe (r)
 i,2 ei (r) =  eχ(E9)ϕi,3 (r)ϕei,4 (r)and ∗
−ϕei,2 (r) ϕe∗ (r)bχeī (r).
ei (r) =  e∗i,1  = χK
−ϕi,4 (r)ϕe∗i,3 (r)(E10)

It is easy to prove that the tensor product of 4-component NOs contains singlet and triplet contributions, but
we can introduce some constraints to remove the triplet contribution. To that end, let us assume that the socalled small component terms of the NOs are negligible (i.e. ϕei,3 (r) and ϕei,4 (r) → 0) [142]. Then, adding theeeefollowing constraint ϕei,2 (r) = aϕei,1 (r) with a ∈ C we mayhwrite thei components as ϕi,1 (r) =h ℜϕi,1 (r)i + iℑϕi,1 (r)and ϕei,2 (r) = aℜϕei,1 (r) + aiℑϕei,1 (r), where ℜϕei,1 (r) = Real ϕei,1 (r) and ℑϕei,1 (r) = Imaginary ϕei,1 (r) . Thence, wearrive to the following expression for the NOs
ei (r) =χℜϕei,1 (r) + iℑϕei,1 (r)aℜϕei,1 (r) + aiℑϕei,1 (r)

!
(E11)

and
eī (r) =χ!
−a∗ ℜϕei,1 (r) + a∗ iℑϕei,1 (r),ℜϕei,1 (r) − iℑϕei,1 (r)(E12)

where a∗ is the complex-conjugate of a and we have already omitted the small component terms.
Assuming that the scalar spinors are real (i.e. ℑϕei,1 (r) = 0) and inserting the NOs in Eq. (E5) we arrive to [143]

29

1
2e−
Ψ+
(r1 , r2 ) = √

2

"
X
i

ci

!
!
!
!#
1
−a∗
1
−a∗
e
e
e
e
ℜϕi,1 (r1 )ℜϕi,1 (r2 )− ℜϕi,1 (r2 )ℜϕi,1 (r1 )
,
a
1
a
1

which in the a = 0 case allows us to introduce the usual nonrelativistic spin functions α =!
1and β =
0

(E13)

0
1

!
[144].

Inserting the spin functions may rewrite the wavefunction in a simplified notation as
1 X
ci ℜϕei,1 (r1 )ℜϕei,1 (r2 )(αβ − βα)Ψ2e−+ (r1 , r2 ) = √
2 i

(E14)

that corresponds to the singlet wavefunction in the nonrelativistic limit with the {ℜϕei,1 } being the scalar NOs.
Indeed, from Eq. (E14) the nonrelativistic Fixed-Phases functional can be readily obtained; therefore, the relativistic
Fixed-Phases functional presented in this work contains its nonrelativistic counterpart.

[1] J. Laerdahl, K. Faegri Jr, L. Visscher, and T. Saue, J. Chem. Phys. 109, 10806 (1998).
[2] P. Tecmer, A. S. P. Gomes, U. Ekström, and L. Visscher, Phys. Chem. Chem. Phys. 13, 6249 (2011).
[3] D. Samultsev, Y. Y. Rusakov, and L. Krivdin, Russ. Chem. Bull. 64, 2756 (2015).
[4] A. C. Castro, H. Fliegl, M. Cascella, T. Helgaker, M. Repisky, S. Komorovsky, M. Á. Medrano, A. G. Quiroga, and
M. Swart, Dalton Trans. 48, 8076 (2019).
[5] E. Tiesinga, J. Klos, M. Li, A. Petrov, and S. Kotochigova, New J. Phys. 23, 085007 (2021).
[6] E. Valderrama, E. V. Ludeña, and J. Hinze, J. Chem. Phys. 106, 9227 (1997).
[7] E. Valderrama, E. V. Ludeña, and J. Hinze, J. Chem. Phys. 110, 2343 (1999).
[8] E. Valderrama, J. Mercero, and J. Ugalde, J. Phys. B: Atom. Molec. Phys. 34, 275 (2001).
[9] J. Cioslowski, Phys. Rev. A 43, 1223 (1991).
[10] E. Ramos-Cordoba, P. Salvador, and E. Matito, Phys. Chem. Chem. Phys. 18, 24015 (2016).
[11] C. L. Benavides-Riveros, N. N. Lathiotakis, and M. A. Marques, Phys. Chem. Chem. Phys. 19, 12655 (2017).
[12] E. Ramos-Cordoba and E. Matito, J. Chem. Theory Comput. 13, 2705 (2017).
[13] M. Via-Nadal, M. Rodrı́guez-Mayorga, and E. Matito, Phys. Rev. A 96, 050501 (2017).
[14] M. Via-Nadal, M. Rodrı́guez-Mayorga, E. Ramos-Cordoba, and E. Matito, J. Phys. Chem. Lett. 10, 4032 (2019).
[15] L. Gagliardi and B. O. Roos, Chem. Phys. Lett. 331, 229 (2000).
[16] L. Gagliardi, B. Schimmelpfennig, L. Maron, U. Wahlgren, and A. Willetts, Chem. Phys. Lett. 344, 207 (2001).
[17] L. Gagliardi and B. O. Roos, Nature 433, 848 (2005).
[18] L. Gagliardi and B. O. Roos, Chem. Soc. Rev. 36, 893 (2007).
[19] S. Knecht, Ö. Legeza, and M. Reiher, J. Chem. Phys. 140, 041101 (2014).
[20] S. Battaglia, S. Keller, and S. Knecht, J. Chem. Theory Comput. 14, 2353 (2018).
[21] T. Fleig, L. K. Sørensen, and J. Olsen, Theor. Chim. Acta (Berlin) 118, 347 (2007).
[22] A. Ghosh, R. K. Chaudhuri, and S. Chattopadhyay, J. Chem. Phys. 145, 124303 (2016).
[23] E. Eliav and U. Kaldor, in Recent Progress in Coupled Cluster Methods (Springer, 2010) pp. 113–144.
[24] A. V. Oleynichenko, A. Zaitsevskii, L. V. Skripnikov, and E. Eliav, Symmetry 12, 1101 (2020).
[25] Y.-B. Tang, B.-Q. Lou, and T.-Y. Shi, Phys. Rev. A 96, 022513 (2017).
[26] A. Shee, T. Saue, L. Visscher, and A. Severo Pereira Gomes, J. Chem. Phys. 149, 174113 (2018).
[27] J. F. H. Lew-Yee, M. Piris, and J. M. del Campo, J. Chem. Phys. 154, 064102 (2021).
[28] J. Wang and E. J. Baerends, Phys. Rev. Lett. 128, 013001 (2022).
[29] R. López-Sandoval and G. Pastor, Phys. Rev. B 66, 155118 (2002).
[30] N. Lathiotakis and M. A. Marques, J. Chem. Phys. 128, 184103 (2008).
[31] S. Sharma, J. Dewhurst, S. Shallcross, and E. Gross, Phys. Rev. Lett. 110, 116403 (2013).
[32] K. Pernal, New J. Phys. 17, 111001 (2015).
[33] I. Mitxelena, M. Piris, and M. Rodrı́guez-Mayorga, J. Phys. Condens. Matter 29, 425602 (2017).
[34] I. Mitxelena and M. Piris, J. Phys. Condens. Matter 32, 17LT01 (2020).
[35] I. Mitxelena and M. Piris, J. Chem. Phys. 152, 064108 (2020).
[36] A. J. Cohen, P. Mori-Sánchez, and W. Yang, Science 321, 792 (2008).
[37] N. Q. Su, C. Li, and W. Yang, PNAS 115, 9678 (2018).
[38] M. Seidl, P. Gori-Giorgi, and A. Savin, Phys. Rev. A 75, 042511 (2007).
[39] A. D. Becke, J. Chem. Phys. 138, 074109 (2013).
[40] A. D. Becke, J. Chem. Phys. 138, 161101 (2013).
[41] J. Schmidt, C. L. Benavides-Riveros, and M. A. Marques, Phys. Rev. B 99, 224502 (2019).

30[42] C. L. Benavides-Riveros, J. Wolff, M. A. Marques, and C. Schilling, Phys. Rev. Lett. 124, 180603 (2020).
[43] T. Ohsaku, S. Yamanaka, D. Yamaki, and K. Yamaguchi, Int. J. Quant. Chem. 90, 273 (2002).
[44] T. L. Gilbert, Phys. Rev. B 12, 2111 (1975).
[45] J. Toulouse, SciPost Chem. 1, 002 (2021).
[46] P. Chaix and D. Iracane, Journal of Physics B: Atomic, Molecular and Optical Physics 22, 3791 (1989).
[47] P. Chaix, D. Iracane, and P.-L. Lions, Journal of Physics B: Atomic, Molecular and Optical Physics 22, 3815 (1989).
[48] M. H. Mittleman, Phys. Rev. A 24, 1167 (1981).
[49] T. Saue and L. Visscher, in Theoretical chemistry and physics of heavy and superheavy elements (Springer, 2003) pp.
211–267.
[50] H. A. Kramers, Proc. Acad. Amst 33 (1930).
[51] W. Kutzelnigg, Theor. Chem. Acc. 103, 182 (2000).
[52] W. Kutzelnigg, K. Shamasundarb, and D. Mukherjee, Mol. Phys. 108, 433 (2010).
PP[53] The spinor-components annihilation field operators are defined as ψbµ (r) = I b
bI ϕI,µ (r) + R dbR† ϕR,µ (r); the creationPP
field operators are written as ψbµ† (r) = I b
bI† ϕI,µ (r) + R dbR ϕR,µ (r).
[54] P. A. M. Dirac, Proc. R. Soc. Lond. A. 126, 360 (1930).
[55] K. G. Dyall and K. Fægri Jr, Introduction to relativistic quantum chemistry (Oxford University Press, 2007).
[56] We exploit the fact that any unitary matrix can be written in terms of an auxiliary antihermitian matrix κ.
[57] J. Valatin, Il Nuovo Cimento (1955-1965) 7, 843 (1958).
[58] N. Bogoljubov, Il Nuovo Cimento (1955-1965) 7, 794 (1958).
[59] Also known as the single Slater determinant wavefunction.
[60] T. Helgaker, P. Jorgensen, and J. Olsen, Molecular Electronic Structure Theory (Wiley, Chichester, 2000).
[61] L. Visscher, P. Aerts, O. Visser, and W. Nieuwpoort, Int. J. Quant. Chem. 40, 131 (1991).
[62] V. Shabaev, I. Tupitsyn, V. Yerokhin, G. Plunien, and G. Soff, Physical review letters 93, 130405 (2004).
[63] K. G. Dyall, Chem. Phys. 395, 35 (2012).
b and the interacting Hamiltonian operator (H)b commute, which is only[64] Assuming that the opposite charge operator (Q)valid for an appropriate definition of the fermion-fermion interaction and with the adequate definition for the normal
ordering.
[65] Could also be call as the natural spinor representation.
[66] M. Levy, Proc. Natl. Acad. Sci. USA 76, 6062 (1979).
[67] E. H. Lieb, Int. J. Quant. Chem. 24, 243 (1983).
[68] S. M. Valone, Physical Review B 44, 1509 (1991).
[69] We currently have no easy characterization of the set of N -representable 1-RDMs DQ . We at least know that the traceof these 1-RDMs should be equal to the charge Q. Other necessary conditions will be investigated in future work.
[70] Whether we have a true minimum, is an open question, but we suspect it to be the case like in the nonrelativistic limit [67].
[71] Recall that this vacuum state is given using some orthonormal basis ψA (r) that is the initial guess. Also notice that theoperators bbI† , dbR† , b
bI , and dbR are the initial ones.
[72] E. H. Lieb and H. Siedentop, Commun. Math. Phys. 213, 673 (2000).
[73] C. Hainzl, M. Lewin, and É. Séré, Commun. Math. Phys. 257, 515 (2005).
[74] C. Hainzl, M. Lewin, and E. Séré, J. Phys. A Math. Theor. 38, 4483 (2005).
[75] C. Hainzl, M. Lewin, and J. P. Solovej, Comm. Pure Appl. Math. 60, 546 (2007).
[76] C. Hainzl, M. Lewin, E. Séré, and J. P. Solovej, Phys. Rev. A 76, 052104 (2007).
[77] M. Lewin, in Mathematical Results In Quantum Physics: (With DVD-ROM) (World Scientific, 2011) pp. 45–59.
[78] P. Gravejat, M. Lewin, and É. Séré, Commun. Math. Phys. 306, 1 (2011).
[79] C. Brouder, EPJ direct 4, 1 (2002).
[80] J. Sucher, Phys. Rev. A 22, 348 (1980).
[81] J. H. Hubbell, Radiat. Phys. Chem. 75, 614 (2006).
[82] Let us mention that NS are not populated within the np approximation, but they cannot be neglected when optimizing
the PS (i.e. NS play a similar role as the secondary space in nonrelativistic complete active space calculations).
[83] C. Garrod and J. K. Percus, J. Math. Phys. 5, 1756 (1964).
[84] J. M. Herbert and J. E. Harriman, J. Chem. Phys. 118, 10835 (2003).
[85] L. Visscher, O. Visser, P. J. Aerts, H. Merenga, and W. Nieuwpoort, Comput. Phys. Commun. 81, 120 (1994).
[86] J. A. Gaunt, Proc. R. Soc. A 112, 513 (1929).
[87] P.-O. Löwdin and H. Shull, Phys. Rev. 101, 1730 (1956).
[88] M. Piris, Int. J. Quant. Chem. 113, 620 (2013).
[89] M. Piris, J. Chem. Phys. 141, 044107 (2014).
[90] M. Piris, Phys. Rev. A 98, 022504 (2018).
[91] M. Piris, J. Math. Chem. 25, 47 (1999).
[92] Let us comment that the formation of pairs of degenerated spinor is not unique, but we have adopted the one given by
Kramers’ pairs in this work.
[93] The same kind of partitioning of the spinor states can be applied to the NS.
s1 s2[94] The uppercase indices that enter in the definition of x
bIK,JLare readily used in the creation and annihilation operators,
but they are used as (lowercase) unbarred indices in Tij operator. Actually, the Tij is only defined with unbarred indices.
[95] M. Piris, Phys. Rev. Lett. 119, 063002 (2017).

31[96] M. Rodrı́guez-Mayorga, I. Mitxelena, F. Bruneval, and M. Piris, J. Chem. Theory Comput. 17, 7562 (2021).
[97] L. M. Mentel, R. van Meer, O. V. Gritsenko, and E. J. Baerends, J. Chem. Phys. 140, 214105 (2014).
[98] M. Rodrı́guez-Mayorga, E. Ramos-Cordoba, M. Via-Nadal, M. Piris, and E. Matito, Phys. Chem. Chem. Phys. 19, 24029(2017).
[99] S. Sun, T. F. Stetina, T. Zhang, H. Hu, E. F. Valeev, Q. Sun, and X. Li, J. Chem. Theory Comput. 17, 3388 (2021).
[100] P. Hafner, J. Phys. B 13, 3297 (1980).
[101] A. M. K. Müller, Phys. Lett. 105A, 446 (1984).
[102] M. A. Buijse, Thesis: Electron Correlation. Fermi and Coulomb holes, dynamical and nondynamical correlation, Ph.D.
thesis, Vrije Universiteit, Amsterdam, The Netherlands (1991).
[103] M. A. Buijse and E. J. Baerends, Mol. Phys. 100, 401 (2002).
[104] J. Cioslowski and K. Pernal, J. Chem. Phys. 111, 3396 (1999).
[105] S. Sharma, J. K. Dewhurst, N. N. Lathiotakis, and E. K. U. Gross, Phys. Rev. B 78, 201103 (2008).
[106] J. Cioslowski and K. Pernal, Phys. Rev. A 61, 034503 (2000).
[107] S. Goedecker and C. J. Umrigar, Phys. Rev. Lett. 81, 866 (1998).
[108] M. A. L. Marques and N. N. Lathiotakis, Phys. Rev. A 77, 032509 (2008).
[109] E. H. Lieb and S. Oxford, Int. J. Quant. Chem. 19, 427 (1981).
[110] E. Kamil, R. Schade, T. Pruschke, and P. E. Blöchl, Phys. Rev. B 93, 085141 (2016).
[111] P.-O. Löwdin, Phys. Rev. 97, 1474 (1955).
[112] M. Piris, J. Matxain, and X. Lopez, J. Chem. Phys. 139, 234109 (2013).
[113] I. Mitxelena, M. Rodrı́guez-Mayorga, and M. Piris, Eur. Phys. J. B 91, 109 (2018).
[114] A variant proposed for PNOF7 functional that is used to define NOF-MP2 method [90].
[115] The relativistic version of PNOF6 functional is not included because one of us proved that beyond the perfect-pairing
approach this functional does not obey the sum rule [98].
[116] M. Piris, Phys. Rev. Lett. 127, 233001 (2021).
[117] H. Shull and P.-O. Löwdin, J. Chem. Phys. 30, 617 (1959).
[118] K. Pernal, Comput. Theor. Chem. 1003, 127 (2013).
[119] A. J. Coleman, Rev. Mod. Phys. 35, 668 (1963).
[120] D. W. Smith, Phys. Rev. 147, 896 (1966).
[121] R. Borland and K. Dennis, J. Phys. B 5, 7 (1972).
[122] A. Coleman, Int. J. Quant. Chem. 13, 67 (1978).
[123] A. A. Klyachko, in Journal of Physics: Conference Series, Vol. 36 (IOP Publishing, 2006) p. 72.
[124] M. Altunbulak and A. Klyachko, Commun. Math. Phys. 282, 287 (2008).
[125] D. A. Mazziotti, in Many-Electron Densities and Reduced Density Matrices (Springer, 2000) pp. 139–163.
[126] C. Schilling, M. Altunbulak, S. Knecht, A. Lopes, J. D. Whitfield, M. Christandl, D. Gross, and M. Reiher, Phys. Rev.
A 97, 052503 (2018).
[127] I. Theophilou, N. N. Lathiotakis, M. A. L. Marques, and N. Helbig, J. Chem. Phys. 142, 154108 (2015).
[128] R. Chakraborty and D. A. Mazziotti, Phys. Rev. A 89, 042505 (2014).
[129] R. Chakraborty and D. A. Mazziotti, Int. J. Quant. Chem. 116, 784 (2016).
[130] D. A. Mazziotti, Phys. Rev. A 94, 032516 (2016).
[131] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
[132] R. Pino, O. Bokanowski, E. V. Ludeña, and R. L. Boada, Theor. Chim. Acta (Berlin) 118, 557 (2007).
[133] M. Penz and R. van Leeuwen, arXiv preprint arXiv:2106.15370 (2021).
[134] Notice that the λīj̄,k̄l̄ terms are obtained by replacing the unbarred indices by barred ones.
[135] For i = j, the elements are 0 by definition and they fulfill the condition.
[136] M. Piris, X. Lopez, F. Ruipérez, J. M. Matxain, and J. M. Ugalde, J. Chem. Phys. 134, 164102 (2011).
[137] Recall that the ‘spin’ information is contained in the 4-component spinors.
[138] B. Carlson and J. M. Keller, Phys. Rev. 121, 659 (1961).
[139] E. Schmidt, Math. Ann. 63, 433 (1907).
b in the many electron case leads to the complex conjugate of the CI coefficients and transforms all spinors in the tensor[140] K
product [55, 141].
[141] H. Jørgen A. Jensen, K. G. Dyall, T. Saue, and K. Fægri Jr, J. Chem. Phys. 104, 4083 (1996).
[142] This is normally the case in the nonrelativistic limit [55].
ei,1 (r) = 0) we would arrive to similar conclusions.
[143] Retaining only the imaginary part of the scalar spinors (setting ℜϕei,1 (r) = ℑϕei,1 (r) to remove the triplet
Nevertheless, to retain the real and the imaginary parts we need to impose ℜϕcontribution.
[144] In the a = 1 case the spin functions are those of the Sbx operator. On the other hand, taking a = i we arrive to spinfunctions of the Sby operator.

32

FIG. 2. Graphical abstract.