# Relativistic reduced density matrix functional theory.

M. Rodríguez-Mayorga<EMAIL_ADDRESS>Department of Theoretical Chemistry,
Vrije Universiteit Amsterdam, De Boelelaan 1083, 1081 HV Amsterdam, The
Netherlands.  K.J.H. Giesbertz  Department of Theoretical Chemistry, Vrije
Universiteit Amsterdam, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands.
L. Visscher  Department of Theoretical Chemistry, Vrije Universiteit
Amsterdam, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands.

###### Abstract

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 first-order 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 the functional 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 versions of
some common RDMFT approximations that are used in the nonrelativistic context
and discuss their properties.

###### Contents

  1. I Introduction
  2. II ReRDMFT including electron-positron pair creation and annihilation processes.
    1. II.1 The free particle Dirac equation and quantization of the Dirac field
    2. II.2 The Hamiltonian operator including an external nonlocal potential, charge sectors in Fock space, and Bogoliubov transformations
    3. II.3 Vacuum polarization
    4. II.4 The interacting particle Hamiltonian and the configuration-interaction vacuum
    5. II.5 The 1-RDM and its natural orbital representation
    6. II.6 Introducing ReRDMFT through the constrained-search formalism
  3. III The no-pair approximation and ReRDMFT
    1. III.1 Vacuum polarization effects at fermion-fermion interacting case
    2. III.2 The no-pair approximation including vacuum polarization ReRDFMT approach
    3. III.3 The no-pair approximation ReRDFMT approach
    4. III.4 The $\langle\Psi_{+}|\widehat{\widetilde{W}}|\Psi_{+}\rangle\;$ term as an explicit functional of the second-order reduced density matrix.
    5. III.5 np-ReNOFT functional approximations by imposing Kramers’ symmetry
      1. III.5.1 The Dirac–Hartree–Fock approximation
      2. III.5.2 $f(n_{I},n_{J})$-functional approximations
      3. III.5.3 Approximations built imposing $N$-representability conditions
    6. III.6 Properties of rel-GNOF/PNOF$x$ functionals.
  4. IV Closing remarks
  5. A Commutation relation between the Hamiltonian and the opposite charge operator
  6. B Extending Gilbert’s theorem to the relativistic domain.
  7. C The npvp energy functional
  8. D The $N$-representability conditions within the np approximation.
  9. E The relativistic Fixed-Phases functional.

##  I Introduction

When relativistic effects play a role in the quantum world, the Schrö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, 2, 3, 3, 4, 5] as electrons reach high velocities due to
their strongly attractive nuclear potentials. When such elements 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, 7, 8, 9, 10, 11, 12, 13,
14]. While methods like complete active space including second-order
perturbation theory (CASPT2) [15, 16, 17, 18], DMRG [19, 20] and multi-
reference coupled cluster theory [21, 22, 23, 24, 25, 26] can be employed they
all exhibit a high computational scaling with system size 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, 30, 31, 32, 33, 34, 35].
Indeed, the popular working horse of physicists and chemists (i.e. the use of
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 exploited and 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 et al.
(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-pair approximation [48, 49].

This work is organized as follows: 1) we introduce the most general ReRDMFT
approach, where creation and annihilation of electron-positron pairs is
allowed. To that end, we first discuss the single-particle problem subject to
a nonlocal external potential; then, we present its many-particle
generalization and introduce the Fock space. Next, we analyze 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-positron pairs 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 where these 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. Finally, we discuss
some of the properties of the functional approximations. Before proceeding,
let us present Table 1 where we have collected the index conventions that will
be employed in this work.

Table 1: Indices used throughout this work. Indices | labeling
---|---
$I$,$J$,$K$,$L$ | positive energy spinors $\left(\textrm{PS}\right)$
$R$,$S$ | negative energy spinors $\left(\textrm{NS}\right)$
$A$,$B$ | all spinors ($\text{NS}\cup\text{PS}$)
$i$,$j$,$k$,$l$ | half of the PS
(not related by Kramers’ symmetry)
$\bar{i}$,$\bar{j}$,$\bar{k}$,$\bar{l}$ | half of the PS
(Kramers partners of the unbarred spinors)
$\mu$, $\nu$, $\tau$, $\eta$ | scalar orbitals (components of a 4-spinor)
$p$,$q$ | electron pairs

##  II ReRDMFT including electron-positron pair creation and annihilation
processes.

###  II.1 The free particle Dirac equation and quantization of the Dirac field

Let us start by defining the time-independent free particle Dirac equation

$\displaystyle\widehat{{\bf T}}_{D}({\bf r}){\bm{\psi}}_{A}({\bf r})$
$\displaystyle=\left[-\textrm{i}c({\bm{\alpha}}\cdot{\bm{\nabla}}_{{\bf
r}})+c^{2}m{\bm{\beta}}\right]{\bm{\psi}}_{A}({\bf r})$
$\displaystyle=E_{A}{\bm{\psi}}_{A}({\bf r}),$ (1)

where $\widehat{{\bf T}}_{D}({\bf r})$ is the usual first-quantized $4\times
4$ Dirac kinetic + rest mass operator, $\textrm{i}=\sqrt{-1}$, $c=137.036$
a.u. is the speed of light, $m=1$ a.u. is the electron mass,

$\displaystyle{\bm{\alpha}}$
$\displaystyle=({\bm{\alpha}}_{x},{\bm{\alpha}}_{y},{\bm{\alpha}}_{z})$
$\displaystyle=\left({\begin{pmatrix}{\bm{0}}_{2}&{\bm{\sigma}}_{x}\\\
{\bm{\sigma}}_{x}&{\bm{0}}_{2}\end{pmatrix}},{\begin{pmatrix}{\bm{0}}_{2}&{\bm{\sigma}}_{y}\\\
{\bm{\sigma}}_{y}&{\bm{0}}_{2}\end{pmatrix}},{\begin{pmatrix}{\bm{0}}_{2}&{\bm{\sigma}}_{z}\\\
{\bm{\sigma}}_{z}&{\bm{0}}_{2}\end{pmatrix}}\right),$ (2)

${\bm{0}}_{2}$ is the 2$\times$2 null matrix,

$\displaystyle{\bm{\sigma}}_{x}$ $\displaystyle={\begin{pmatrix}0&1\\\
1&0\end{pmatrix}},$ $\displaystyle{\bm{\sigma}}_{y}$
$\displaystyle={\begin{pmatrix}0&-\textrm{i}\\\ \textrm{i}&0\end{pmatrix}},$
$\displaystyle{\bm{\sigma}}_{z}$ $\displaystyle={\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}};$ (3)
${\bm{\beta}}=\begin{pmatrix}\mathbb{I}_{2}&{\bm{0}}_{2}\\\
{\bm{0}}_{2}&-\mathbb{I}_{2}\end{pmatrix},$ (4)

$\mathbb{I}_{2}$ is the 2$\times$2 unit matrix. Solutions of the free Dirac
equation (${\bm{\psi}}_{A}$) are 4-component-spinor orbitals

${\bm{\psi}}_{A}({\bf r})=\begin{pmatrix}\phi_{A,1}({\bf r})\\\
\phi_{A,2}({\bf r})\\\ \phi_{A,3}({\bf r})\\\ \phi_{A,4}({\bf
r})\end{pmatrix},$ (5)

whose conjugate-transpose form reads

$\displaystyle{\bm{\psi}}^{\dagger}_{A}({\bf r})$
$\displaystyle=\begin{pmatrix}\phi^{*}_{A,1}({\bf r})&\phi_{A,2}^{*}({\bf
r})&\phi_{A,3}^{*}({\bf r})&\phi_{A,4}^{*}({\bf r})\end{pmatrix}.$ (6)

Let us remark that the scalar functions $(\phi)$ represent spatial orbitals
because the spin is accounted by their arrangement in 4-component spinors.
These solutions can be partitioned into a set of positive-energy 4-component
spinors ($E_{A}>0$) and a set of negative-energy ones ($E_{A}<0$), i.e.
$\\{{\bm{\psi}}_{A}\\}=\\{{\bm{\psi}}_{R}\\}\cup\\{{\bm{\psi}}_{I}\\}$. From
now on, we will designate the “positive energy spinors” (“negative energy
spinors”) as PS (NS). Hence, the Dirac field is quantized as

${\widehat{\bm{\psi}}^{\vphantom{\dagger}}}({\bf
r})=\sum_{\mathclap{A}}{\widehat{a}^{\vphantom{\dagger}}}_{A}{\bm{\psi}}_{A}({\bf
r})=\sum_{I}{\widehat{b}^{\vphantom{\dagger}}}_{I}{\bm{\psi}}_{I}({\bf
r})+\sum_{R}{\widehat{d}^{\;\dagger}}_{R}{\bm{\psi}}_{R}({\bf r}),$ (7)

where the sum has been decomposed into contributions involving
electron(positron) annihilation(creation) operators. These operators obey the
usual anticommutation relations

$\\{{\widehat{a}^{\vphantom{\dagger}}}_{A},{\widehat{a}^{\;\dagger}}_{B}\\}=\delta_{AB}\;\;\textrm{and}\;\;\\{{\widehat{a}^{\;\dagger}}_{A},{\widehat{a}^{\;\dagger}}_{B}\\}=\\{{\widehat{a}^{\vphantom{\dagger}}}_{A},{\widehat{a}^{\vphantom{\dagger}}}_{B}\\}=0.$
(8)

###  II.2 The Hamiltonian operator including an external nonlocal potential,
charge sectors in Fock space, and Bogoliubov transformations

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

${\bf v}^{\textrm{nl}}_{\textrm{ext}}({\bf r},{\bf
r}^{\prime})=\begin{pmatrix}v^{\textrm{nl}}_{\textrm{ext},1,1}({\bf r},{\bf
r}^{\prime})&v^{\textrm{nl}}_{\textrm{ext},1,2}({\bf r},{\bf
r}^{\prime})&v^{\textrm{nl}}_{\textrm{ext},1,3}({\bf r},{\bf
r}^{\prime})&v^{\textrm{nl}}_{\textrm{ext},1,4}({\bf r},{\bf r}^{\prime})\\\
v^{\textrm{nl}}_{\textrm{ext},2,1}({\bf r},{\bf
r}^{\prime})&v^{\textrm{nl}}_{\textrm{ext},2,2}({\bf r},{\bf
r}^{\prime})&v^{\textrm{nl}}_{\textrm{ext},2,3}({\bf r},{\bf
r}^{\prime})&v^{\textrm{nl}}_{\textrm{ext},2,4}({\bf r},{\bf r}^{\prime})\\\
v^{\textrm{nl}}_{\textrm{ext},3,1}({\bf r},{\bf
r}^{\prime})&v^{\textrm{nl}}_{\textrm{ext},3,2}({\bf r},{\bf
r}^{\prime})&v^{\textrm{nl}}_{\textrm{ext},3,3}({\bf r},{\bf
r}^{\prime})&v^{\textrm{nl}}_{\textrm{ext},3,4}({\bf r},{\bf r}^{\prime})\\\
v^{\textrm{nl}}_{\textrm{ext},4,1}({\bf r},{\bf
r}^{\prime})&v^{\textrm{nl}}_{\textrm{ext},4,2}({\bf r},{\bf
r}^{\prime})&v^{\textrm{nl}}_{\textrm{ext},4,3}({\bf r},{\bf
r}^{\prime})&v^{\textrm{nl}}_{\textrm{ext},4,4}({\bf r},{\bf
r}^{\prime})\end{pmatrix},$ (9)

with the matrix elements obeying the constraint
$v^{\text{nl}}_{\textrm{ext},\mu,\nu}({\bf r},{\bf
r}^{\prime})=\left[v^{\text{nl}}_{\textrm{ext},\nu,\mu}({\bf r}^{\prime},{\bf
r})\right]^{*}$ to ensure the hermiticity of the operator.

The bound particle Dirac equation that describes the states of a single
particle (electron or positron) subject to a nonlocal external potential as in
(9) reads

$\widehat{\bf T}_{D}({\bf r}){\bm{\psi}}_{A}({\bf r})+\int d{\bf
r}^{\prime}{\bf v}_{\textrm{ext}}^{\textrm{nl}}({\bf r},{\bf
r}^{\prime}){\bm{\psi}}_{A}({\bf r}^{\prime})=E_{A}{\bm{\psi}}_{A}({\bf r}).$
(10)

The states ${\bm{\psi}}$ depend on the choice of nonlocal external potential
and form an orthonormal set ($\int d{\bf r}{\bm{\psi}}^{\dagger}_{A}({\bf
r}){\bm{\psi}}_{B}({\bf r})=\delta_{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 as

$\displaystyle\widehat{H}_{0}^{v}$
$\displaystyle=\widehat{T}_{D}+\widehat{V}^{\textrm{nl}}_{\textrm{ext}}$
$\displaystyle=\int d{\bf r}d{\bf r}^{\prime}\delta({\bf r}-{\bf
r}^{\prime})\textrm{Tr}\left[\widehat{{\bf T}}_{D}({\bf r})\widehat{{\bf
n}}_{1}({\bf r},{\bf r}^{\prime})\right]+\int d{\bf r}d{\bf
r}^{\prime}\textrm{Tr}\left[{\bf v}_{\textrm{ext}}^{\textrm{nl}}({\bf
r}^{\prime},{\bf r})\widehat{{\bf n}}_{1}({\bf r},{\bf r}^{\prime})\right]$
(11) $\displaystyle=\sum_{\mu,\tau}\int d{\bf r}d{\bf
r}^{\prime}\left[\delta({\bf r}-{\bf r}^{\prime})T_{D,\mu,\tau}({\bf
r})+{v}_{\textrm{ext},\mu,\tau}^{\textrm{nl}}({\bf r}^{\prime},{\bf
r})\right]\widehat{n}_{1,\tau,\mu}({\bf r},{\bf r}^{\prime}),$

where we have introduced the one-particle density matrix operator whose
elements read as

$\displaystyle\widehat{n}_{1,\tau,\mu}({\bf r},{\bf r}^{\prime})$
$\displaystyle=\mathcal{N}\left[{\widehat{\psi}^{\;\dagger}}_{\mu}({\bf
r}^{\prime}){\widehat{\psi}^{\vphantom{\dagger}}}_{\tau}({\bf r})\right]$
$\displaystyle=\sum_{I,J}{\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{J}\phi^{*}_{I,\mu}({\bf
r}^{\prime})\phi_{J,\tau}({\bf
r})+\sum_{I}\sum_{S}{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{S}\phi^{*}_{I,\mu}({\bf
r}^{\prime})\phi_{S,\tau}({\bf r})$ (12)
$\displaystyle+{}\sum_{R}\sum_{J}{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{J}\phi^{*}_{R,\mu}({\bf
r}^{\prime})\phi_{J,\tau}({\bf
r})-\sum_{R,S}{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{R}\phi^{*}_{R,\mu}({\bf
r}^{\prime})\phi_{S,\tau}({\bf r})$

that is defined using creation and annihilation Dirac field operators with
normal ordering [51, 52] (using Wick’s theorem)
$\mathcal{N}\left[\dotsb\right]$ of the elementary operators 111The spinor-
components annihilation field operators are defined as
${\widehat{\psi}^{\vphantom{\dagger}}}_{\mu}({\bf
r})=\sum_{I}{\widehat{b}^{\vphantom{\dagger}}}_{I}\phi_{I,\mu}({\bf
r})+\sum_{R}{\widehat{d}^{\;\dagger}}_{R}\phi_{R,\mu}({\bf r})$; the creation
field operators are written as ${\widehat{\psi}^{\;\dagger}}_{\mu}({\bf
r})=\sum_{I}{\widehat{b}^{\;\dagger}}_{I}\phi_{I,\mu}({\bf
r})+\sum_{R}{\widehat{d}^{\vphantom{\dagger}}}_{R}\phi_{R,\mu}({\bf r})$.
${\widehat{b}^{\;\dagger}}_{I}$, ${\widehat{d}^{\;\dagger}}_{R}$,
${\widehat{b}^{\vphantom{\dagger}}}_{I}$, and
${\widehat{d}^{\vphantom{\dagger}}}_{R}$ w.r.t. the effective vacuum state
($|0_{v}\rangle$). The formal definition of the effective vacuum state
requires a detailed discussion that we briefly summarize in the following
lines. First we note that any Hamiltonian leads to a set of PS and NS. One
could then define a bare vacuum state in which all spinors are unoccupied.
This definition of the vacuum is problematic, however. Since NS are much lower
in energy 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
(by electrons) NS are redefined as empty positronic spinors. Thus, the
effective vacuum state, which sets the zero of the energy scale, corresponds
to the state where all PS do not contain electrons and all NS do not contain
positrons. In the particular case when ${\bf
v}_{\textrm{ext}}^{\textrm{nl}}={\bf 0}_{4\times 4}$
($\widehat{H}_{0}^{v}=\widehat{H}_{0}^{0}\equiv\widehat{T}_{D}$), we refer to
its effective vacuum state as $|0_{0}\rangle$ and to its spinor basis as
${\bm{\psi}}^{0}_{A}({\bf r})$.

Let us assume that a given $\widehat{H}_{0}^{v}$ (with ${\bf
v}_{\textrm{ext}}^{\textrm{nl}}\neq{\bf 0}_{4\times 4}$) is built initially in
the free particle ${\bm{\psi}}^{0}_{A}({\bf r})$ basis, with the normal
ordering taken w.r.t. free particle vacuum $|0_{0}\rangle$. In this basis
$\widehat{H}_{0}^{v}$ is not diagonal and it does not commute with the
electron (positron) _number operators_
$\widehat{N}_{e}=\sum_{I}{\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{I}$($\widehat{N}_{p}=\sum_{R}{\widehat{d}^{\;\dagger}}_{R}{\widehat{d}^{\vphantom{\dagger}}}_{R}$).
Moreover, $\widehat{H}_{0}^{v}$ does not commute with the total number of
particles operator $\widehat{N}=\widehat{N}_{e}+\widehat{N}_{p}$. A useful
operator that does commute with $\widehat{H}_{0}^{v}$ is, however, the
_opposite charge operator_ [55, 49] that is defined as

$\widehat{Q}=\widehat{N}_{e}-\widehat{N}_{p}=\sum_{I}{\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{I}-\sum_{R}{\widehat{d}^{\;\dagger}}_{R}{\widehat{d}^{\vphantom{\dagger}}}_{R}$
(13)

(see Appendix A for more details). We then consider using the eigensolutions
${\bm{\psi}}^{v}_{A}$ that provide the diagonal representation
$\widehat{\widetilde{H^{v}_{0}}}$ with the tilde symbol indicating that the
Hamiltonian is transformed with respect to its original representation and has
normal ordering w.r.t. its own vacuum $|0_{v}\rangle$. For consistency we will
also indicate this vacuum as $|\widetilde{0}_{v}\rangle$ with the tilde symbol
indicating that these quantities are transformed with respect to their
original representations. In the current case
$|\widetilde{0}_{v}\rangle=|0_{v}\rangle$, but this will change when two-
particle interactions are considered below. In the
$\widetilde{\bm{\psi}}_{A}({\bf r})={\bm{\psi}}^{v}_{A}({\bf r})$ basis the
Hamiltonian commutes with the number of electrons (positrons) operators, the
number of particles operator, and the opposite charge operator. Due to the
diagonal representation, the pair creation and annihilation processes vanish;
thence, the Hamiltonian conserves the particle number. Since both basis are
orthonormal, they are related by a unitary transformation ${\bf
V}=\textrm{e}^{\bm{\kappa}}$ 222We exploit the fact that any unitary matrix
can be written in terms of an auxiliary antihermitian matrix ${\bm{\kappa}}$,
where the antihermitian matrix ${\bm{\kappa}}$ (formed by the parameters
$\kappa_{AB}\in\mathbb{C}$) can be chosen to contain all its diagonal entries
equal to zero to avoid redundant phase-shifts of the spinors. Therefore, the
relationship between the two bases reads as

$\displaystyle\widetilde{\bm{\psi}}_{A}({\bf r})$
$\displaystyle=\sum_{B}{\bm{\psi}}^{0}_{B}({\bf r})V_{BA}.$ (14)

The $\widehat{\widetilde{H^{v}_{0}}}$ operator is written in terms of the
transformed operators (${\widehat{\widetilde{b}}^{\;\dagger}}_{I}$,
${\widehat{\widetilde{d}}^{\;\dagger}}_{R}$,
${\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{I}$, and
${\widehat{\widetilde{d}}^{\vphantom{\dagger}}}_{R}$) that are given by

$\displaystyle{\widehat{\widetilde{a}}^{\vphantom{\dagger}}}_{A}$
$\displaystyle=\textrm{e}^{\widehat{\kappa}}{\widehat{a}^{\vphantom{\dagger}}}_{A}\textrm{e}^{-\widehat{\kappa}}=\sum_{\mathclap{B}}{\widehat{a}^{\vphantom{\dagger}}}_{B}{V}_{BA}^{*},$
(15)

where the spinor rotation operator ($e^{\widehat{\kappa}}$) is expressed as
the exponential of an antihermitian operator $\widehat{\kappa}$ with

$\displaystyle\widehat{\kappa}={}$
$\displaystyle\sum_{\mathclap{A,B}}\kappa_{AB}{\widehat{a}^{\;\dagger}}_{A}{\widehat{a}^{\vphantom{\dagger}}}_{B}$
$\displaystyle={}$
$\displaystyle\sum_{\mathclap{I,J}}\kappa_{IJ}{\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{J}+\sum_{I}\sum_{S}\kappa_{IS}{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{S}$
(16)
$\displaystyle{}+\sum_{R}\sum_{J}\kappa_{RJ}{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{J}+\sum_{\mathclap{R,S}}\kappa_{RS}{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{d}^{\;\dagger}}_{S}.$

Actually, the spinor rotation operator $\textrm{e}^{\widehat{\kappa}}$
corresponds to a Bogoliubov transformation mixing electron annihilation and
positron creation operators [57, 58], producing modified electron and positron
creation and annihilation operators that 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

$\mathcal{F}=\underset{(N_{e},N_{p})=(0,0)}{\stackrel{{\scriptstyle(\infty,\infty)}}{{\bigoplus}}}\mathclap{\mathcal{H}^{(N_{e},N_{p})}}\qquad\;,$
(17)

where $\bigoplus$ designates the direct sum. Alternatively, the Fock space can
also be decomposed into charge ($Q=N_{e}-N_{p}$) sectors

$\mathcal{F}=\underset{Q=-\infty}{\stackrel{{\scriptstyle\infty}}{{\bigoplus}}}\mathclap{\mathcal{H}_{Q}}\;\;\;.$
(18)

For a given Hamiltonian $\widehat{H}_{0}^{v}$, the wavefunction of a system
containing a fixed number of electron ($N_{e}$) and positrons ($N_{p}$) can be
written as a single determinant (SD) 333Also known as the single Slater
determinant wavefunction.. The SD wavefunction reads as

$\lvert\Phi\rangle={\widehat{b}^{\;\dagger}}_{I_{1}}\dotsb{\widehat{b}^{\;\dagger}}_{I_{N_{e}}}{\widehat{d}^{\;\dagger}}_{R_{1}}\dotsb{\widehat{d}^{\;\dagger}}_{R_{N_{p}}}\lvert
0_{0}\rangle,$ (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 $\widehat{V}$ as

$\displaystyle\lvert\widetilde{\Phi}\rangle$
$\displaystyle=\widehat{V}\lvert\Phi\rangle=\textrm{e}^{\widehat{\kappa}}\lvert\Phi\rangle=\textrm{e}^{\widehat{\kappa}}\;{\widehat{b}^{\;\dagger}}_{I_{1}}\dotsb{\widehat{b}^{\;\dagger}}_{I_{N_{e}}}{\widehat{d}^{\;\dagger}}_{R_{1}}\dotsb{\widehat{d}^{\;\dagger}}_{R_{N_{p}}}\lvert
0_{0}\rangle$
$\displaystyle=\textrm{e}^{\widehat{\kappa}}\;{\widehat{b}^{\;\dagger}}_{I_{1}}\textrm{e}^{-\widehat{\kappa}}\dotsb\textrm{e}^{\widehat{\kappa}}{\widehat{b}^{\;\dagger}}_{I_{N_{e}}}\textrm{e}^{-\widehat{\kappa}}\textrm{e}^{\widehat{\kappa}}{\widehat{d}^{\;\dagger}}_{R_{1}}\textrm{e}^{-\widehat{\kappa}}\dotsb\textrm{e}^{\widehat{\kappa}}{\widehat{d}^{\;\dagger}}_{R_{N_{p}}}\textrm{e}^{-\widehat{\kappa}}\textrm{e}^{\widehat{\kappa}}\lvert
0_{0}\rangle$
$\displaystyle={\widehat{\widetilde{b}}^{\;\dagger}}_{I_{1}}\dotsb{\widehat{\widetilde{b}}^{\;\dagger}}_{I_{N_{e}}}{\widehat{\widetilde{d}}^{\;\dagger}}_{R_{1}}\dotsb{\widehat{\widetilde{d}}^{\;\dagger}}_{R_{N_{p}}}|\widetilde{0}_{v}\rangle,$
(20)

with $|\widetilde{0}_{v}\rangle=e^{\widehat{\kappa}}\lvert 0_{0}\rangle$. Let
us also rewrite the annihilation Dirac field operators in the
$\widetilde{\bm{\psi}}_{A}({\bf r})$ basis as

${\widehat{\bm{\psi}}^{\vphantom{\dagger}}}({\bf
r})=\sum_{\mathclap{A}}{\widehat{\widetilde{a}}^{\vphantom{\dagger}}}_{A}\widetilde{\bm{\psi}}_{A}({\bf
r})=\sum_{\mathclap{I}}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{I}\widetilde{\bm{\psi}}_{I}({\bf
r})+\sum_{\mathclap{R}}{\widehat{\widetilde{d}}^{\;\dagger}}_{R}\widetilde{\bm{\psi}}_{R}({\bf
r}),$ (21)

a similar expression can be written for the creation Dirac field operator
($\widehat{\bm{\psi}}^{\dagger}({\bf 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 given nonlocal external
potential. Finally, taking the Taylor expansion of
$\textrm{e}^{\widehat{\kappa}}$, it is easy to recognize that basis
$\widetilde{\bm{\psi}}_{A}({\bf r})$ in which the Hamiltonian is diagonal is
given by the set of $\\{\kappa_{AB}\\}$ parameters that lead to a stationary
point [60] (i.e. a minimum, a maximum, or a saddle point) of the energy w.r.t.
variations of these parameters.

###  II.3 Vacuum polarization

From now on, we will assume that an initially given orthonormal basis
${\bm{\psi}}_{A}({\bf r})$ is employed to express the Hamiltonian operator,
this could be the ${\bm{\psi}}^{0}_{A}({\bf r})$ basis, but any orthonormal
basis can in principle be used as long as it is possible to distinguish
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, 62, 63] between the bases for the
large and small components of the spinors. The initial effective vacuum state
is now defined as $\lvert 0\rangle$ 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 that the Hamiltonian is non-diagonal in this
basis. In the previous section we introduced the rotation operator
($\textrm{e}^{\widehat{\kappa}}$) and defined the transformed effective vacuum

$\lvert\widetilde{0}_{v}\rangle=e^{\widehat{\kappa}}\lvert 0\rangle$ (22)

that provides a new reference for normal ordering the transformed operators.
As we previously mentioned, the initially defined $\widehat{H}_{0}^{v}$
operator can be rewritten in another basis ($\widehat{\widetilde{H^{v}_{0}}}$)
employing the operator $\widehat{\widetilde{n}}_{1,\tau,\mu}({\bf r},{\bf
r}^{\prime})$ (i.e. expressing the one-body density matrix operator in terms
of the elementary operators ${\widehat{\widetilde{b}}^{\;\dagger}}_{I}$,
${\widehat{\widetilde{d}}^{\;\dagger}}_{R}$,
${\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{I}$, and
${\widehat{\widetilde{d}}^{\vphantom{\dagger}}}_{R}$). Recalling that the
normal ordering is introduced to guarantee that the energy of the effective
vacuum is zero, i.e.,

$\langle 0\lvert\widehat{H}_{0}^{v}\lvert
0\rangle=\langle\widetilde{0}_{v}\lvert\widehat{\widetilde{H^{v}_{0}}}\lvert\widetilde{0}_{v}\rangle=0,$
(23)

when a Hamiltonian is not evaluated with its reference effective vacuum state
it leads to nonzero contributions
$\langle\widetilde{0}_{v}|\widehat{H}_{0}^{v}|\widetilde{0}_{v}\rangle\neq
0\neq\langle 0|\widehat{\widetilde{H^{v}_{0}}}|0\rangle$. Actually, the
relationship between the noninteracting Hamiltonian operators can simply be
written as

$\widehat{{H}}_{0}^{v}=\widehat{\widetilde{H^{v}_{0}}}+\widetilde{E}_{0}^{0},$
(24)

where

$\widetilde{E}_{0}^{0}=\langle\widetilde{0}_{v}|\widehat{{H}}_{0}^{v}|\widetilde{0}_{v}\rangle=\int
d{\bf r}d{\bf r}^{\prime}\textrm{Tr}\left[\left(\delta({\bf r}-{\bf
r}^{\prime})\widehat{{\bf T}}_{D}({\bf r})+{\bf
v}_{\textrm{ext}}^{\textrm{nl}}({\bf r}^{\prime},{\bf
r})\right)\widetilde{{\bf n}}^{\textrm{vp}}_{1}({\bf r},{\bf
r}^{\prime})\right],$ (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 as

$\widetilde{{\bf n}}^{\textrm{vp}}_{1}({\bf r},{\bf
r}^{\prime})=\sum_{R}\widetilde{\bm{\psi}}_{R}({\bf
r})\widetilde{\bm{\psi}}^{\dagger}_{R}({\bf
r}^{\prime})-\sum_{R}{\bm{\psi}}_{R}({\bf r}){\bm{\psi}}^{\dagger}_{R}({\bf
r}^{\prime})$ (26)

whose components read as

$\displaystyle\widetilde{n}^{\textrm{vp}}_{1,\eta,\nu}({\bf r},{\bf
r}^{\prime})$
$\displaystyle=\langle\widetilde{0}_{v}|\widehat{{n}}_{1,\eta,\nu}({\bf
r},{\bf r}^{\prime})|\widetilde{0}_{v}\rangle$ (27)
$\displaystyle=\sum_{R}\widetilde{\phi}_{R,\eta}({\bf
r})\widetilde{\phi}^{*}_{R,\nu}({\bf r}^{\prime})-\sum_{R}{\phi}_{R,\eta}({\bf
r}){\phi}^{*}_{R,\nu}({\bf r}^{\prime}).$

Before proceeding, let us remark that vp effects are present whenever spinor
rotations mixing the positive and negative 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.

###  II.4 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 take a 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 widely used interaction at present; it improves the description of the
interaction w.r.t. nonrelativistic Coulomb interaction (i.e. $\lvert{\bf
r}_{1}-{\bf r}_{2}\rvert^{-1}$). Thus, the approximate interacting Hamiltonian
reads as

$\displaystyle\widehat{H}$ $\displaystyle=\widehat{H}_{0}^{v}+\widehat{W}$
$\displaystyle=\widehat{H}_{0}^{v}+\frac{1}{2}\int d{\bf r}_{1}d{\bf
r}_{2}\textrm{Tr}\left[{\bf W}({\bf r}_{1},{\bf r}_{2})\widehat{{\bf
n}}_{2}({\bf r}_{1},{\bf r}_{2})\right]$
$\displaystyle=\widehat{H}_{0}^{v}+\frac{1}{2}\sum_{\mu,\nu,\tau,\eta}\int
d{\bf r}_{1}d{\bf r}_{2}W_{\mu,\nu,\tau,\eta}({\bf r}_{1},{\bf
r}_{2})\widehat{n}_{2,\tau,\eta,\mu,\nu}({\bf r}_{1},{\bf r}_{2}),$ (28)

where

$W_{\mu,\nu,\tau,\eta}({\bf r}_{1},{\bf
r}_{2})=\frac{1}{r_{12}}\left[\delta_{\mu,\tau}\delta_{\nu,\eta}-\frac{1}{2}\left({\bm{\alpha}}_{\mu,\tau}\cdot{\bm{\alpha}}_{\nu,\eta}\right)-\frac{\left({\bm{\alpha}}_{\mu,\tau}\cdot{\bf
r}_{12}\right)\left({\bm{\alpha}}_{\nu,\eta}\cdot{\bf
r}_{12}\right)}{2r_{12}^{2}}\right]$ (29)

with ${\bf r}_{12}={\bf r}_{1}-{\bf r}_{2}$ and $r_{12}=\lvert{\bf
r}_{12}\rvert$; and the pair-density matrix operator is defined using creation
and annihilation Dirac field operators

$\widehat{n}_{2,\tau,\eta,\mu,\nu}({\bf r}_{1},{\bf
r}_{2})=\mathcal{N}\left[{\widehat{\psi}^{\;\dagger}}_{\nu}({\bf
r}_{2}){\widehat{\psi}^{\;\dagger}}_{\mu}({\bf
r}_{1}){\widehat{\psi}^{\vphantom{\dagger}}}_{\tau}({\bf
r}_{1}){\widehat{\psi}^{\vphantom{\dagger}}}_{\eta}({\bf r}_{2})\right],$ (30)

where the normal ordering is taken w.r.t. the effective vacuum state
$|0\rangle$. Due to the interaction, a SD wavefunction (Eqs. (19) and (II.2))
is no-longer an eigenstate of $\widehat{H}$ and representations of
$\widehat{H}$ are in general non-diagonal.

Let us start the search for the wavefunctions of $\widehat{H}$ by noticing
that the Fock space can be split into charge sectors, such that the ground
state energy of $\widehat{H}$ can be written as 444Assuming that the opposite
charge operator ($\widehat{Q}$) and the interacting Hamiltonian operator
($\widehat{H}$) commute, which is only valid for an appropriate definition of
the fermion-fermion interaction and with the adequate definition for the
normal ordering.

$E=\underset{|\Psi\rangle\in\mathcal{H}_{Q}}{\textrm{min}}\langle\Psi|\widehat{H}_{0}^{v}+\widehat{W}|\Psi\rangle.$
(31)

Thence, the wavefunctions must also be eigenstates of $\widehat{Q}$ with
eigenvalue $Q$ (i.e. $\widehat{Q}|\Psi\rangle=Q|\Psi\rangle$). Consequently,
$|\Psi\rangle$ can be constrained to have a given charge $Q$ (i.e. $\int d{\bf
r}\langle\Psi|\textrm{Tr}\left[\widehat{{\bf n}}_{1}({\bf r},{\bf
r})\right]|\Psi\rangle=Q$). Assuming that $\Psi$ is normalized to 1, a state
that belongs to $\mathcal{H}_{Q}$ for a particular charge sector $Q\geq 0$ can
be written (parameterized) as

$\displaystyle\lvert\Psi\rangle={}$
$\displaystyle\left(\sum_{I_{1},\dotsc,I_{Q}}c_{I_{1}\dotsc
I_{Q}}{\widehat{b}^{\;\dagger}}_{I_{1}}\dotsb{\widehat{b}^{\;\dagger}}_{I_{Q}}+\sum_{I_{1},\dotsc,I_{Q},I_{Q+1}}\;\sum_{R_{1}}c_{I_{1}\dotsc
I_{Q}I_{Q+1}R_{1}}{\widehat{b}^{\;\dagger}}_{I_{1}}\dotsb{\widehat{b}^{\;\dagger}}_{I_{Q}}{\widehat{b}^{\;\dagger}}_{I_{Q+1}}{\widehat{d}^{\;\dagger}}_{R_{1}}\right.$
$\displaystyle{}+\left.\sum_{I_{1},\dotsc,I_{Q},I_{Q+1},I_{Q+2}}\;\sum_{R_{1},R_{2}}c_{I_{1}\dotsc
I_{Q}I_{Q+1}I_{Q+2}R_{1}R_{2}}{\widehat{b}^{\;\dagger}}_{I_{1}}\dotsb{\widehat{b}^{\;\dagger}}_{I_{Q}}{\widehat{b}^{\;\dagger}}_{I_{Q+1}}{\widehat{b}^{\;\dagger}}_{I_{Q+2}}{\widehat{d}^{\;\dagger}}_{R_{1}}{\widehat{d}^{\;\dagger}}_{R_{2}}+\cdot\cdot\cdot\right)\lvert
0\rangle.$ (32)

Since 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\geq 0$ ones.
Thus, without any loss of generality and for practical reasons (i.e. as we are
usually interested in electronic wavefunctions more than in positronic ones),
from now 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\rangle$) 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 ($\lvert
0^{\textrm{CI}}\rangle$) as a linear combination of states with arbitrary
number of electron-positron pairs that belongs to the charge sector $Q=0$.
Thus, the CI vacuum wavefunction reads as

$\lvert
0^{\textrm{CI}}\rangle=\left(c_{0}+\sum_{I_{1}}\;\sum_{R_{1}}c_{I_{1}R_{1}}{\widehat{b}^{\;\dagger}}_{I_{1}}{\widehat{d}^{\;\dagger}}_{R_{1}}+\sum_{I_{1},I_{2}}\;\sum_{R_{1},R_{2}}c_{I_{1}I_{2}R_{1}R_{2}}{\widehat{b}^{\;\dagger}}_{I_{1}}{\widehat{b}^{\;\dagger}}_{I_{2}}{\widehat{d}^{\;\dagger}}_{R_{1}}{\widehat{d}^{\;\dagger}}_{R_{2}}+\cdot\cdot\cdot\right)\lvert
0\rangle,$ (33)

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
$\lvert 0^{\textrm{CI}}\rangle$ vacuum state. The vp effects (i.e. $\langle
0^{\textrm{CI}}\lvert\widehat{H}\lvert 0^{\textrm{CI}}\rangle=0$) can then be
omitted. In the usual situation, when the external potential is 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 latter do not contribute to the wavefunction. Moreover, in that limit the
external potential does not affect the definition of the vacuum state; thence,
the effective vacuum can be used independently of the interactions considered
for the electrons and for any different external potential employed.

###  II.5 The 1-RDM and its natural orbital representation

Let us comment on some properties of the 1-RDM obtained from a wavefunction
$|\Psi\rangle$ given by Eq. (II.4). Using the one-particle density matrix
operator Eq. (II.2) we define the matrix elements of the 1-RDM in the spinor
basis as

$\displaystyle n_{1,\nu,\mu}({\bf r},{\bf r}^{\prime})$
$\displaystyle={}\langle\Psi|\widehat{n}_{1,\nu,\mu}({\bf r},{\bf
r}^{\prime})|\Psi\rangle$
$\displaystyle={}\sum_{I,J}{}^{1}D_{I}^{J}\phi^{*}_{I,\mu}({\bf
r}^{\prime})\phi_{J,\nu}({\bf
r})+\sum_{I}\sum_{R}{}^{1}D_{I}^{R}\phi^{*}_{I,\mu}({\bf
r}^{\prime})\phi_{R,\nu}({\bf r})$
$\displaystyle{}+\sum_{R}\sum_{I}{}^{1}D_{R}^{I}\phi^{*}_{R,\mu}({\bf
r}^{\prime})\phi_{I,\nu}({\bf
r})+\sum_{R,S}{}^{1}D_{S}^{R}\phi^{*}_{R,\mu}({\bf
r}^{\prime})\phi_{S,\nu}({\bf r}),$ (34)

where we have introduced the 1-RDM coefficients
${}^{1}D_{I}^{J}=\langle\Psi|{\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{J}|\Psi\rangle$,
${}^{1}D_{I}^{R}=\langle\Psi|{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R}|\Psi\rangle$,
${}^{1}D_{R}^{I}=\langle\Psi|{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I}|\Psi\rangle$,
and
${}^{1}D_{S}^{R}=-\langle\Psi|{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{R}|\Psi\rangle$.
In a nonrelativistic context, the coefficients arising from the expectation
values
$\langle\Psi|{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R}|\Psi\rangle$
and
$\langle\Psi|{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I}|\Psi\rangle$
do not contribute since the wavefunction is formed by a fixed number of
particles and due to the absence of positrons in this limit, the last
contribution
$\langle\Psi|{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{R}|\Psi\rangle$
then trivially vanishes.

Collecting the matrix elements we form the corresponding 4$\times$4 matrix,
and after trivial algebra we obtain the full 1-RDM as

$\displaystyle{\bf n}_{1}({\bf r},{\bf r}^{\prime})$
$\displaystyle={}\sum_{I,J}{}^{1}D_{I}^{J}\begin{pmatrix}\phi^{*}_{I,1}({\bf
r}^{\prime})\phi_{J,1}({\bf r})&\phi^{*}_{I,2}({\bf
r}^{\prime})\phi_{J,1}({\bf r})&\phi^{*}_{I,3}({\bf
r}^{\prime})\phi_{J,1}({\bf r})&\phi^{*}_{I,4}({\bf
r}^{\prime})\phi_{J,1}({\bf r})\\\ \phi^{*}_{I,1}({\bf
r}^{\prime})\phi_{J,2}({\bf r})&\phi^{*}_{I,2}({\bf
r}^{\prime})\phi_{J,3}({\bf r})&\phi^{*}_{I,3}({\bf
r}^{\prime})\phi_{J,2}({\bf r})&\phi^{*}_{I,4}({\bf
r}^{\prime})\phi_{J,2}({\bf r})\\\ \phi^{*}_{I,1}({\bf
r}^{\prime})\phi_{J,3}({\bf r})&\phi^{*}_{I,2}({\bf
r}^{\prime})\phi_{J,3}({\bf r})&\phi^{*}_{I,3}({\bf
r}^{\prime})\phi_{J,3}({\bf r})&\phi^{*}_{I,4}({\bf
r}^{\prime})\phi_{J,3}({\bf r})\\\ \phi^{*}_{I,1}({\bf
r}^{\prime})\phi_{J,4}({\bf r})&\phi^{*}_{I,2}({\bf
r}^{\prime})\phi_{J,4}({\bf r})&\phi^{*}_{I,3}({\bf
r}^{\prime})\phi_{J,4}({\bf r})&\phi^{*}_{I,4}({\bf
r}^{\prime})\phi_{J,4}({\bf r})\end{pmatrix}$
$\displaystyle+\sum_{I}\sum_{R}{}^{1}D_{I}^{R}\left(\cdot\cdot\cdot\right)_{4\times
4}+\sum_{R}\sum_{I}{}^{1}D_{R}^{I}\left(\cdot\cdot\cdot\right)_{4\times
4}+\sum_{R,S}{}^{1}D_{S}^{R}\left(\cdot\cdot\cdot\right)_{4\times 4}$
$\displaystyle={}\sum_{I,J}{}^{1}D_{I}^{J}{\bm{\psi}}_{J}({\bf
r}){\bm{\psi}}^{\dagger}_{I}({\bf
r}^{\prime})+\sum_{I}\sum_{R}{}^{1}D_{I}^{R}{\bm{\psi}}_{R}({\bf
r}){\bm{\psi}}^{\dagger}_{I}({\bf
r}^{\prime})+\sum_{R}\sum_{I}{}^{1}D_{R}^{I}{\bm{\psi}}_{I}({\bf
r}){\bm{\psi}}^{\dagger}_{R}({\bf
r}^{\prime})+\sum_{R,S}{}^{1}D_{S}^{R}{\bm{\psi}}_{S}({\bf
r}){\bm{\psi}}^{\dagger}_{R}({\bf r}^{\prime}).$ (35)

The trace of the 1-RDM is equal to the charge $Q$, i.e.
$\sum_{A}{}^{1}D_{A}^{A}=Q$. Moreover, the coefficients ${}^{1}D_{A}^{B}$ form
an Hermitian matrix that can be diagonalized, which leads to the natural
orbital 555Could also be call as the natural spinor representation. (NO)
representation of the 1-RDM

${\bf n}_{1}({\bf r},{\bf r}^{\prime})=\sum_{A}n_{A}{\bm{\chi}}_{A}({\bf
r}){\bm{\chi}}^{\dagger}_{A}({\bf r}^{\prime}),$ (36)

where $\\{n_{A}\\}$ are known as the occupation numbers (ONs) that take values
in the interval between $\left[-1,1\right]$ and the NOs are written as
4-component spinors ${\bm{\chi}}_{A}({\bf r})=\sum_{B}{\bm{\psi}}_{B}({\bf
r})U_{BA}$ (the columns of the ${\bf U}$ matrix contain the eigenvectors
obtained from the diagonalization of the ${}^{1}{\bf D}$ matrix).

In the particular case of a SD wavefunction (19) the (canonical) ${\bm{\psi}}$
spinors coincide with the NOs with ONs that are then either 0 or $\pm$1\.
Consequently, the 1-RDM for $|\Phi\rangle$ reads as

${\bf n}^{\textrm{SD}}_{1}({\bf r},{\bf
r}^{\prime})=\sum_{I}^{N_{e}}{\bm{\psi}}_{I}({\bf
r}){\bm{\psi}}^{\dagger}_{I}({\bf
r}^{\prime})-\sum_{R}^{N_{p}}{\bm{\psi}}_{R}({\bf
r}){\bm{\psi}}^{\dagger}_{R}({\bf r}^{\prime}).$ (37)

###  II.6 Introducing ReRDMFT through the constrained-search formalism

Considering the constrained-search formalism [66, 67, 68] we define the
universal functional of the 1-RDM $W\bigl{[}{\bf n}_{1}\bigr{]}$ for
$N$-representable 1-RDMs (${\bf n}_{1}\in\mathcal{D}_{Q}$) that come from a
state $|\Psi\rangle\in\mathcal{H}_{Q}$ 666We currently have no easy
characterization of the set of $N$-representable 1-RDMs $\mathcal{D}_{Q}$. We
at least know that the trace of these 1-RDMs should be equal to the charge
$Q$. Other necessary conditions will be investigated in future work.,

$\displaystyle W\left[{\bf n}_{1}\right]$
$\displaystyle=\underset{|\Psi\rangle\in\mathcal{H}_{Q}({\bf
n}_{1})}{\textrm{min}}\langle\Psi|\widehat{W}|\Psi\rangle$
$\displaystyle=\langle\Psi\left[{\bf n}_{1}\right]|\widehat{W}|\Psi\left[{\bf
n}_{1}\right]\rangle,$ (38)

where $\mathcal{H}_{Q}({\bf n}_{1})$ is the set of states
$|\Psi\rangle\in\mathcal{H}_{Q}$ that yield a constrained 1-RDM (${\bf
n}_{1}$), and $|\Psi\left[{\bf n}_{1}\right]\rangle$ designates the state that
minimizes this energy contribution 777Whether we have a true minimum, is an
open question, but we suspect it to be the case like in the nonrelativistic
limit [67].. Thence, we may write the ground state energy functional of the
1-RDM as

$E_{Q}=\underset{{\bf n}_{1}\in\mathcal{D}_{Q}}{\textrm{min}}\left[W\left[{\bf
n}_{1}\right]+\int d{\bf r}d{\bf r}^{\prime}\textrm{Tr}\left[\left(\delta({\bf
r}-{\bf r}^{\prime})\widehat{\bf T}_{D}({\bf r})+{\bf
v}_{\textrm{ext}}^{\textrm{nl}}({\bf r}^{\prime},{\bf r})\right){\bf
n}_{1}({\bf r},{\bf r}^{\prime})\right]\right].$ (39)

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 density

$n({\bf r})=\int d{\bf r}^{\prime}\delta({\bf r}-{\bf
r}^{\prime})\textrm{Tr}\left[{\bf n}_{1}({\bf r},{\bf r}^{\prime})\right].$
(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.,

$\displaystyle F\left[n\right]$
$\displaystyle=\underset{|\Psi\rangle\in\mathcal{H}_{Q}(n)}{\textrm{min}}\langle\Psi|\widehat{T}_{D}+\widehat{W}|\Psi\rangle$
$\displaystyle=\langle\Psi\left[n\right]|\widehat{T}_{D}+\widehat{W}|\Psi\left[n\right]\rangle,$
(41)

where $\mathcal{H}_{Q}(n)$ is the set of states
$|\Psi\rangle\in\mathcal{H}_{Q}$ that yield a constrained opposite charge
density ($n$), and $|\Psi\left[n\right]\rangle$ designates the state that
minimizes this energy contribution. Consequently, the energy functional within
the constrained search formalism for $N$-representable opposite charge
densities ($n\in\mathcal{D}^{\textrm{dens}}_{Q}$) reads as

$E_{Q}=\underset{n\in\mathcal{D}^{\textrm{dens}}_{Q}}{\textrm{min}}\left[F\left[n\right]+\int
d{\bf r}v_{\textrm{ext}}({\bf r})n({\bf r})\right].$ (42)

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 ${\bf n}_{1}$. For instance, adding
an energy contribution arising from an external vector field (${\bf
A}^{\textrm{ext}}$) to the above expression is straight forward: only the term
$\delta({\bf r}-{\bf r}^{\prime})c{\bm{\alpha}}\cdot{\bf
A}^{\textrm{ext}}({\bf r})$ needs to be added and fits the general form of
${\bf v}^{\textrm{nl}}_{\textrm{ext}}({\bf r}^{\prime},{\bf r})$ in (9). Thus,
the functional of the 1-RDM contains more information than the functional of
the current (i.e. ${\bf j}({\bf r})=c{\bm{\psi}}^{\dagger}({\bf
r}){\bm{\alpha}}{\bm{\psi}}({\bf r})=\textrm{Tr}\left[c{\bm{\alpha}{\bf
n}_{1}({\bf r},{\bf r})}\right]$), which clearly shows the generality of the
present theory and makes it interesting also for the description of magnetism
in molecular 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 the constrained-search formalism is presented in the Appendix B,
where we discuss the extension of Gilbert’s theorem [44] to the relativistic
domain.

##  III The no-pair approximation and ReRDMFT

###  III.1 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\rangle$888Recall that this vacuum state is
given using some orthonormal basis ${\bm{\psi}}_{A}({\bf r})$ that is the
initial guess. Also notice that the operators ${\widehat{b}^{\;\dagger}}_{I}$,
${\widehat{d}^{\;\dagger}}_{R}$, ${\widehat{b}^{\vphantom{\dagger}}}_{I}$, and
${\widehat{d}^{\vphantom{\dagger}}}_{R}$ are the initial ones., which leads to
$\langle 0|\widehat{H}|0\rangle=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

$\widehat{{H}}=\widehat{\widetilde{H}}+\widehat{\widetilde{V}}^{\textrm{vp}}+\widetilde{E}_{0},$
(43)

where

$\displaystyle\widehat{\widetilde{V}}^{\textrm{vp}}={}$
$\displaystyle\widehat{\widetilde{V}}_{H}^{\textrm{vp}}+\widehat{\widetilde{V}}_{\textrm{x}}^{\textrm{vp}}$
$\displaystyle={}$ $\displaystyle\sum_{\mu,\tau}\int d{\bf
r}_{1}\left[\sum_{\nu,\eta}\int d{\bf r}_{2}W_{\mu,\nu,\tau,\eta}({\bf
r}_{1},{\bf r}_{2})\widetilde{n}^{\textrm{vp}}_{1,\eta,\nu}({\bf r}_{2},{\bf
r}_{2})\right]\widehat{\widetilde{n}}_{1,\tau,\mu}({\bf r}_{1},{\bf r}_{1})$
(44) $\displaystyle-$ $\displaystyle\sum_{\mu,\nu,\tau,\eta}\int\int d{\bf
r}_{1}d{\bf r}_{2}W_{\mu,\nu,\tau,\eta}({\bf r}_{1},{\bf
r}_{2})\widetilde{n}^{\textrm{vp}}_{1,\eta,\mu}({\bf r}_{2},{\bf
r}_{1})\widehat{\widetilde{n}}_{1,\tau,\nu}({\bf r}_{1},{\bf r}_{2}),$
$\displaystyle\widetilde{E}_{0}={}$
$\displaystyle\langle\widetilde{0}_{v}|\widehat{{H}}|\widetilde{0}_{v}\rangle=\widetilde{E}_{0}^{0}+\frac{1}{2}\int
d{\bf r}_{1}d{\bf r}_{2}\textrm{Tr}\left[{\bf W}({\bf r}_{1},{\bf
r}_{2})\widetilde{{\bf n}}^{\textrm{vp}}_{2}({\bf r}_{1},{\bf r}_{2})\right],$
(45)

where we have introduced the vp pair-density matrix (written in terms of its
components) as

$\displaystyle\widetilde{{n}}^{\textrm{vp}}_{2,\tau,\eta,\mu,\nu}({\bf
r}_{1},{\bf r}_{2})=\widetilde{n}^{\textrm{vp}}_{1,\eta,\nu}({\bf r}_{2},{\bf
r}_{2})\widetilde{n}^{\textrm{vp}}_{1,\tau,\mu}({\bf r}_{1},{\bf
r}_{1})-\widetilde{n}^{\textrm{vp}}_{1,\tau,\nu}({\bf r}_{1},{\bf
r}_{2})\widetilde{n}^{\textrm{vp}}_{1,\eta,\mu}({\bf r}_{2},{\bf r}_{1}).$
(46)

Let us highlight that the definition of the vacuum as an effective vacuum
($|{0}\rangle$ or $|\widetilde{0}_{v}\rangle$, i.e. at the mean-field level),
leads to only Hartree ($\widehat{\widetilde{V}}^{\textrm{vp}}_{H}$) and
exchange ($\widehat{\widetilde{V}}^{\textrm{vp}}_{x}$) like terms for the
interaction between the fermions (electrons or positrons) present in the
$\lvert\Psi\rangle$ and the polarization of the vacuum. Indeed, the
$\widetilde{\bf n}_{1}^{\textrm{vp}}$ can be regarded as the difference
between the 1-RDM obtained from a SD wavefunction built from all positronic
$\widetilde{\bm{\psi}}$ states, i.e.,

$\lvert\widetilde{\Phi}\rangle={\widehat{\widetilde{d}}^{\;\dagger}}_{R_{1}}\dotsb{\widehat{\widetilde{d}}^{\;\dagger}}_{R_{M}}\lvert\widetilde{0}_{v}\rangle,$
(47)

and the same kind of function build with all positronic ${\bm{\psi}}$ states,
i.e.,

$\lvert\Phi\rangle={\widehat{d}^{\;\dagger}}_{R_{1}}\dotsb{\widehat{d}^{\;\dagger}}_{R_{M}}\lvert
0\rangle,$ (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 are present.

The effective vacuum state $|\widetilde{0}_{v}\rangle$ that minimizes Eq. 45
can actually produce an energy that diverges to $-\infty$ 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, 73, 74, 75,
76, 77, 78, 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 other
words, 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, which
makes 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, the use 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.

###  III.2 The no-pair approximation including vacuum polarization ReRDFMT
approach

The no-pair (np) approximation [80, 48] 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 $\sim
150$keV) than required to create electron-positron pairs ($>1$MeV) [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. $N_{e}=N$ and $N_{p}=0$). Within the np approximation only the PS are
populated while the NS remain empty. The np vp ground state energy reads as

$E_{N}^{\textrm{npvp}}=\underset{|\Psi_{+}\rangle\in\widetilde{\mathcal{H}}^{(N,0)}}{\textrm{min}}\langle\Psi_{+}|\widehat{H}|\Psi_{+}\rangle,$
(49)

where the minimization is over normalized states in the set
$\widetilde{\mathcal{H}}^{(N,0)}=\textrm{e}^{\widehat{\kappa}}{\mathcal{H}}^{(N,0)}$
that is the set of states generated by all spinor rotations of $N$-electron
states. A state $|\Psi_{+}\rangle\in\widetilde{\mathcal{H}}^{(N,0)}$ can be
written as

$\displaystyle|\Psi_{+}\rangle$
$\displaystyle=\textrm{e}^{\widehat{\kappa}}\sum_{\mathclap{I_{1},\dotsc,I_{N}}}c_{I_{1}\dotsc
I_{N}}{\widehat{b}^{\;\dagger}}_{I_{1}}\cdot\cdot\cdot{\widehat{b}^{\;\dagger}}_{I_{N}}|0\rangle$
$\displaystyle=\sum_{\mathclap{I_{1},\dotsc,I_{N}}}c_{I_{1}\dotsc
I_{N}}{\widehat{\widetilde{b}}^{\;\dagger}}_{I_{1}}\cdot\cdot\cdot{\widehat{\widetilde{b}}^{\;\dagger}}_{I_{N}}|\widetilde{0}_{v}\rangle.$
(50)

Since the state $|\Psi_{+}\rangle$ can also be written as
$|\Psi_{+}\rangle=\widehat{\widetilde{P}}_{+}|\Psi\rangle$, where
$|\Psi\rangle\in\mathcal{H}_{N}$ is an arbitrary state constrained to have $N$
negative charges, and $\widehat{\widetilde{P}}_{+}$ projects onto the
$N$-electron Hilbert space constructed by the
$\\{{\widehat{\widetilde{b}}^{\;\dagger}}_{I}\\}$ operators. Thence, 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 $N_{e}$ through
Eqs. (III.1) and (45). Lastly, using the constrained search formalism [66, 67,
68] and Eq. (49) we can readily introduce npvp-ReRDMFT, where vp effects are
accounted at the effective QED level 999Let 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)..

For making npvp-ReRDMFT more explicit, let us denote the np 1-RDM as

${\bf n}_{1}^{+}({\bf r},{\bf
r}^{\prime})=\sum_{I,J}{}^{1}\widetilde{D}_{I}^{J}\widetilde{\bm{\psi}}_{J}({\bf
r})\widetilde{\bm{\psi}}^{\dagger}_{I}({\bf r}^{\prime}),$ (51)

where
${}^{1}\widetilde{D}_{I}^{J}=\langle\Psi_{+}|{\widehat{\widetilde{b}}^{\;\dagger}}_{I}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{J}|\Psi_{+}\rangle$.
Making use of ${\bf n}_{1}^{+}$ and $\widetilde{\bf n}_{1}^{\textrm{vp}}$ all
one-body interactions of Eq. (49) (i.e. the electronic and the vp effects) can
be readily evaluated as

$\langle\Psi_{+}|\widehat{T}_{D}+\widehat{V}^{\textrm{nl}}_{\textrm{ext}}|\Psi_{+}\rangle=\int
d{\bf r}d{\bf r}^{\prime}\textrm{Tr}\left[\left(\delta({\bf r}-{\bf
r}^{\prime})\widehat{\bf T}_{D}({\bf r})+{\bf
v}_{\textrm{ext}}^{\textrm{nl}}({\bf r}^{\prime},{\bf r})\right)\left({\bf
n}_{1}^{\textrm{+}}({\bf r},{\bf r}^{\prime})+\widetilde{\bf
n}_{1}^{\textrm{vp}}({\bf r},{\bf r}^{\prime})\right)\right].$ (52)

Also, the explicit vp contributions in terms of ${\bf n}_{1}^{\textrm{+}}$ and
$\widetilde{\bf n}_{1}^{\textrm{vp}}$ to
$\langle\Psi_{+}|\widehat{W}|\Psi_{+}\rangle$ is known through Eqs. (III.1)
and (45). Indeed, only the pure electronic contribution to
$\langle\Psi_{+}|\widehat{\widetilde{W}}|\Psi_{+}\rangle$ in terms of ${\bf
n}_{1}^{+}$ is unknown and it needs to be approximated (see Appendix C for
more details).

###  III.3 The no-pair approximation ReRDFMT approach

The vp contribution to the total energy (i.e. $\widetilde{E}_{0}$ and
$\langle\Psi_{+}|\widehat{\widetilde{V}}^{\textrm{vp}}|\Psi_{+}\rangle$) is
usually neglected in practical applications of 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

$E_{N}^{\textrm{np}}=\underset{|\Psi_{+}\rangle\in{\mathcal{\widetilde{H}}}^{(N,0)}}{\textrm{min}}\langle\Psi_{+}|\widehat{\widetilde{H}}|\Psi_{+}\rangle,$
(53)

where the Hamiltonian is written in normal ordering w.r.t. the floating vacuum
state $|\widetilde{0}_{v}\rangle$. From Eq. (53) we may write, within the
constrained-search formalism, the functional expression for $N$-representable
np 1-RDMs (${\bf n}_{1}^{+}\in\mathcal{D}_{N}^{+}$) that come from a state
$|\Psi_{+}\rangle\in{\mathcal{\widetilde{H}}}^{(N,0)}$

$\displaystyle W^{\textrm{np}}\left[{\bf n}_{1}^{+}\right]$
$\displaystyle=\underset{|\Psi_{+}\rangle\in{\mathcal{\widetilde{H}}}^{(N,0)}({\bf
n}_{1}^{+})}{\textrm{min}}\langle\Psi_{+}|\widehat{\widetilde{W}}|\Psi_{+}\rangle$
$\displaystyle=\langle\Psi_{+}\left[{\bf
n}_{1}^{+}\right]|\widehat{\widetilde{W}}|\Psi_{+}\left[{\bf
n}_{1}^{+}\right]\rangle,$ (54)

where ${\mathcal{\widetilde{H}}}^{(N,0)}({\bf n}_{1}^{+})$ is the set of
states $|\Psi_{+}\rangle\in\widetilde{H}^{(N,0)}$ that yield a constrained
1-RDM (${\bf n}_{1}^{+}$), and $|\Psi_{+}\left[{\bf n}_{1}^{+}\right]\rangle$
designates 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

$E_{N}^{\textrm{np}}=\underset{{\bf
n}_{1}^{+}\in\mathcal{D}_{N}^{+}}{\textrm{min}}\left[{W}^{\textrm{np}}\left[{\bf
n}_{1}^{+}\right]+\int d{\bf r}d{\bf
r}^{\prime}\textrm{Tr}\left[\left(\delta({\bf r}-{\bf r}^{\prime})\widehat{\bf
T}_{D}({\bf r})+{\bf v}_{\textrm{ext}}^{\textrm{nl}}({\bf r}^{\prime},{\bf
r})\right){\bf n}_{1}^{+}({\bf r},{\bf r}^{\prime})\right]\right].$ (55)

The np approximation resembles the nonrelativistic result because the trace of
the ${\bf n}_{1}^{+}$ matrix is equal to the number of electrons, i.e.
$\textrm{Tr}\left[{}^{1}\widetilde{\bf D}\right]=N_{e}$. The ${\bf n}_{1}^{+}$
matrix can also be expressed in the NO representation,

${\bf n}_{1}^{+}({\bf r},{\bf
r}^{\prime})=\sum_{I}n_{I}\widetilde{\bm{\chi}}_{I}({\bf
r})\widetilde{\bm{\chi}}_{I}^{\dagger}({\bf r}^{\prime}),$ (56)

with $\sum_{I}n_{I}=N_{e}$ which allows us to introduce np-ReNOFT.

To conclude this section, as it is mentioned in Ref. 45, the np approximation
regains the concept of an $N$-electron wavefunction,

$\Psi_{+}({\bf r}_{1},{\bf r}_{2},\dotsc,{\bf
r}_{N})=\sum_{\mathclap{I_{1}<\dotsc<I_{N}}}c_{I_{1}\dotsc
I_{N}}\widetilde{\bm{\psi}}_{I_{1}}({\bf
r}_{1})\wedge\dotsb\wedge\widetilde{\bm{\psi}}_{I_{N}}({\bf r}_{N}),$ (57)

where $\wedge$ denotes the normalized antisymmetrized tensor product.

###  III.4 The $\langle\Psi_{+}|\widehat{\widetilde{W}}|\Psi_{+}\rangle\;$
term as an explicit functional of the second-order reduced density matrix.

To conclude this section, let us comment on some properties of the expectation
value of the $\widehat{\widetilde{W}}$ operator. In the np approximation this
term reads as

$\langle\Psi_{+}|\widehat{\widetilde{W}}|\Psi_{+}\rangle=\sum_{I,J,K,L}{}^{2}D_{IJ}^{KL}\int
d{\bf r}_{1}d{\bf r}_{2}\textrm{Tr}\left[{\bf W}({\bf r_{1}},{\bf
r}_{2})(\widetilde{{\bm{\psi}}}_{L}({\bf
r}_{1})\otimes\widetilde{{\bm{\psi}}}_{K}({\bf
r}_{2}))(\widetilde{{\bm{\psi}}}^{\dagger}_{I}({\bf
r}_{2})\otimes\widetilde{{\bm{\psi}}}^{\dagger}_{J}({\bf r}_{1}))\right]$ (58)

where $\otimes$ denotes the tensor product, and
${}^{2}D_{IJ}^{KL}=\frac{1}{2}\langle\Psi_{+}|{\widehat{\widetilde{b}}^{\;\dagger}}_{I}{\widehat{\widetilde{b}}^{\;\dagger}}_{J}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{L}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{K}|\Psi_{+}\rangle$
is the np second-order reduced density matrix (2-RDM) element. Clearly, the
$\langle\Psi_{+}|\widehat{\widetilde{W}}|\Psi_{+}\rangle$ contribution is an
explicit functional of the 2-RDM. In the np approximation, the properties of
the np 2-RDM are the same as in the nonrelativistic case. In particular, the
trace of this matrix reads

$\textrm{Tr}\left[{}^{2}{\bf
D}\right]=\sum_{I,J}{}^{2}D_{IJ}^{IJ}=\frac{N_{e}(N_{e}-1)}{2}.$ (59)

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, practical applications of relativistic quantum
chemistry/physics normally retain only part of the contributions to
$\widehat{\widetilde{W}}$. Indeed, the Coulomb–Breit interaction (29) consists
of two physically distinct contributions: one representing the magnetic
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 Coulomb
interaction to yield

$W_{\mu,\nu,\tau,\eta}({\bf r}_{1},{\bf
r}_{2})=\frac{1}{r_{12}}\left[\delta_{\mu,\tau}\delta_{\nu,\eta}-({\bm{\alpha}}_{\mu,\tau}\cdot{\bm{\alpha}}_{\nu,\eta})\right],$
(60)

which will be used in the rest of this work (i.e. inserting it in Eq. (58)).

Before proceeding, let us stress that both energy expressions
$E_{N}^{\textrm{npvp}}$ and $E_{N}^{\textrm{np}}$ can be minimized since the
Hamiltonian is written using normal ordering. Hence, they are bounded from
below and the state $|\Psi_{+}\rangle$ can actually be considered as a
projected state from the $|\Psi\rangle\in\mathcal{H}_{N}$.

###  III.5 np-ReNOFT functional approximations by imposing Kramers’ symmetry

Since the np approximation is the most common initial starting point for
relativistic calculations, we propose the construction of functional
approximations in this framework.

The explicit dependence of the $W\bigl{[}{\bf n}^{+}_{1}\bigr{]}$ functional
in terms of the ${\bf n}_{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, 88, 89, 90] are built by imposing the so-called
$N$-representability conditions [83, 84, 91] and using only up to two
different indices to define 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 not included [85],
which leads to a restricted formulation of spin magnetic moments. Then, a pair
of NOs forms a Kramers’ pair ($i,\bar{i}$) if they transform as
$\widehat{\mathcal{K}}\widetilde{{\bm{\chi}}}_{i}=\widetilde{{\bm{\chi}}}_{\bar{i}}$
and
$\widehat{\mathcal{K}}\widetilde{{\bm{\chi}}}_{\bar{i}}=-\widetilde{{\bm{\chi}}}_{i}$,
where we have defined the time-reversal operator as

$\widehat{\mathcal{K}}=-\textrm{i}{\begin{pmatrix}{\bm{\sigma}}_{y}&{\bm{0}}_{2}\\\
{\bm{0}}_{2}&{\bm{\sigma}}_{y}\end{pmatrix}}\widehat{\mathcal{K}}_{0},$ (61)

and where $\widehat{\mathcal{K}}_{0}$ is the complex conjugation operator.
From now on we form two subsets of PS by splitting all spinors into Kramers’
pairs 101010Let 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.. Actually, the proper definition of Kramers’ pairs is not unique because
the pair ($i,\bar{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 spin quantization 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 111111The same kind of
partitioning of the spinor states can be applied to the NS..

Imposing the Kramers’ symmetry we may rewrite the $\widehat{\widetilde{H}}$
operator (including only the Coulomb and the Gaunt interactions) in the NOs
representation as [55]

$\displaystyle\widehat{\widetilde{H}}^{\textrm{KR}}={}$
$\displaystyle\sum_{i}h_{ii}\widehat{X}^{+}_{ii}+\frac{1}{2}\sum_{i,j,k,l}\left[\langle
ij|kl\rangle\widehat{x}^{++}_{ik,jl}-\langle
ij|\bm{\alpha}\odot\bm{\alpha}|kl\rangle\widehat{x}^{--}_{ik,jl}+\langle\bar{i}j|kl\rangle\widehat{x}^{++}_{\bar{i}k,jl}-\langle\bar{i}j|\bm{\alpha}\odot\bm{\alpha}|kl\rangle\widehat{x}^{--}_{\bar{i}k,jl}\right.$
$\displaystyle{}+\left.\langle
ij|\bar{k}l\rangle\widehat{x}^{++}_{i\bar{k},jl}-\langle
ij|\bm{\alpha}\odot\bm{\alpha}|\bar{k}l\rangle\widehat{x}^{--}_{i\bar{k},jl}\right]+\frac{1}{4}\sum_{i,j,k,l}\left[\langle\bar{i}j|k\bar{l}\rangle\widehat{x}^{++}_{\bar{i}k,j\bar{l}}-\langle\bar{i}j|\bm{\alpha}\odot\bm{\alpha}|k\bar{l}\rangle\widehat{x}^{--}_{\bar{i}k,j\bar{l}}\right]$
(62)
$\displaystyle{}+\frac{1}{8}\sum_{i,j,k,l}\left[\langle\bar{i}\bar{j}|kl\rangle\widehat{x}^{++}_{\bar{i}k,\bar{j}l}-\langle\bar{i}\bar{j}|\bm{\alpha}\odot\bm{\alpha}|kl\rangle\widehat{x}^{--}_{\bar{i}k,\bar{j}l}+\langle
ij|\bar{k}\bar{l}\rangle\widehat{x}^{++}_{i\bar{k},j\bar{l}}-\langle
ij|\bm{\alpha}\odot\bm{\alpha}|\bar{k}\bar{l}\rangle\widehat{x}^{--}_{i\bar{k},j\bar{l}}\right],$

where
$\bm{\alpha}\odot\bm{\alpha}=\bm{\alpha}_{x}\otimes\bm{\alpha}_{x}+\bm{\alpha}_{y}\otimes\bm{\alpha}_{y}+\bm{\alpha}_{z}\otimes\bm{\alpha}_{z}$
and the integrals are defined as

$\displaystyle h_{II}$ $\displaystyle=\int d{{\bf r}}d{{\bf
r}^{\prime}}\widetilde{{\bm{\chi}}}^{\dagger}_{I}({\bf
r}^{\prime})\left(\delta({\bf r}-{\bf r}^{\prime})\widehat{\bf T}_{D}({\bf
r})+{\bf v}^{\textrm{nl}}_{\textrm{ext}}({\bf r},{\bf
r}^{\prime})\right)\widetilde{{\bm{\chi}}}_{I}({\bf r}),$ (63a)
$\displaystyle\langle IJ|KL\rangle$ $\displaystyle=\int d{{\bf r}}d{{\bf
r}^{\prime}}\frac{\bigl{(}\widetilde{{\bm{\chi}}}^{\dagger}_{I}({\bf
r})\otimes\widetilde{{\bm{\chi}}}_{J}^{\dagger}({\bf
r}^{\prime})\bigr{)}\,\bigl{(}\widetilde{{\bm{\chi}}}_{K}({\bf
r})\otimes\widetilde{{\bm{\chi}}}_{L}({\bf r}^{\prime})\bigr{)}}{\lvert{\bf
r}^{\prime}-{\bf r}\rvert},$ (63b) $\displaystyle\langle
IJ|\bm{\alpha}\odot\bm{\alpha}|KL\rangle$ $\displaystyle=\int d{{\bf r}}d{{\bf
r}^{\prime}}\frac{\bigl{(}\widetilde{{\bm{\chi}}}^{\dagger}_{I}({\bf
r})\otimes\widetilde{{\bm{\chi}}}_{J}^{\dagger}({\bf
r}^{\prime})\bigr{)}\,\bigl{[}{\bm{\alpha}}\odot{\bm{\alpha}}\bigr{]}\,\bigl{(}\widetilde{{\bm{\chi}}}_{K}({\bf
r})\otimes\widetilde{{\bm{\chi}}}_{L}({\bf r}^{\prime})\bigr{)}}{\lvert{\bf
r}^{\prime}-{\bf r}\rvert},$ (63c)

$\widehat{X}^{s}_{ij}=(1+s\widehat{T}_{ij}){\widehat{\widetilde{b}}^{\;\dagger}}_{i}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{j}$,
$\widehat{x}^{s_{1}s_{2}}_{IK,JL}=(1+s_{1}\widehat{T}_{ik})(1+s_{2}\widehat{T}_{jl}){\widehat{\widetilde{b}}^{\;\dagger}}_{I}{\widehat{\widetilde{b}}^{\;\dagger}}_{J}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{L}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{K}=\widehat{x}^{s_{2}s_{1}}_{JL,IK}$
with
$\widehat{T}_{ij}{\widehat{\widetilde{b}}^{\;\dagger}}_{i}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{j}={\widehat{\widetilde{b}}^{\;\dagger}}_{\bar{j}}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{\bar{i}}$
and
$\widehat{T}_{ij}{\widehat{\widetilde{b}}^{\;\dagger}}_{\bar{i}}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{j}=-{\widehat{\widetilde{b}}^{\;\dagger}}_{\bar{j}}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{i}$
121212The uppercase indices that enter in the definition of
${\widehat{x}^{s_{1}s_{2}}_{IK,JL}}$ are readily used in the creation and
annihilation operators, but they are used as (lowercase) unbarred indices in
$T_{ij}$ operator. Actually, the $T_{ij}$ is only defined with unbarred
indices.. Let us remark that the $+$ and $-$ signs placed in Eq. (III.5) do
not refer to positive or negative energy spinors; they only enter in the
definition of $\widehat{X}^{s}_{IJ}$ and $\widehat{x}^{s_{1}s_{2}}_{IJ,KL}$.
Using the above Hamiltonian, the np energy reads as
$E^{\textrm{np}}_{N}=\langle\Psi_{+}|\widehat{\widetilde{H}}{}^{\textrm{KR}}|\Psi_{+}\rangle$.
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 in the 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 [95, 90, 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

$\displaystyle E_{N}^{\textrm{np}}$
$\displaystyle\approx{}\sum_{i}h_{ii}\left(n_{i}+n_{\bar{i}}\right)+\sum_{i,j}\left({}^{2}D_{ij}^{ij}+{}^{2}D_{\bar{i}j}^{\bar{i}j}+{}^{2}D_{i\bar{j}}^{i\bar{j}}+{}^{2}D_{\bar{i}\bar{j}}^{\bar{i}\bar{j}}\right)J_{ij}$
$\displaystyle{}-\sum_{i,j}\left({}^{2}D_{ij}^{ij}-{}^{2}D_{\bar{i}j}^{\bar{i}j}-{}^{2}D_{i\bar{j}}^{i\bar{j}}+{}^{2}D_{\bar{i}\bar{j}}^{\bar{i}\bar{j}}\right)J^{G}_{ij}+\sum_{i,j}\left[\left({}^{2}D_{ij}^{ji}+{}^{2}D_{\bar{j}\bar{i}}^{\bar{i}\bar{j}}\right)\left(K_{ij}-K^{G}_{ij}\right)\right]$
$\displaystyle{}+\frac{1}{2}\sum_{i,j}\left[\left({}^{2}D_{\bar{i}j}^{j\bar{i}}+{}^{2}D_{\bar{j}i}^{i\bar{j}}+{}^{2}D_{j\bar{i}}^{\bar{i}j}+{}^{2}D_{i\bar{j}}^{\bar{j}i}\right)\left(L_{ij}-L^{G}_{ij}\right)\right]+\sum_{i\neq
j}\left[\left({}^{2}D_{i\bar{i}}^{j\bar{j}}+{}^{2}D_{\bar{j}j}^{\bar{i}i}\right)\left(K_{ij}+K^{G}_{ij}\right)\right]$
$\displaystyle{}-\frac{1}{2}\sum_{i\neq
j}\left[\left({}^{2}D_{\bar{j}j}^{i\bar{i}}+{}^{2}D_{\bar{i}i}^{j\bar{j}}+{}^{2}D_{j\bar{j}}^{\bar{i}i}+{}^{2}D_{i\bar{i}}^{\bar{j}j}\right)\left(L_{ij}+L^{G}_{ij}\right)\right],$
(64)

where $J_{ij}=\langle ij|ij\rangle$, $J^{G}_{ij}=\langle
ij|{\bm{\alpha}}\odot{\bm{\alpha}}|ij\rangle$, $K_{ij}=\langle ij|ji\rangle$,
$K^{G}_{ij}=\langle ij|{\bm{\alpha}}\odot{\bm{\alpha}}|ji\rangle$,
$L_{ij}=\langle\bar{i}j|j\bar{i}\rangle$ (notice that $L_{ii}=0$ [85]), and
$L^{G}_{ij}=\langle\bar{i}j|{\bm{\alpha}}\odot{\bm{\alpha}}|j\bar{i}\rangle$.
Interestingly, some NOFT approximations (i.e. GNOF/PNOF$x$ approximations)
employed 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 similarly in
the evaluation of the $N$-representability conditions, etc.); hence, we have
adopted the same strategy that is 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$\left[\textrm{np(vp)-ReRDMFT}\right]$). But, it also facilitates
the search for analytic constraints on the approximated 2-RDM matrix elements
upon evaluation of 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.

####  III.5.1 The Dirac–Hartree–Fock approximation

The np Dirac–Hartree–Fock (DHF) approximations uses a SD wavefunction to
approximate $|\Psi_{+}\rangle$ (occupying only PS). Actually, the DHF energy
can be readily written from Eq. (III.5) for ‘spin-compensated’ systems where
$n_{i}=n_{\bar{i}}$ and defining the matrix elements of the 2-RDM with the SD
approximation (${}^{2}{\bf D}^{\textrm{SD}}$)

$({}^{2}D_{IJ}^{KL})^{\textrm{SD}}=\frac{n_{I}n_{J}}{2}(\delta_{IK}\delta_{JL}-\delta_{IL}\delta_{JK}).$
(65)

This approximation produces the correct symmetry and antisymmetry properties
of the 2-RDM [98]. Thence, the np DHF energy reads as

$\displaystyle E^{\textrm{np,DHF}}$
$\displaystyle=2\sum_{i}h_{ii}n_{i}+\sum_{i,j}{n_{i}n_{j}}\left[2J_{ij}-\left(K_{ij}-K^{G}_{ij}\right)-\left(L_{ij}-L^{G}_{ij}\right)\right]$
$\displaystyle=2\sum^{N_{e}/2}_{i}h_{ii}+\sum^{N_{e}/2}_{i,j}\left[2J_{ij}-\left(K_{ij}-K^{G}_{ij}\right)-\left(L_{ij}-L^{G}_{ij}\right)\right].$
(66)

In the last expression we have use $n_{i}=n_{\bar{i}}=1$ for the PS that form
the SD wavefunction, and 0 otherwise. An efficient algorithm 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_{\mu,\nu,\tau,\eta}({\bf r}_{1},{\bf
r}_{2})=\frac{\delta_{\mu,\tau}\delta_{\nu,\eta}}{r_{12}}$) 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 ($n_{i}=n_{\bar{i}}$). The first
family of approximations is based on mimicking the exchange-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}D_{IJ}^{KL}={}^{2}D_{IJ}^{KL}(n_{I},n_{J},n_{K},n_{L})$ (67)

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

####  III.5.2 $f(n_{I},n_{J})$-functional approximations

The family of nonrelativistic $f(n_{I},n_{J})$-functional approximations is
introduced as an attempt to separate the electron-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 ${n_{I}n_{J}}$ in the second term on the r.h.s. of the DHF energy
by a $f(n_{I},n_{J})$ 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(n_{I},n_{J})$-functional
approximations (i.e. the 2RDM elements) reads as

$\left({}^{2}D_{IJ}^{KL}\right)^{\textrm{approx.}}=\frac{n_{I}n_{J}}{2}\delta_{IK}\delta_{JL}-\frac{f(n_{I},n_{J})}{2}\delta_{IL}\delta_{JK},$
(68)

In Table 2 we have collected some $f(n_{I},n_{J})$ functions that lead to some
representative functional approximations.

Table 2: Representative $f(n_{I},n_{J})$-functional approximations. Functional | $f(n_{I},n_{J})$ | case | Ref.
---|---|---|---
MBB | $\sqrt{n_{I}n_{J}}$ | all | 101, 102, 103
Power | $\left({n_{I}n_{J}}\right)^{\alpha}$ | all | 104, 105, 106
GU | $\sqrt{n_{I}n_{J}}$ | $I\neq J$ | 107
$n_{I}n_{J}$ | otherwise
MLSIC | $n_{I}n_{J}\frac{a_{0}+a_{1}n_{I}n_{J}}{1+b_{1}n_{I}n_{J}}$ | $I\neq J$ | 108
$n_{I}n_{J}$ | otherwise
$a_{0}=1298.78$, |
$a_{1}=35114.4$, and |
| $b_{1}=36412.2$ |  |

The functionals given by Eq. (68) retrieve their nonrelativistic counterparts.
Although, in nonrelativistic NOFT approximations the $f(n_{I},n_{J})$
functions only affect the ‘same-spin’ terms (i.e. ${}_{X}{}^{2}D_{ij}^{ji}$
and ${}_{X}{}^{2}D_{\bar{i}\bar{j}}^{\bar{j}\bar{i}}$); therefore, the
‘opposite-spin’ ones (i.e. ${}_{X}{}^{2}D_{i\bar{j}}^{\bar{j}i}$ and
${}_{X}{}^{2}D_{\bar{i}j}^{j\bar{i}}$) could be defined as
$-\frac{n_{i}n_{j}}{2}$ (or using any other function), which may lead to an
extended definition of $f(n_{I},n_{J})$-functional approximations in the
relativistic scenario.

Let us briefly introduce the origin of the functionals presented in Table 2
that have their origin in the nonrelativistic context. The Müller, Buijse and
Baerends (MBB) functional was introduced independently by Müller and by Buijse
and Baerends [101, 102, 103], it produces an approximated functional which
fulfills the sum rule, (i.e. $\textrm{Tr}\left[{}^{2}{\bf
D}\right]=\frac{N_{e}(N_{e}-1)}{2}$). 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 $\alpha$ parameter was first proven
to have to be $\alpha\geq 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 solid state
physics and, according to Sharma et al. [105], the MBB functional
overcorrelates the electrons and the Power functional mediates with the
overcorrelation of the electrons by using a parameter $\alpha>1/2$. In some
applications, the bound $\alpha\geq 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 approximations developed to remove the so-call
self-interaction [107, 108] terms (i.e. they remove the nonphysical
${}^{2}D_{ii}^{ii}$ and ${}^{2}D_{\bar{i}\bar{i}}^{\bar{i}\bar{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 Padé
approximant expression. The MLSIC parameters were optimized to reproduce the
G2 test.

The energy expression for the $f(n_{I},n_{J})$ functionals reads as

$E^{\textrm{rel-}x}=2\sum_{i}h_{ii}n_{i}+\sum_{i,j}{2n_{i}n_{j}}J_{ij}-\sum_{i,j}f(n_{i},n_{j})\left[\left(K_{ij}-K^{G}_{ij}\right)+\left(L_{ij}-L^{G}_{ij}\right)\right]$
(69)

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

####  III.5.3 Approximations built imposing $N$-representability conditions

A second family of functionals approximations is based on the reconstruction
of the 2RDM elements imposing the D, Q, G $N$-representability conditions
[111, 83, 88, 98].

From Eq. (III.5), 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

${}^{2}{\bf D}=\begin{pmatrix}\left({}^{2}{\bf D}^{ij}_{ij},{}^{2}{\bf
D}^{ij}_{ji}\right)&0&0&0\\\ 0&\left({}^{2}{\bf
D}^{\bar{i}i}_{\bar{j}j},{}^{2}{\bf
D}^{\bar{i}j}_{\bar{i}j}\biggr{\rvert}_{i\neq j}\right)&\left({}^{2}{\bf
D}^{\bar{i}i}_{j\bar{j}},{}^{2}{\bf
D}^{\bar{i}j}_{j\bar{i}}\biggr{\rvert}_{i\neq j}\right)&0\\\
0&\left({}^{2}{\bf D}^{i\bar{i}}_{\bar{j}j},{}^{2}{\bf
D}^{j\bar{i}}_{\bar{i}j}\biggr{\rvert}_{i\neq j}\right)&\left({}^{2}{\bf
D}^{i\bar{i}}_{j\bar{j}},{}^{2}{\bf
D}^{j\bar{i}}_{j\bar{i}}\biggr{\rvert}_{i\neq j}\right)&0\\\
0&0&0&\left({}^{2}{\bf D}^{\bar{i}\bar{j}}_{\bar{i}\bar{j}},{}^{2}{\bf
D}^{\bar{i}\bar{j}}_{\bar{j}\bar{i}}\right)\end{pmatrix}.$ (70)

We may identify three decoupled blocks containing different spinor
combinations. The $\left({}^{2}{\bf D}^{ij}_{ij},{}^{2}{\bf
D}^{ij}_{ji}\right)$ and $\left({}^{2}{\bf
D}^{\bar{i}\bar{j}}_{\bar{i}\bar{j}},{}^{2}{\bf
D}^{\bar{i}\bar{j}}_{\bar{j}\bar{i}}\right)$ blocks resemble the
$\alpha\alpha$ and $\beta\beta$ blocks of the nonrelativistic case, while the
large middle block corresponds to the $\alpha\beta\alpha\beta$,
$\alpha\beta\beta\alpha$, 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}{\bf
D}_{j\bar{j}}^{\bar{i}i}=0$). Proposing functional approximations is
facilitated by the introduction of the so-called cumulant matrix, which can be
defined as

${\bm{\Lambda}}={}^{2}{\bf D}-{}^{2}{\bf D}^{\textrm{SD}}.$ (71)

Then, the auxiliar matrices ${\bm{\Delta}}$ and ${\bm{\Pi}}$ 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 $\Omega$ into subspaces mutually disjoint $\Omega_{p}$ (see Fig. 1 for
more details) we arrive to the relativistic version of the GNOF/PNOF$x$
($x=5,7$) functional approximations [112, 95, 90] that can be written as

$E^{\textrm{rel-GNOF/PNOF}x}=\sum_{p=1}^{N_{e}/2}E_{p}+\sum_{q\neq
p}^{N_{e}/2}E_{qp}.$ (72)

The first sum accounts for all intra-subspace contributions and reads as

$\displaystyle E_{p}={}$
$\displaystyle\sum_{i\in\Omega_{p}}n_{i}(2h_{ii}+J_{ii}+J^{G}_{ii}+L^{G}_{ii})$
$\displaystyle{}+\sum_{\begin{subarray}{c}i,j\in\Omega_{p}\\\ i\neq
j\end{subarray}}\Pi^{\textrm{intra}}_{i,j}(K_{ij}+K^{G}_{ij}+L_{ij}+L^{G}_{ij}),$
(73)

where

$\Pi^{\textrm{intra}}_{i,j}=\begin{cases}-\sqrt{n_{i}n_{j}},&i\;\textrm{or}\;j\leq
N_{e}/2\\\ +\sqrt{n_{i}n_{j}},&i,j>N_{e}/2,\\\ \end{cases}$ (74)

the second sum accounts for inter-subspace contributions ($E_{ba}$) that can
be defined as

$\displaystyle E_{qp}={}$
$\displaystyle\\!\\!\sum_{i\in\Omega_{q}}\sum_{j\in\Omega_{p}}\\!n_{i}n_{j}\left[2J_{ij}-(K_{ij}-K^{G}_{ij})-(L_{ij}-L^{G}_{ij})\right]$
$\displaystyle{}+\sum_{i\in\Omega_{q}}\sum_{j\in\Omega_{p}}\Pi^{\textrm{inter}}_{i,j;p,q}(K_{ij}+K^{G}_{ij}+L_{ij}+L^{G}_{ij}),$
(75)

where $\Pi^{\textrm{inter}}_{i,j;p,q}=0$ in rel-PNOF5,
$\Pi^{\textrm{inter}}_{i,j;p,q}=-\sqrt{n_{i}h_{i}n_{j}h_{j}}$ in rel-PNOF7
[113], $\Pi^{\textrm{inter}}_{i,j;p,q}=-4n_{i}h_{i}n_{j}h_{j}$ in rel-PNOF7s
131313A variant proposed for PNOF7 functional that is used to define NOF-MP2
method [90]. 141414The 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], and

$\Pi^{\textrm{inter}}_{i,j;p,q}\begin{cases}n_{i}^{d}n_{j}^{d}-\sqrt{n_{i}h_{i}n_{j}h_{j}}-\sqrt{n_{i}^{d}n_{j}^{d}},&i>q,j=p\;\textrm{or}\;i=q,j>p\\\
n_{i}^{d}n_{j}^{d}-\sqrt{n_{i}h_{i}n_{j}h_{j}}+\sqrt{n_{i}^{d}n_{j}^{d}},&i>q,j>p\\\
0,&i=q,j=p\\\ \end{cases}$ (76)

with $n_{i}^{d}=\frac{n_{i}h_{p}^{d}}{h_{p}}$, $h_{p}=1-n_{p}$, and
$h_{p}^{d}=h_{p}\exp{\left[-(h_{p}/(0.02\sqrt{2}))^{2}\right]}$ in GNOF [116].
Each electron pair is described by an even number of PS forming Kramers’
pairs. Moreover, the rel-GNOF/PNOF$x$ functionals presented in this work
produce 2-RDM elements that are real, but complex elements could be introduced
by taking other functions to approximate the $\Pi^{\textrm{intra}}$ and
$\Pi^{\textrm{inter}}$ sub-matrices.

Figure 1: Partition of the spinor space $\Omega$ into subspaces mutually
disjoint $\Omega_{p}$. While arrows refer to $\alpha$ or $\beta$ spin in the
nonrelativistic context, in this work they correspond to electrons occupying
spinors that are Kramers’ pairs. Let us remark that some NOs may remain
uncoupled (a.k.a. virtual spinors) like the orange ones in this scheme.

As we did for the $f(n_{I},n_{J})$-functionals, we briefly introduce the
origin of the nonrelativistic PNOF$x$ approximations whose relativistic
adaptations are presented in this work for completeness. PNOF5 was developed
to produce the correct electron distribution upon bond cleavage
processes[112]. This functional was initially introduced in the perfect-
pairing approach, where each sub-space $\Omega_{p}$ contained only two NOs.
Although, it was later extended beyond the perfect-pairing approach to allow
more orbitals in each sub-space. The lack of inter-subspace interaction in
PNOF5 lead to the later versions of PNOF$x$ 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 [113, 34, 35]. 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
correlation energy. 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 (III.5.1) in this limit.

###  III.6 Properties of rel-GNOF/PNOF$x$ functionals.

For ‘spin-compensated’ two-electron systems (see Appendix E) the rel-
GNOF/PNOF$x$ functionals introduced in this work reduce to

$\displaystyle E^{\textrm{rel-GNOF/PNOF}x}_{2e-}$
$\displaystyle=\sum_{i}n_{i}(2h_{ii}+J_{ii}+J_{ii}^{G}+L_{ii}^{G}){}-2\sum_{\begin{subarray}{c}i\\\
i\neq 1\end{subarray}}\sqrt{n_{1}n_{i}}(K_{1i}+K_{1i}^{G}+L_{1i}+L_{1i}^{G})$
$\displaystyle+\sum_{\mathclap{\begin{subarray}{c}i,j\\\ i\neq j,i\neq 1,j\neq
1\end{subarray}}}\sqrt{n_{i}n_{j}}(K_{ij}+K_{ij}^{G}+L_{ij}+L_{ij}^{G}),$ (77)

where 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 be proven to be the relativistic extension of the Fixed-Phases
functional [117], which is based on the NO representation of 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 the energy 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

$\Psi_{+}^{\textrm{APSG}}({\bf r}_{1},{\bf r}_{2},...,{\bf
r}_{N})=\widehat{\mathcal{A}}\prod^{N_{e}/2}_{p=1}\Xi_{p}({\bf r}_{2p-1},{\bf
r}_{2p}),$ (78)

where

$\Xi_{p}({\bf r},{\bf
r}^{\prime})=\frac{1}{\sqrt{2}}\sum_{i\in\Omega_{p}}c_{i}\left[\widetilde{\bm{\chi}}_{i}({\bf
r})\otimes\widetilde{\bm{\chi}}_{\bar{i}}({\bf
r}^{\prime})-\widetilde{\bm{\chi}}_{i}({\bf
r}^{\prime})\otimes\widetilde{\bm{\chi}}_{\bar{i}}({\bf r})\right]$ (79)

is a geminal wavefunction written in the NO basis including phase factors
(i.e. $c_{i}=\sqrt{n_{i}}\textrm{e}^{\textrm{i}\zeta_{i}}$, see Appendix E for
more details), and the $\widehat{\mathcal{A}}$ operator produces the
antisymmetrized product of geminal wavefunctions. The geminal wavefunctions
are normalized $\int d{\bf r}d{\bf r}^{\prime}\Xi^{\dagger}_{p}({\bf r},{\bf
r}^{\prime})\Xi_{p}({\bf r},{\bf r}^{\prime})=1$ for each electron pair, where

$\Xi^{\dagger}_{p}({\bf r},{\bf
r}^{\prime})=\frac{1}{\sqrt{2}}\sum_{i\in\Omega_{p}}c_{i}^{*}\left[\widetilde{\bm{\chi}}^{\dagger}_{i}({\bf
r})\otimes\widetilde{\bm{\chi}}^{\dagger}_{\bar{i}}({\bf
r}^{\prime})-\widetilde{\bm{\chi}}^{\dagger}_{i}({\bf
r}^{\prime})\otimes\widetilde{\bm{\chi}}^{\dagger}_{\bar{i}}({\bf r})\right].$
(80)

Imposing the strong-orthogonality condition

$\forall_{pq}\;\int d{\bf r}_{1}d{\bf r}_{2}\Xi^{\dagger}_{p}({\bf r}_{1},{\bf
r}_{2})\Xi_{q}({\bf r}_{1},{\bf r}_{2})=\delta_{pq},$ (81)

we arrive to the following expression for the 1-RDM

${\bf n}_{1}^{+,\textrm{APSG}}({\bf r},{\bf
r}^{\prime})=\sum_{p=1}^{N_{e}/2}\sum_{i\in\Omega_{p}}n_{i}\left[\widetilde{\bm{\chi}}_{i}({\bf
r})\widetilde{\bm{\chi}}^{\dagger}_{i}({\bf
r}^{\prime})+\widetilde{\bm{\chi}}_{\bar{i}}({\bf
r})\widetilde{\bm{\chi}}^{\dagger}_{\bar{i}}({\bf r}^{\prime})\right]$ (82)

Thus, the relativistic APSG wavefunction permits us to write the energy
functional

$\displaystyle
E^{\textrm{APSG}}\left[\\{\zeta_{i}\\},\\{n_{i}\\},\\{\widetilde{\bm{\chi}}_{i}\\}\right]$
$\displaystyle=2\sum_{i}n_{i}h_{ii}+\sum^{N_{e}/2}_{p\neq
q}\sum_{\begin{subarray}{c}i,j\\\
i\in\Omega_{p},j\in\Omega_{q}\end{subarray}}n_{i}n_{j}\left[2J_{ij}-(K_{ij}-K^{G}_{ij})-(L_{ij}-L^{G}_{ij})\right]$
(83) $\displaystyle+\sum_{p=1}^{N_{e}/2}\left[\sum_{{i\neq
j\in\Omega_{p}}}\textrm{e}^{\textrm{i}(\zeta_{j}-\zeta_{i})}\sqrt{n_{i}n_{j}}(K_{ij}+K^{G}_{ij}+L_{ij}+L^{G}_{ij})\right]$

which can be optimized subject to the conditions: a)
$\forall_{i}\;n_{i}\in[0,1]$, and b) $N_{e}=2\sum_{i}n_{i}$. 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 phases
$\textrm{e}^{\textrm{i}\zeta_{i}}=\pm 1$ are fixed (see Eq. (74)) while in
$E^{\textrm{APSG}}$ they are parameters, we can readily conclude that in
$E^{\textrm{APSG}}$ there is more flexibility during the optimization
procedure that can lead to a lower energy.

##  IV Closing remarks

In 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 (${\bf
n}_{1}({\bf r},{\bf r}^{\prime})$) is the charge $Q$; the minimum of the
energy ($E_{Q}$) is attained for the CI ground state wavefunction
($|\Psi_{Q}\rangle$). 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. this vacuum,
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 be
valid 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 $|\widetilde{0}_{v}\rangle$ and how 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. At this level, two 1-RDMs play a crucial role (i.e. ${\bf
n}^{\textrm{+}}_{1}({\bf r},{\bf r}^{\prime})$ and $\widetilde{\bf
n}^{\textrm{vp}}_{1}({\bf r},{\bf r}^{\prime})$). 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 because all vp effects as well
as electron-vacuum interactions are explicit functionals of ${\bf
n}^{\textrm{+}}_{1}({\bf r},{\bf r}^{\prime})$ and $\widetilde{\bf
n}^{\textrm{vp}}_{1}({\bf r},{\bf r}^{\prime})$. Only the functional of the
electron-electron interaction in terms of ${\bf n}^{\textrm{+}}_{1}({\bf
r},{\bf r}^{\prime})$ 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 the concept of an
$N$-electron wavefunction is recovered. When the Hamiltonian preserves time-
reversal symmetry it is possible 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, 120, 121, 122, 123, 124, 125, 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 are
particularly important in systems that are not ‘spin-compensated’ (i.e. that
in our case do not preserve time-reversal symmetry [128, 129, 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 in strongly-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
Skłodowska-Curie Individual Fellowship (891647-ReReDMFT). M.R.-M. would also
like to acknowledge Prof. M. Piris for the fruitful discussions.

##  Appendix A Commutation relation between the Hamiltonian and the opposite
charge operator

In order to prove the commutation relation between $\widehat{H}^{v}_{0}$ and
the opposite charge operator $\widehat{Q}$ let us first rewrite the
Hamiltonian operator in the 4-component spinor basis as

$\displaystyle\widehat{H}^{v}_{0}$
$\displaystyle=\sum_{I,J}{\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{J}\mathcal{H}_{IJ}+\sum_{I}\sum_{R}{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R}\mathcal{H}_{IR}$
(84)
$\displaystyle=\sum_{R,I}{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I}\mathcal{H}_{RI}-\sum_{R,S}{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{R}\mathcal{H}_{RS}$

where

$\mathcal{H}_{AB}=\int d{\bf r}d{\bf r}^{\prime}{\bm{\psi}}^{\dagger}_{A}({\bf
r}^{\prime})\left[\delta({\bf r}-{\bf r}^{\prime})\widehat{\bf T}_{D}+{\bf
v}_{\textrm{ext}}^{\textrm{nl}}({\bf r}^{\prime},{\bf
r})\right]{\bm{\psi}}_{B}({\bf r}).$ (85)

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

$\displaystyle\sum_{I,J,K}\mathcal{H}_{IJ}\left[{\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{J},{\widehat{b}^{\;\dagger}}_{K}{\widehat{b}^{\vphantom{\dagger}}}_{K}\right]-\sum_{I,J}\mathcal{H}_{IJ}\sum_{R}\left[{\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{J},{\widehat{d}^{\;\dagger}}_{R}{\widehat{d}^{\vphantom{\dagger}}}_{R}\right]$
(86)
$\displaystyle=\sum_{I,J,K}\mathcal{H}_{IJ}\left({\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{J}{\widehat{b}^{\;\dagger}}_{K}{\widehat{b}^{\vphantom{\dagger}}}_{K}-{\widehat{b}^{\;\dagger}}_{K}{\widehat{b}^{\vphantom{\dagger}}}_{K}{\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{J}\right)$
$\displaystyle=\sum_{I,J,K}\mathcal{H}_{IJ}\left({\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{J}{\widehat{b}^{\;\dagger}}_{K}{\widehat{b}^{\vphantom{\dagger}}}_{K}-{\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{J}{\widehat{b}^{\;\dagger}}_{K}{\widehat{b}^{\vphantom{\dagger}}}_{K}+{\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{K}\delta_{JK}-{\widehat{b}^{\;\dagger}}_{K}{\widehat{b}^{\vphantom{\dagger}}}_{J}\delta_{IK}\right)$
$\displaystyle=0$

the commutator
$\left[{\widehat{b}^{\;\dagger}}_{I}{\widehat{b}^{\vphantom{\dagger}}}_{J},{\widehat{d}^{\;\dagger}}_{R}{\widehat{d}^{\vphantom{\dagger}}}_{R}\right]=0$
because positronic and electronic creation and annihilation operators do
commute. The proof that the fourth term of the Hamiltonian (i.e.
$-\sum_{R,S}{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{R}\mathcal{H}_{RS}$)
also commutes with the opposite charge operator follows closely the one
presented for the first term; we therefore omit its detailed description.
Next, the commutation relation between the second term of the Hamiltonian and
the opposite charge operator reads as

$\displaystyle\sum_{I,J}\sum_{R}\mathcal{H}_{IR}\left[{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R},{\widehat{b}^{\;\dagger}}_{J}{\widehat{b}^{\vphantom{\dagger}}}_{J}\right]-\sum_{I}\sum_{R,S}\mathcal{H}_{IR}\left[{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R},{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{S}\right]$
(87)
$\displaystyle\sum_{I,J}\sum_{R}\mathcal{H}_{IR}\left({\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R}{\widehat{b}^{\;\dagger}}_{J}{\widehat{b}^{\vphantom{\dagger}}}_{J}-{\widehat{b}^{\;\dagger}}_{J}{\widehat{b}^{\vphantom{\dagger}}}_{J}{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R}\right)-\sum_{I}\sum_{R,S}\mathcal{H}_{IR}\left({\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R}{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{S}-{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{S}{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R}\right)$
$\displaystyle\sum_{I,J}\sum_{R}\mathcal{H}_{IR}\left({\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R}{\widehat{b}^{\;\dagger}}_{J}{\widehat{b}^{\vphantom{\dagger}}}_{J}-{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R}{\widehat{b}^{\;\dagger}}_{J}{\widehat{b}^{\vphantom{\dagger}}}_{J}-{\widehat{b}^{\;\dagger}}_{J}{\widehat{d}^{\;\dagger}}_{R}\delta_{IJ}\right)-\sum_{I}\sum_{R,S}\mathcal{H}_{IR}\left({\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R}{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{S}-{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R}{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{S}-{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{S}\delta_{SR}\right)$
$\displaystyle=\sum_{I}\sum_{R}\mathcal{H}_{IR}\left(-{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R}+{\widehat{b}^{\;\dagger}}_{I}{\widehat{d}^{\;\dagger}}_{R}\right)=0,$

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

$\displaystyle\sum_{I,J}\sum_{R}\mathcal{H}_{RI}\left[{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I},{\widehat{b}^{\;\dagger}}_{J}{\widehat{b}^{\vphantom{\dagger}}}_{J}\right]-\sum_{I}\sum_{R,S}\mathcal{H}_{RI}\left[{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I},{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{S}\right]$
(88)
$\displaystyle\sum_{I,J}\sum_{R}\mathcal{H}_{RI}\left({\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I}{\widehat{b}^{\;\dagger}}_{J}{\widehat{b}^{\vphantom{\dagger}}}_{J}-{\widehat{b}^{\;\dagger}}_{J}{\widehat{b}^{\vphantom{\dagger}}}_{J}{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I}\right)-\sum_{I}\sum_{R,S}\mathcal{H}_{RI}\left({\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I}{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{S}-{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I}\right)$
$\displaystyle\sum_{I,J}\sum_{R}\mathcal{H}_{RI}\left({\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I}{\widehat{b}^{\;\dagger}}_{J}{\widehat{b}^{\vphantom{\dagger}}}_{J}-{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I}{\widehat{b}^{\;\dagger}}_{J}{\widehat{b}^{\vphantom{\dagger}}}_{J}+{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{J}\delta_{IJ}\right)-\sum_{I}\sum_{R,S}\mathcal{H}_{RI}\left({\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I}{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{S}-{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I}{\widehat{d}^{\;\dagger}}_{S}{\widehat{d}^{\vphantom{\dagger}}}_{S}+{\widehat{d}^{\vphantom{\dagger}}}_{S}{\widehat{b}^{\vphantom{\dagger}}}_{I}\delta_{RS}\right)$
$\displaystyle=\sum_{I}\sum_{R}\mathcal{H}_{RI}\left({\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I}-{\widehat{d}^{\vphantom{\dagger}}}_{R}{\widehat{b}^{\vphantom{\dagger}}}_{I}\right)=0,$

that completes the demonstration of the commutation relation between the
Hamiltonian and the opposite charge operator.

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

In 1975, Gilbert extended the Hohenberg and Kohn theorems [131, 132, 133] for
external nonlocal potentials. In this work 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
$\lvert\Psi\rangle\in\mathcal{H}_{Q}$ (see Eq. (II.4)) 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.

###### Theorem 1

Given a relativistic Hamiltonian
$\widehat{H}=\widehat{T}_{D}+\widehat{V}^{\textrm{nl}}_{\textrm{ext}}+\widehat{W}$,
whose nondegenerated ground state wavefunction reads as $\Psi_{Q}$ (for a
particular charge sector $Q$). There exists a one-to-one correspondence
between $\Psi_{Q}$ and its 1-RDM (${\bf n}_{1}({\bf r},{\bf r^{\prime}})$).

Proof. Assume that two-different nonlocal external potentials
$\widehat{V}^{1,\textrm{nl}}_{\textrm{ext}}$ and
$\widehat{V}^{2,\textrm{nl}}_{\textrm{ext}}$ that differ in more than a
constant lead to the nondegenerated ground states $\Psi^{1}_{Q}$ and
$\Psi^{2}_{Q}$ (for a charge sector $Q$). That is to say,
$\widehat{H}_{1}\Psi^{1}_{Q}=(\widehat{T}_{D}+\widehat{V}^{1,\textrm{nl}}_{\textrm{ext}}+\widehat{W})\Psi^{1}_{Q}=E_{1}\Psi^{1}_{Q}$
and $\widehat{H}_{2}\Psi^{2}_{Q}=E_{2}\Psi^{2}_{Q}$. The 1-RDMs of
$\Psi^{1}_{Q}$ and $\Psi^{2}_{Q}$ are ${\bf n}_{1}^{1}({\bf r},{\bf
r}^{\prime})$ and ${\bf n}_{1}^{2}({\bf r},{\bf r}^{\prime})$, respectively.
If the two Hamiltonians only differ in the nonlocal external potential (i.e.
$\widehat{H}_{2}-\widehat{H}_{1}=\widehat{V}^{2,\textrm{nl}}_{\textrm{ext}}-\widehat{V}^{1,\textrm{nl}}_{\textrm{ext}}$),
by Rayleigh-Ritz variational principle we can write

$E_{1}=\langle\Psi^{1}_{Q}|\widehat{H}_{1}|\Psi^{1}_{Q}\rangle<\langle\Psi^{2}_{Q}|\widehat{H}_{1}|\Psi^{2}_{Q}\rangle=E_{1}^{2}$
(89)

and

$E_{2}=\langle\Psi^{2}_{Q}|\widehat{H}_{2}|\Psi^{2}_{Q}\rangle<\langle\Psi^{1}_{Q}|\widehat{H}_{2}|\Psi^{1}_{Q}\rangle=E_{2}^{1},$
(90)

so that

$\Delta E=(E_{1}^{2}-E_{1})+(E_{2}^{1}-E_{2})>0$ (91)

and

$\displaystyle\Delta E={}$ $\displaystyle(E_{1}^{2}-E_{1})+(E_{2}^{1}-E_{2})$
$\displaystyle={}$
$\displaystyle\langle\Psi^{1}_{Q}|\widehat{H}_{2}-\widehat{H}_{1}|\Psi^{1}_{Q}\rangle+\langle\Psi^{2}_{Q}|\widehat{H}_{1}-\widehat{H}_{2}|\Psi^{2}_{Q}\rangle$
$\displaystyle={}$ $\displaystyle-\int d{\bf r}d{\bf
r}^{\prime}\textrm{Tr}\left[\left({\bf v}^{\textrm{nl}}_{2,\textrm{ext}}({\bf
r}^{\prime},{\bf r})-{\bf v}^{\textrm{nl}}_{1,\textrm{ext}}({\bf
r}^{\prime},{\bf r})\right)\right.$ $\displaystyle{}\times\left(\left.{\bf
n}^{2}_{1}({\bf r},{\bf r}^{\prime})-{\bf n}^{1}_{1}({\bf r},{\bf
r}^{\prime})\right)\right]>0$ (92)

Since we imposed that $\widehat{V}^{1,\textrm{nl}}_{\textrm{ext}}$ and
$\widehat{V}^{2,\textrm{nl}}_{\textrm{ext}}$ are different, the corresponding
1-RDMs (${\bf n}^{1}_{1}({\bf r},{\bf r^{\prime}})$ and ${\bf n}^{2}_{1}({\bf
r},{\bf r^{\prime}})$) must also be different. Therefore, there exists a one-
to-one correspondence between the 1-RDM and the wavefunction.

$\Psi_{Q}\leftrightarrow{\bf n}_{1}({\bf r},{\bf r}^{\prime}).\qquad\square$
(93)

Consequently, there exists a functional of the energy

$E_{Q}=E\left[{\bf n}_{1}\right]\qquad\textrm{for}\qquad{\bf n}_{1}({\bf
r},{\bf r}^{\prime})\in\mathcal{D}_{Q}.$ (94)

This functional retrieves its nonrelativistic counterpart for a fixed number
of electrons (i.e. $Q=N_{e}$, 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
$\lvert\Psi_{+}\rangle$ and ${\bf n}^{+}_{1}({\bf r},{\bf r}^{\prime})$ 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. (III.2) (or Eq. (57)), similar arguments to the above presented lead to
the desired one-to-one correspondence (i.e.
$\lvert\Psi_{+}\rangle\leftrightarrow{\bf n}^{+}_{1}({\bf r},{\bf
r}^{\prime})$). Finally, neglecting vacuum polarization effects, we may obtain
the energy functional

$E_{N_{e}}=E\left[{\bf n}^{+}_{1}\right]\qquad\textrm{for}\qquad{\bf
n}_{1}\in\mathcal{D}_{N_{e}}^{+}$ (95)

that allows us to introduce np-ReRDMFT and clearly retrieves the
nonrelativistic limit.

##  Appendix C The npvp energy functional

The $E^{\textrm{npvp}}$ functional of the np 1-RDM and the vp 1-RDM reads as

$\displaystyle E^{\textrm{npvp}}\left[{\bf n}_{1}^{\textrm{+}},{\widetilde{\bf
n}_{1}^{\textrm{vp}}}\right]$ $\displaystyle=\int d{\bf r}d{\bf
r}^{\prime}\textrm{Tr}\left[\left(\delta({\bf r}-{\bf r}^{\prime})\widehat{\bf
T}_{D}({\bf r})+{\bf v}_{\textrm{ext}}^{\textrm{nl}}({\bf r}^{\prime},{\bf
r})\right)\left({\bf n}_{1}^{\textrm{+}}({\bf r},{\bf
r}^{\prime})+{\widetilde{\bf n}_{1}^{\textrm{vp}}({\bf r},{\bf
r}^{\prime})}\right)\right]$ (96) $\displaystyle+\sum_{\mu,\tau}\int d{\bf
r}_{1}\left[\sum_{\nu,\eta}\int d{\bf r}_{2}W_{\mu,\nu,\tau,\eta}({\bf
r}_{1},{\bf r}_{2}){\widetilde{n}^{\textrm{vp}}_{1,\eta,\nu}({\bf r}_{2},{\bf
r}_{2})}\right]{\widetilde{n}}^{+}_{1,\tau,\mu}({\bf r}_{1},{\bf r}_{1})$
$\displaystyle-\sum_{\mu,\nu,\tau,\eta}\int\int d{\bf r}_{1}d{\bf
r}_{2}W_{\mu,\nu,\tau,\eta}({\bf r}_{1},{\bf
r}_{2}){\widetilde{n}^{\textrm{vp}}_{1,\eta,\mu}({\bf r}_{2},{\bf
r}_{1})}{\widetilde{n}}^{+}_{1,\tau,\nu}({\bf r}_{1},{\bf r}_{2})$
$\displaystyle+\frac{1}{2}\int d{\bf r}_{1}d{\bf r}_{2}\textrm{Tr}\left[{\bf
W}({\bf r}_{1},{\bf r}_{2}){\widetilde{{\bf n}}^{\textrm{vp}}_{2}({\bf
r}_{1},{\bf r}_{2})}\right]+{{\widetilde{W}}\left[{\bf
n}_{1}^{\textrm{+}}\right]},$

where the explicit form of ${{\widetilde{W}}\left[{\bf
n}_{1}^{\textrm{+}}\right]}=\langle\Psi_{+}\left[{\bf
n}_{1}^{+}\right]|\widehat{\widetilde{W}}|\Psi_{+}\left[{\bf
n}_{1}^{+}\right]\rangle$ is unknown and needs to be approximated. Finally,
the usual np approximation corresponds to $E^{\textrm{np}}\left[{\bf
n}_{1}^{\textrm{+}}\right]=E^{\textrm{npvp}}\left[{\bf
n}_{1}^{\textrm{+}},{\bf 0}\right]$.

##  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}{\bf
D}\leftrightarrow\lvert\Psi_{+}\rangle$). 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 151515Notice that
the $\lambda_{\bar{i}\bar{j},\bar{k}\bar{l}}$ terms are obtained by replacing
the unbarred indices by barred ones.

$\displaystyle\lambda_{ij,kl}$
$\displaystyle=-\frac{\Delta_{ij}}{2}\delta_{ik}\delta_{jl}+\frac{\Delta_{ij}}{2}\delta_{il}\delta_{jk},$
(97a) and for the middle block $\displaystyle\lambda_{\bar{i}j,\bar{k}l}$
$\displaystyle=-\frac{\Delta_{\bar{i}j}}{2}\delta_{ik}\delta_{jl}+\frac{\Pi_{\bar{i}i,\bar{k}k}}{2}\delta_{ij}\delta_{kl},$
(97b) $\displaystyle\lambda_{i\bar{j},k\bar{l}}$
$\displaystyle=-\frac{\Delta_{i\bar{j}}}{2}\delta_{ik}\delta_{jl}+\frac{\Pi_{i\bar{i},k\bar{k}}}{2}\delta_{ij}\delta_{kl},$
(97c) $\displaystyle\lambda_{i\bar{j},\bar{k}l}$
$\displaystyle=\frac{\Delta_{i\bar{j}}}{2}\delta_{il}\delta_{jk}-\frac{\Pi_{i\bar{i},\bar{k}k}}{2}\delta_{ij}\delta_{kl},$
(97d) and $\displaystyle\lambda_{\bar{i}j,k\bar{l}}$
$\displaystyle=\frac{\Delta_{\bar{i}j}}{2}\delta_{il}\delta_{jk}-\frac{\Pi_{\bar{i}i,k\bar{k}}}{2}\delta_{ij}\delta_{kl}.$
(97e)

In the $i\neq j$ case, the ${\Delta}$ sub-matrices with a fixed $i$ and $j$
values must be Hermitian; thence, $\Delta_{i\bar{j}}^{*}=\Delta_{\bar{j}i}$.
From the antisymmetry properties of the 2-RDM, we obtain that
$\Delta_{\bar{j}i}=\Delta_{i\bar{j}}$, which makes these matrix elements be
real. Thus, the ${\bm{\Delta}}$ sub-matrices must be symmetric. The
hermiticity of the 2-RDM makes us to impose
$\Pi_{\bar{i}i,k\bar{k}}=\Pi_{k\bar{k},\bar{i}i}^{*}$. And, the antisymmetry
of the 2-RDM imposes
$\Pi_{\bar{i}i,k\bar{k}}=\Pi_{\bar{i}i,\bar{k}k}=\Pi_{i\bar{i},\bar{k}k}=\Pi_{i\bar{i},k\bar{k}}$.

Using the above definitions for the SD contribution and the cumulant matrix
elements, the 2-RDM elements can be approximated as

${}^{2}D_{ij}^{kl}$
$\displaystyle=\frac{n_{i}n_{j}-\Delta_{ij}}{2}\delta_{ik}\delta_{jl}-\frac{n_{i}n_{j}-\Delta_{ij}}{2}\delta_{il}\delta_{jk},$
(98a) ${}^{2}D_{\bar{i}\bar{j}}^{\bar{k}\bar{l}}$
$\displaystyle=\frac{n_{\bar{i}}n_{\bar{j}}-\Delta_{\bar{i}\bar{j}}}{2}\delta_{{i}{k}}\delta_{{j}{l}}-\frac{n_{\bar{i}}n_{\bar{j}}-\Delta_{\bar{i}\bar{j}}}{2}\delta_{{i}{l}}\delta_{{j}{k}},$
(98b) ${}^{2}D_{i\bar{j}}^{k\bar{l}}$
$\displaystyle=\frac{n_{i}n_{\bar{j}}-\Delta_{i\bar{j}}}{2}\delta_{ik}\delta_{{j}{l}}+\frac{\Pi_{i\bar{i},k\bar{k}}}{2}\delta_{ij}\delta_{kl},$
(98c) ${}^{2}D_{\bar{j}i}^{\bar{l}k}$
$\displaystyle=\frac{n_{i}n_{\bar{j}}-\Delta_{\bar{j}i}}{2}\delta_{ik}\delta_{jl}+\frac{\Pi_{\bar{i}i,\bar{k}k}}{2}\delta_{ij}\delta_{kl},$
(98d) ${}^{2}D_{\bar{j}i}^{k\bar{l}}$
$\displaystyle=-\frac{n_{i}n_{\bar{j}}-\Delta_{\bar{j}i}}{2}\delta_{ik}\delta_{jl}-\frac{\Pi_{\bar{i}i,k\bar{k}}}{2}\delta_{ij}\delta_{kl},$
(98e) and ${}^{2}D_{i\bar{j}}^{\bar{l}k}$
$\displaystyle=-\frac{n_{i}n_{\bar{j}}-\Delta_{i\bar{j}}}{2}\delta_{ik}\delta_{jl}-\frac{\Pi_{i\bar{i},\bar{k}k}}{2}\delta_{ij}\delta_{kl}.$
(98f)

The antisymmetry properties of the 2-RDM elements implies that Eqs. (98a),
(98b), and (98c) provide enough information to build the rest of 2-RDM matrix
elements. In contrast to the nonrelativistic approach, the 2-RDM matrix
elements given by Eqs. (98e) and (98f) contribute to the energy; we must
approximate them in terms of the ${\bm{\Delta}}$ and ${{\bm{\Pi}}}$ 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, 98 for more details)

$D_{IJ}^{KL}=\frac{1}{2}\langle\Psi_{+}|{\widehat{\widetilde{b}}^{\;\dagger}}_{I}{\widehat{\widetilde{b}}^{\;\dagger}}_{I}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{L}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{K}|\Psi_{+}\rangle\,,$
(99)
$Q_{IJ}^{KL}=\frac{1}{2}\langle\Psi_{+}|{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{I}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{J}{\widehat{\widetilde{b}}^{\;\dagger}}_{L}{\widehat{\widetilde{b}}^{\;\dagger}}_{K}|\Psi_{+}\rangle\,,$
(100)
$G_{IJ}^{KL}=\frac{1}{2}\langle\Psi_{+}|{\widehat{\widetilde{b}}^{\;\dagger}}_{I}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{J}{\widehat{\widetilde{b}}^{\;\dagger}}_{L}{\widehat{\widetilde{b}}^{\vphantom{\dagger}}}_{K}|\Psi_{+}\rangle\,.$
(101)

In the nonrelativistic approach, these matrices can be split in terms of
different spin-blocks for the evaluation of the positive semi-definite
character.

  * •

The D-Condition:

Eq. (99) tests the positive semidefinite character of the 2-RDM. From the
block structure of the 2-RDM (see Eq. (70)), let us first focus on the block
$({}^{2}{\bf D}^{ij}_{ij},{}^{2}{\bf D}^{ij}_{ji})$, which using Eq. (98a)
gives us the condition that for $i\neq j$ there is a set of
$[(M/2)\times(M/2-1)]/2$ two-by-two sub-blocks of the form161616For $i=j$, the
elements are 0 by definition and they fulfill the condition.

$\begin{pmatrix}\frac{(n_{i}n_{j}-\Delta_{ij})}{2}&-\frac{(n_{i}n_{j}-\Delta_{ij})}{2}\\\
-\frac{(n_{i}n_{j}-\Delta_{ij})}{2}&\frac{(n_{i}n_{j}-\Delta_{ij})}{2}\\\
\end{pmatrix},$ (102)

which upon diagonalization produce the eigenvalues 0 and
$2(n_{i}n_{j}-\Delta_{ij})$. Since the ${\bf D}$ matrix must be positive
semidefinite, we can conclude that

$\Delta_{ij}\leq n_{i}n_{j}.$ (103)

Similarly, the third block produces the same type of eigenvalues and
conditions (i.e. $\Delta_{\bar{i}\bar{j}}\leq n_{\bar{i}}n_{\bar{j}}$).
Finally, for the middle block and $i\neq j$ we have a set of
$[(M/2)\times(M/2)]$ two-by-two blocks of the form

$\begin{pmatrix}\frac{(n_{i}n_{\bar{j}}-\Delta_{i{\bar{j}}})}{2}&-\frac{(n_{i}n_{\bar{j}}-\Delta_{i{\bar{j}}})}{2}\\\
-\frac{(n_{i}n_{\bar{j}}-\Delta_{i{\bar{j}}})}{2}&\frac{(n_{i}n_{\bar{j}}-\Delta_{i{\bar{j}}})}{2}\\\
\end{pmatrix},$ (104)

whose eigenvalues are 0 and $2(n_{i}n_{\bar{j}}-\Delta_{i{\bar{j}}})$, and to
fulfill the D-condition the following constraint must hold

$\Delta_{i{\bar{j}}}\leq n_{i}n_{\bar{j}}.$ (105)

Also, for the middle block we have one large $[(M/2)\times(M/2)]$ block whose
off-diagonal elements are given by all the ${\bm{\Pi}}$ matrix off-diagonal
elements. The diagonal part of this block contains the $i=j$ case, which
includes the ${\bm{\Delta}}$ 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.

  * •

The Q-Condition:

Using the anticommutation rules for the creation and annihilation operators,
we may rewrite Eq. (100) as

$\displaystyle Q_{IJ}^{KL}$
$\displaystyle={}^{2}D_{KL}^{IJ}+\frac{1}{2}\left[\delta_{JL}(\delta_{IK}-{}^{1}D_{K}^{I})-\delta_{IL}(\delta_{JK}-{}^{1}D_{K}^{J})+\delta_{JK}{}^{1}D_{L}^{I}-\delta_{IK}{}^{1}D_{L}^{J}\right]$
(106)
$\displaystyle={}^{2}D_{KL}^{IJ}+\frac{1}{2}\left[\delta_{JL}\delta_{IK}h_{K}-\delta_{IL}\delta_{JK}h_{K}+\delta_{JK}\delta_{LI}h_{L}-\delta_{IK}\delta_{JL}h_{L}\right],$
(107)

where we have defined $h_{I}=1-n_{I}$, and taken the 1-RDM in the NO
representation in the last expression. It is straightforward to recognize that
the ${\bm{Q}}$ matrix shows the same block structure as the ${\bm{D}}$ matrix.
Hence, we have one large block that depends on the ${\bm{\Pi}}$ and
${\bm{\Delta}}$ matrices (and on the basis set), and two-by-two blocks for
$i\neq j$ of the form

$\begin{pmatrix}\frac{(1-n_{i}-n_{j}+n_{i}n_{j}-\Delta_{ij})}{2}&-\frac{(1-n_{i}-n_{j}+n_{i}n_{j}-\Delta_{ij})}{2}\\\
-\frac{(1-n_{i}-n_{j}+n_{i}n_{j}-\Delta_{ij})}{2}&\frac{(1-n_{i}-n_{j}+n_{i}n_{j}-\Delta_{ij})}{2}\\\
\end{pmatrix},$ (108)
$\begin{pmatrix}\frac{(1-n_{i}-n_{\bar{j}}+n_{i}n_{\bar{j}}-\Delta_{i{\bar{j}}})}{2}&-\frac{(1-n_{i}-n_{\bar{j}}+n_{i}n_{\bar{j}}-\Delta_{i{\bar{j}}})}{2}\\\
-\frac{(1-n_{i}-n_{\bar{j}}+n_{i}n_{\bar{j}}-\Delta_{i{\bar{j}}})}{2}&\frac{(1-n_{i}-n_{\bar{j}}+n_{i}n_{\bar{j}}-\Delta_{i{\bar{j}}})}{2}\\\
\end{pmatrix},$ (109)

that upon diagonalization lead us to the following conditions

$\Delta_{ij}\leq h_{i}h_{j}\qquad\textrm{and}\qquad\Delta_{i{\bar{j}}}\leq
h_{i}h_{\bar{j}}.$ (110)

  * •

The G-Condition:

As we did for the Q-condition, we may rewrite Eq. (101), using the
anticommutation relations of the creation and annihilation operators, as

$\displaystyle G_{IJ}^{KL}$
$\displaystyle=-{}^{2}D_{IL}^{KJ}+\frac{1}{2}\delta_{JL}{}^{1}D_{I}^{K}$ (111)
$\displaystyle=-{}^{2}D_{IL}^{KJ}+\frac{1}{2}\delta_{JL}\delta_{IK}n_{I},$
(112)

where in the last expression we have used the 1-RDM in the NO representation.
For $i\neq j$ the first block $\left({}^{2}{\bf D}^{ij}_{ij},{}^{2}{\bf
D}^{ij}_{ji}\right)$ (and also for the third block) the $\bf G$ matrix is
formed by a large basis set dependent block that does not provide any extra
analytical constraint. But it also contains $1\times 1$ blocks with elements
of the form $-{}^{2}D_{ij}^{ij}+\frac{1}{2}n_{i}$ that lead to a condition of
the form

$\Delta_{ij}\geq-h_{i}n_{j},$ (113)

which is easy to satisfy with ONs between 0 and 1 and $\Delta_{ij}\geq 0$.
Using the matrix elements of the central block of the 2-RDM let us introduce
the ${\bf G}^{\textrm{central}}$ matrix that reads as

${\bf G}^{\textrm{central}}=\begin{pmatrix}\left({\bf
G}^{\bar{i}j}_{\bar{j}i},{\bf G}^{\bar{i}j}_{\bar{i}j}\biggr{\rvert}_{i\neq
j}\right)&\left({\bf G}^{\bar{i}\bar{j}}_{ji},{\bf
G}^{\bar{i}\bar{i}}_{jj}\biggr{\rvert}_{i\neq j}\right)\\\ \left({\bf
G}^{ij}_{\bar{j}\bar{i}},{\bf G}^{jj}_{\bar{i}\bar{i}}\biggr{\rvert}_{i\neq
j}\right)&\left({\bf G}^{i\bar{j}}_{j\bar{i}},{\bf
G}^{j\bar{i}}_{j\bar{i}}\biggr{\rvert}_{i\neq j}\right)\\\ \end{pmatrix},$
(114)

whose elements can be built using Eq. (111). It is easy to recognize that this
block can be split into three disjoint sub-blocks (i.e. $\left({\bf
G}^{\bar{i}j}_{\bar{j}i},{\bf G}^{\bar{i}j}_{\bar{i}j}\biggr{\rvert}_{i\neq
j}\right)$, $\left({\bf G}^{i\bar{j}}_{j\bar{i}},{\bf
G}^{j\bar{i}}_{j\bar{i}}\biggr{\rvert}_{i\neq j}\right)$, and the rest). The
latter sub-block is basis set dependent and does 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

$\frac{-(n_{i}n_{\bar{i}}-\Delta_{i\bar{i}})-\Pi_{\bar{i}i,\bar{i}i}+n_{\bar{i}}}{2},$
(115)

which fulfill the $N$-representability condition if we let for example
$\Pi_{\bar{i}i,\bar{i}i}=n_{\bar{i}}$ and
$\Delta_{i\bar{i}}=n_{i}n_{\bar{i}}$. Nevertheless, this sub-block also
produces for $i\neq j$ two-by-two sub-blocks of the form

$\begin{pmatrix}G^{\bar{i}j}_{\bar{i}j}&G^{\bar{i}j}_{\bar{j}i}\\\
G^{\bar{j}i}_{\bar{i}j}&G^{\bar{j}i}_{\bar{j}i}\\\
\end{pmatrix}=\frac{1}{2}\begin{pmatrix}\Delta_{\bar{i}j}+n_{\bar{i}}h_{j}&-\Pi_{\bar{j}j,\bar{i}i}\\\
-\Pi_{\bar{i}i,\bar{j}j}&\Delta_{\bar{j}i}+n_{\bar{j}}h_{i}\\\ \end{pmatrix},$
(116)

whose eigenvalues obtained upon diagonalization read as

$\displaystyle\frac{1}{4}\left(n_{\bar{i}}h_{j}+\Delta_{{\bar{i}}j}+\Delta_{{\bar{j}}i}+n_{\bar{j}}h_{i}\right)$
$\displaystyle\pm\frac{1}{4}\sqrt{4\Pi^{2}_{\bar{i}i,\bar{j}j}+(n_{\bar{i}}h_{j}-n_{\bar{j}}h_{i})^{2}};$
(117)

where we have used the symmetry properties of the ${\bm{\Pi}}$ matrix. The
eigenvalue with the positive square root is greater or equal to zero and it
satisfies the requirement. On the other hand, the eigenvalue with a negative
square root leads to some bounds for the elements of the ${\bm{\Pi}}$ matrix.
These bounds depend on the ${\bm{\Delta}}$ matrix and they can be written as

$\left(n_{\bar{i}}h_{j}+\Delta_{{\bar{i}}j}+\Delta_{{\bar{j}}i}+n_{\bar{j}}h_{i}\right)^{2}\geq
4\Pi^{2}_{\bar{i}i,\bar{j}j}+(n_{\bar{i}}h_{j}-n_{\bar{j}}h_{i})^{2},$ (118)
$n_{\bar{i}}h_{j}(2\Delta_{{\bar{i}}j}+2\Delta_{{\bar{j}}i}+4n_{\bar{j}}h_{i})+\Delta_{{\bar{i}}j}^{2}+\Delta_{{\bar{j}}i}^{2}+2\Delta_{{\bar{j}}i}n_{\bar{j}}h_{i}+\Delta_{{\bar{i}}j}(2\Delta_{{\bar{j}}i}+2n_{\bar{j}}h_{i})\geq
4\Pi^{2}_{\bar{i}i,\bar{j}j},$ (119)

which in the particular case when $\Delta_{{\bar{i}}j}=\Delta_{{\bar{j}}i}=0$
lead to

$\sqrt{n_{\bar{i}}h_{i}n_{\bar{j}}h_{j}}\geq\lvert\Pi_{\bar{i}i,\bar{j}j}\rvert.$
(120)

Table 3: Definition of the ${\bm{\Pi}}$ matrix elements, where
$n_{i}^{d}={n_{i}h_{p}^{d}/h_{p}}$ with $i\in\Omega_{p}$, $h_{p}=1-n_{p}$, and
$h_{p}^{d}=h_{p}\exp{\left[-(h_{p}/(0.02\sqrt{2}))^{2}\right]}$. functional |  $\Pi^{\textrm{intra}}_{ij}$; $i,j\in\Omega_{p}$, $i\neq j$ |  $\Pi^{\textrm{inter}}_{ij}$; $i\in\Omega_{p}$, $j\in\Omega_{q}$, $p\neq q$ | Ref.
---|---|---|---
PNOF5 | $-\sqrt{n_{i}n_{j}}$ |  $i=p$ or $j=p$ | 0 | all | 136
$\sqrt{n_{i}n_{j}}$ | otherwise
PNOF7 | $-\sqrt{n_{i}n_{j}}$ |  $i=p$ or $j=p$ | $-\sqrt{n_{i}n_{j}h_{i}h_{j}}$ | all | 95,113
$\sqrt{n_{i}n_{j}}$ | otherwise
PNOF7s | $-\sqrt{n_{i}n_{j}}$ |  $i=p$ or $j=p$ | $-4n_{i}n_{j}h_{i}h_{j}$ | all | 90
$\sqrt{n_{i}n_{j}}$ | otherwise
GNOF | $-\sqrt{n_{i}n_{j}}$ |  $i=p$ or $j=p$ | $n_{i}^{d}n_{j}^{d}-\sqrt{n_{i}n_{j}h_{i}h_{j}}-\sqrt{n_{i}^{d}n_{j}^{d}}$ |  $i>p$ and $j=q$ or $i=p$ and $j>q$ | 116
$\sqrt{n_{i}n_{j}}$ | otherwise | $n_{i}^{d}n_{j}^{d}-\sqrt{n_{i}n_{j}h_{i}h_{j}}+\sqrt{n_{i}^{d}n_{j}^{d}}$ | $i>p,j>q$
|  | 0 | otherwise

From this analysis, we conclude that at the np-ReNOFT level the ansatz
proposed in Eqs. (97a)-(97c) provides equivalent eigenvalues for the D-, Q-,
and G-conditions as the nonrelativistic context. Indeed, only the definition
of the auxiliar matrix ${\bm{\Pi}}$ 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) because now the elements of the form
${}^{2}D_{\bar{i}j}^{i\bar{j}}$ contribute to the energy. Nevertheless, when
the ${\bm{\Pi}}$ matrix is real and the ONs for barred and unbarred states are
the same (i.e. $n_{i}=n_{\bar{i}}$), only two indices are needed; thus, the
Eq. (120) can be rewritten as in the nonrelativistic approach

$\sqrt{n_{i}h_{i}n_{j}h_{j}}\geq\lvert\Pi_{i,j}\rvert,$ (121)

with $\Pi_{i,j}=\Pi_{\bar{i}i,\bar{j}j}$.

To ensure that the rel-GNOF/PNOF$x$ functionals used in this work recover
their nonrelativistic counterparts in the nonrelativistic limit, we impose the
following constraints:

  1. 1.

Partition the PS space into subspaces $\\{\Omega_{p}\\}$ (see Fig. 1 for more
details).

  2. 2.

Let $\Delta_{ij}=\Delta_{i\bar{j}}=\Delta_{\bar{i}j}=\Delta_{\bar{i}\bar{j}}$.

  3. 3.

Define the diagonal terms $\Pi_{i,i}=n_{i}$ and $\Delta_{ii}=n_{i}^{2}$ to
also satisfy Eq. (115).

  4. 4.

Set $\Delta_{ij}=n_{i}n_{j}$ if $i,j\in\Omega_{p}$ and $\Delta_{ij}=0$
otherwise (this selection satisfies the above inequalities).

And we arrive to the rel-GNOF/PNOF$x$ functionals presented in this work. The
only missing term is the definition of the ${\bm{\Pi}}$ matrix elements. In
Table 3 we have collected different definitions, which are also employed in
the nonrelativistic context; that (as we have seen) are also valid in the
relativistic np approximation scenario. Let us remark that the ${\bm{\Pi}}$
matrix is separated into two contributions. One contribution formed by indices
belonging to the same $\Omega_{p}$ subspace (named as
${\bm{\Pi}}^{\textrm{intra}}$) that accounts for intra-subspace interactions;
the second contribution where the indices belong to different $\Omega$
subspaces (denoted as ${\bm{\Pi}}^{\textrm{inter}}$). Using these terms in the
definitions of the 2-RDM elements (Eqs. 98a-98f), we arrive to the
approximated (reconstructed) 2-RDM matrix. Inserting this matrix elements in
Eq. (III.5) we arrive to the functionals presented in Eq. (72).

##  Appendix E The relativistic Fixed-Phases functional.

In 1956 Löwdin and Shull proved that the configuration interaction
coefficients become simple functions of the ONs (i.e. $c_{I}=\pm\sqrt{n_{I}}$)
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
171717Recall that the ‘spin’ information is contained in the 4-component
spinors.

$\Psi_{+}^{2e-}({\bf r}_{1},{\bf r}_{2})=\sum_{I,J}C_{IJ}{\bm{\psi}}_{I}({\bf
r}_{1})\otimes{\bm{\psi}}_{J}({\bf r}_{2}),$ (122)

where the coefficients $C_{IJ}\in\mathbb{C}$ and form an antisymmetric matrix
${\bf C}$ (i.e. $C_{IJ}=-C_{JI}$ and $C_{II}=0$). From the normalization of
the wavefunction we have that $\sum_{I,J}\lvert C_{IJ}\rvert^{2}=1$. Applying
the Carlson–Keller expansion [138] (a.k.a. the Schmidt decomposition [139]) to
the wavefunction (122) we arrive to

$\displaystyle\Psi_{+}^{2e-}({\bf r}_{1},{\bf r}_{2})={}$