# Variational minimization scheme for the one-particle reduced density matrix functional theory in the ensemble N-representability domain

Matthieu Vladaj, Quentin Marécat, Bruno Senjean, Matthieu Saubanère

► To cite this version:

Matthieu Vladaj, Quentin Marécat, Bruno Senjean, Matthieu Saubanère. Variational minimization scheme for the one-particle reduced density matrix functional theory in the ensemble N-representability domain. The Journal of Chemical Physics, 2024, 161 (7), 10.1063/5.0219898. hal-04672766

HAL Id: hal-04672766

https://hal.science/hal-04672766v1

Submitted on 19 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Variational minimization scheme for the one-particle reduced density matrix
functional theory in the ensemble N-representability domain

Matthieu Vladaj,¹ Quentin Marécat,¹ Bruno Senjean,¹ and Matthieu Saubanère¹,²,∗

¹ICGM, Université de Montpellier, CNRS, ENSCM, Montpellier, France

²Univ. Bordeaux, CNRS, LOMA, UMR 5798, F-33400 Talence, France

The one-particle reduced density-matrix (1-RDM) functional theory is a promising alternative to density-functional theory (DFT) that uses the 1-RDM rather than the electronic density as a basic variable. However, long-standing challenges such as the lack of Kohn-Sham scheme and the complexity of the pure *N*-representability conditions are still impeding its wild utilization. Fortunately, ensemble *N*-representability conditions derived in the natural orbital basis are known and trivial, such that almost every functionals of the 1-RDM are actually natural orbital functionals which do not perform well for all the correlation regimes. In this work, we propose a variational minimization scheme in the ensemble *N*-representable domain that is not restricted to the natural orbital representation of the 1-RDM. We show that splitting the minimization into the diagonal and off-diagonal part of the 1-RDM can open the way toward the development of functionals of the orbital occupations, which remains a challenge for the generalization of site-occupation functional theory in chemistry. Our approach is tested on the uniform Hubbard model using the Müller and the Töws-Pastor functionals, as well as on the dihydrogen molecule using the Müller functional.

# I. INTRODUCTION

Over the last two decades, density functional theory (DFT) has played a major role in the study of the physico-chemical properties of molecules and materials. The key idea behind DFT, formulated by Hohenberg and Kohn [1], is to replace the *N*-electron wavefunction by the electronic density as the fundamental variable of the many-electron problem. Within the Kohn-Sham (KS) decomposition scheme [2], the electronic repulsion effects are treated implicitly by a Hartree-exchange-correlation (Hxc) functional of the density, thus leading to a good ratio between computational cost and accuracy. However, despite the numerous Hxc functionals developed in the past decades [3, 4], KS-DFT still faces issues when dealing with systems that exhibit strong (or static) correlation effects [5, 6].

Alternatively, replacing the density by the one-particle reduced density-matrix (1-RDM) as the fundamental variable leads to the 1-RDM functional theory (1-RDMFT) [7–10] and has several advantages. First, kinetic and exchange energy contributions are expressed exactly and analytically with respect to the 1-RDM. Second, fractional occupations of the orbitals are naturally obtained within 1-RDMFT, thus allowing for the proper description of bond dissociation where static correlation effects are significant. Despite these advantages, 1-RDMFT does not yet compete with DFT, partially due to the lack of an efficient KS scheme. Indeed, the non-idempotent 1-RDM of an interacting system cannot be mapped to its non-interacting equivalent that is idempotent by construction, as it results from a single-electron Schrödinger equation [11]. Besides, the proper minimization of the energy functional of the 1-RDM is

highly non-trivial as it requires maintaining the pure *N*-representability conditions of the 1-RDM (also known as the generalized Pauli constraints [12–14]) which number increases significantly with the number of electrons and orbitals [15]. Instead, the 1-RDM is usually written in the natural orbital (NO) representation, and the minimization is performed via a Langrangian-based system of differential equations where both the NOs and their occupation numbers are optimized while satisfying the more trivial ensemble *N*-representability conditions [9, 16, 17]. Because these conditions are known in the NO basis, Hxc functionals are almost all NO-functionals (NOFs) [18–27]. The most recent and efficient ones are the Piris natural orbitals functionals (PNOFs) [28–38] based on the cumulants [39]. Note that despite promising results, NOFs do not perform well for all the correlation regimes [40–44], although recent progress are worth mentioning [45].

Solving the aforementioned non-linear differential equations is computationally demanding, so that alternative strategies have been developed to define efficient and robust algorithms to perform the energy functional minimization. Gilbert proposed to determine the NOs using a non-local potential for functionals that are explicitly dependent on the 1-RDM [8], whose existence was later shown by Pernal in the case of NOFs [46]. In this approach, partially occupied orbitals are all degenerate. A faster scheme for the NO optimization is based on a hermitian matrix built from the Lagrange multiplier, as proposed by Piris and Ugalde [33]. One can also use the Hessian and couple the optimization of the NOs and occupation numbers [14]. Another interesting scheme called the local 1-RDMFT was developed by Lathiotakis and coworkers, where the minimization of the NOs is restricted to a domain where the orbitals are eigenfunctions of a single-electron Hamiltonian with a local potential, thus bridging the 1-RDMFT and KS-DFT approaches [47, 48]. Written in the NO representation, these methods as well as methods based on an

<footnote>* matthieu.saubanere@cnrs.fr</footnote>
analogy with the temperature are particularly suited for NOFs [49–52]. Indeed, in that case the energy does only depend explicitly on the occupation numbers while the dependence on the NOs is implicitly accounted to in the Hartree and exchange integrals. Alternatively, minimization schemes based on a different representation than natural orbitals were also proposed, for instance using projected gradient-based algorithms constrained to preserve ensemble N-representability conditions [53].

In this work, we introduce a new variational minimization scheme so-called density-matrix interpolation variational ansatz (DIVA) which relies on remarkable topological properties of the ensemble N-representable set of 1-RDM. In contrast to the usual minimization schemes adapted to NOFs, DIVA aims to minimize the energy functional of the 1-RDM expressed in any orthonormal orbital basis, thus motivating the development of new 1-RDM functionals beyond NOFs. By separating the minimization over the diagonal and off-diagonal part of the 1-RDM, one can also make a connection between 1-RDMFT and DFT in a lattice, also called site-occupation functional theory (SOFT) [54–57]. We show how this separation can be used to estimate Hartree Hxc potentials, thus paving the way toward the development of functionals of the orbital occupations which is the current limitation for the implementation of SOFT in chemistry [58, 59].

The paper is organized as follows. In Sec. II, we give a brief introduction to 1-RDMFT. The separation of the minimization process into that bridges 1-RDMFT and SOFT is then presented in Sec. III. Our variational minimization strategy and its various implementations are described in Sec. IV, and tested on the Hubbard model and the hydrogen molecule in Sec. V. Finally, conclusions and perspectives are provided in Sec. VI. Note that in the following, we will use ‘1-RDM’ and ‘density matrix’ to refer to the same object, i.e. the one-particle reduced density matrix.

## II. REDUCED DENSITY-MATRIX FUNCTIONAL THEORY

Let us consider the second-quantized electronic Hamiltonian written in an orthonormal basis set, referred to as the computational basis set and composed of N spin-orbitals {$\phi_i(\mathbf{x})$} that are different from the set of natural spin-orbitals {$\phi_i(\mathbf{x})$},

$$ \begin{aligned} \hat{H} &= \hat{T} + \hat{U} \\ &= \sum_{ij}^{N} t_{ij} \hat{c}_i^\dagger \hat{c}_j + \frac{1}{2} \sum_{ijkl} U_{ijkl} \hat{c}_i^\dagger \hat{c}_j^\dagger \hat{c}_k \hat{c}_l, \end{aligned} \quad (1) $$

where $\hat{c}_i^\dagger$ ($\hat{c}_i$) corresponds to the creation (annihilation) of an electron in $\phi_i(\mathbf{x})$. The one-body integrals

$$ t_{ij} = \int d\mathbf{x} \varphi_i^*(\mathbf{x}) \left( -\frac{\nabla^2}{2} + \sum_I \frac{Z_I}{|\mathbf{r} - \mathbf{R}_I|} \right) \varphi_j(\mathbf{x}) \quad (2) $$

contain the kinetic and external potential contributions, while the two-body integrals

$$ U_{ijkl} = \int d\mathbf{x}_1 d\mathbf{x}_2 \frac{\varphi_i^*(\mathbf{x}_1) \varphi_j^*(\mathbf{x}_2) \varphi_k(\mathbf{x}_2) \varphi_l(\mathbf{x}_1)}{|\mathbf{r}_1 - \mathbf{r}_2|} \quad (3) $$

contain contributions from the electron-electron interactions. The summation in Eq. (2) runs over the number of atoms at positions $\mathbf{R}_I$ and $\mathbf{x}_i = (\mathbf{r}_i, \sigma_i)$ corresponds to the position $\mathbf{r}_i$ and spin $\sigma_i$ of the electron $i$.

In 1-RDMFT, the basic variable $\gamma$ is the 1-RDM given by

$$ \gamma_{ij} = \langle \Psi | \hat{c}_i^\dagger \hat{c}_j | \Psi \rangle, \quad (4) $$

with $|\Psi\rangle$ the ground-state wavefunction of $\hat{H}$. The energy can in principle be written as a functional of the 1-RDM,

$$ E[\gamma] = T[\gamma] + W[\gamma], \quad (5) $$

where

$$ T[\gamma] = \sum_{ij} t_{ij} \gamma_{ij} \quad (6) $$

is an explicit functional of the 1-RDM and $W[\gamma]$ is the energy functional describing the electronic repulsion. The analytical form of this functional remains unknown and is usually approximated, though it can be written using the Levy–Lieb formalism as follows,

$$ W[\gamma] = \min_{\Psi \to \gamma} \langle \Psi | \hat{U} | \Psi \rangle, \quad (7) $$

where the notation $\Psi \to \gamma$ refers to any wavefunction giving the one-particle density matrix $\gamma$. Within the variational principle of 1-RDMFT, the ground-state energy $E_0$ can be obtained as a saddle point,

$$ E_0 = \min_{\gamma} E[\gamma]. \quad (8) $$

More precisely, as a corollary of the Hohenberg–Kohn theorem (later extended by Gilbert [8] to non-local potential), there exists a one-to-one mapping between the ground-state wavefunction and the 1-RDM, thus allowing to determine the ground-state energy variationally. The conditions for $\gamma$ to derive from a $N_e$-electron wavefunction – so-called N-representability conditions – are highly non-trivial [12–14] and cannot be used in practice to define the domain in which the minimization of Eq. (8) takes place. However, the conditions for $\gamma$ to be ensemble N-representable are known and straightforward [9] and are given as follows,

$$ \sum_{k=1}^{N} \eta_k = N_e, \quad 0 \le \eta_k \le 1, \quad (9) $$

where $\eta_k$ is the occupation number of the k-th natural orbital. Hence, it is common to use the natural orbital representation of the density matrix, i.e. $\gamma = (\{\eta_k\}, \{|\phi_k\}\})$. As originally proposed by Gilbert [8], these conditions
are fulfilled during the minimization by introducing a Lagrangian function in the ensemble N-representable domain,

$$ \mathcal{E} = E[\{\eta_k\}, \{|\phi_k\}\}] - \mu \left( \sum_k \cos^2(\theta_k) - N_e \right) \\ - \sum_{jk} \lambda_{jk} \left( \langle \phi_j | \phi_k \rangle - \delta_{jk} \right), \quad (10) $$

where $\mu$ and $\{\lambda_{jk}\}$ are Lagrange multipliers, the first term ($\eta_k = \cos^2(\theta_k)$) fulfills Eq. (9) and the second term the orthonormality of the natural orbitals $\{|\phi_k\}\}$. Minimizing Eq. (10) with respect to the whole set of occupation numbers and NOs leads to two non-linear differential equations,

$$ \frac{\partial \mathcal{E}}{\partial \theta_k} = \left( \mu - \frac{\partial E}{\partial \eta_k} \right) \sin(2\theta_k) = 0, \quad (11) $$

and

$$ \frac{\partial \mathcal{E}}{\partial \langle \phi_k |} = \frac{\partial E}{\partial \langle \phi_k |} - \sum_j \lambda_{jk} |\phi_j\rangle = 0, \quad (12) $$

that are solved successively but with slow convergence [60]. This situation contrasts with the variational equations in KS-DFT and SOFT that are solved efficiently through the definition of a single-particle Schrödinger equation [2, 54].

### III. LINK BETWEEN 1-RDMFT AND SOFT

Rather than decomposing the 1-RDM into natural orbitals and occupation numbers, we decompose it in terms of its diagonal elements, i.e. orbital (or site) occupations, $\mathbf{n} = \{\gamma_{ii}\}_i$ and the $N(N-1)$ off-diagonal elements denoted by $\gamma^*$ such that $\gamma = \mathbf{n} \cdot \mathbf{1} + \gamma^*$, where $\mathbf{1}$ refers to the identity matrix. The energy functional in Eq. (5) therefore reads

$$ E[\gamma] = T[\mathbf{n}, \gamma^*] + W[\mathbf{n}, \gamma^*]. \quad (13) $$

Following the KS decomposition, Eq. (13) can be alternatively written as

$$ E[\gamma] = T_s[\mathbf{n}] + F[\mathbf{n}, \gamma^*], \quad (14) $$

where $T_s[\mathbf{n}]$ represents the one-body (kinetic and potential) energy contributions of the non-interacting system with orbital occupations $\mathbf{n}$ and $F[\mathbf{n}, \gamma^*]$ is a functional containing all non trivial terms defined by

$$ F[\mathbf{n}, \gamma^*] = T[\mathbf{n}, \gamma^*] - T_s[\mathbf{n}] + W[\mathbf{n}, \gamma^*]. \quad (15) $$

In analogy with the natural orbital representation, the minimization of the energy leads to two saddle-point equations. The first saddle-point equation reads

$$ \delta \mathbf{n} \left\{ \frac{\delta T_s[\mathbf{n}]}{\delta \mathbf{n}} + \frac{\delta F[\mathbf{n}, \gamma^*]}{\delta \mathbf{n}} \middle|_{\gamma^*} \right\} = 0, \quad (16) $$

and leads to a KS-like equation:

$$ \hat{h}_{\text{KS}}[\mathbf{n}] |\Phi^{\text{KS}}\rangle = (\hat{T} + \hat{V}_{\text{Hxc}}[\mathbf{n}]) |\Phi^{\text{KS}}\rangle = \mathcal{E}^{\text{KS}} |\Phi^{\text{KS}}\rangle \quad (17) $$

with $V_{\text{Hxc}}[\mathbf{n}] := V_{\text{Hxc}}[\mathbf{n}]|_{\gamma^*} = \delta F[\mathbf{n}, \gamma^*]/\delta \mathbf{n}|_{\gamma^*}$ being the Hxc potential defined at the saddle point of $\gamma^*$ and $\hat{h}_{\text{KS}}[\mathbf{n}]$ is the KS single-particle Hamiltonian. Eq. (17) is determined concomitantly with the other saddle point equation:

$$ \delta \gamma^* \left\{ \frac{\delta F[\mathbf{n}, \gamma^*]}{\delta \gamma^*} \middle|_{\mathbf{n}=\mathbf{n}^{\text{KS}}} \right\} = 0, \quad (18) $$

where $\mathbf{n}^{\text{KS}}$ are the orbital occupations resulting from Eq. (17). In practice, solving Eq. (17) appears numerically easy when $V_{\text{Hxc}}[\mathbf{n}]$ is known, however solving Eq. (18) remains difficult, in particular due to the non-trivial representability conditions that $\gamma^*$ has to fulfill.

Interestingly, one can also use a Levy-Lieb formalism to derive similar equations,

$$ \begin{aligned} E_0 &= \min_{\mathbf{n}} \left\{ \min_{\gamma \to \mathbf{n}} \{ E[\mathbf{n}, \gamma^*] \} \right\} \\ &= \min_{\mathbf{n}} \left\{ \min_{\gamma \to \mathbf{n}} \{ T_s[\mathbf{n}] + F[\mathbf{n}, \gamma^*] \} \right\} \\ &= \min_{\mathbf{n}} \left\{ T_s[\mathbf{n}] + \min_{\gamma \to \mathbf{n}} \{ F[\mathbf{n}, \gamma^*] \} \right\}, \end{aligned} \quad (19) $$

where $\gamma \to \mathbf{n}$ corresponds to density matrices with orbital occupations $\mathbf{n}$. In that framework, the Hxc potential is defined as

$$ V_{\text{Hxc}}[\mathbf{n}] = \frac{\delta}{\delta \mathbf{n}} \min_{\gamma \to \mathbf{n}} \{F[\mathbf{n}, \gamma^*]\}, \quad (20) $$

thus bridging 1-RDMFT and SOFT, where the Hxc potential functional of the orbital occupations can be obtained after solving Eq. (18) for a fixed orbital occupation $\mathbf{n}$.

### IV. DENSITY-MATRIX INTERPOLATION VARIATIONAL ANSATZ

#### A. Topological properties of the ensemble N-representable domain of the 1-RDM

Let us consider the space $M_N(\mathbb{R})$ of $N \times N$ symmetric and real matrices, $N$ being the number of spin-orbitals in the system. In this section we denote by $\Omega_\gamma$ the domain of ensemble $N$-representable density matrices that is included in $M_N(\mathbb{R})$, $\Omega_\gamma \subset M_N(\mathbb{R})$. A significant advantage of working with natural orbitals is that the ensemble $N$-representable conditions are straightforward, i.e. a matrix $\gamma \in M_N(\mathbb{R})$ is ensemble $N$-representable if all the occupation numbers $\eta_k$ fulfill $0 \le \eta_k \le 1$. Similarly, it means that $\gamma$ is ensemble $N$-representable if both $\gamma$ and $1-\gamma$ are positive semi-definite, i.e. both showing positive or null eigenvalues $\eta_k \ge 0$ and $1-\eta_k \ge 0$ for all $1 \le k \le N$.
Hence, the domain of ensemble N-representable density matrices reads

$$ \Omega_{\gamma} = \{\gamma \in \mathcal{M}_N(\mathbb{R}) \mid \eta_k \ge 0 \ \& \ 1 - \eta_k \ge 0, \forall k \}. \quad (21) $$

Interesting topological properties for the ensemble N-representable 1-RDM domain $\Omega_\gamma$ follow from these definitions. The first remarkable property of $\Omega_\gamma$ is the convexity, see proofs in Appendix A. It follows that given a set of ensemble N-representable matrices $\{\gamma^{(i)} \in \Omega_\gamma\}$, all matrices constructed as

$$ \gamma(\mathbf{z}) = \sum_i z_i \gamma^{(i)} \quad (22) $$

also belong to $\Omega_\gamma$ if all $z_i \ge 0$ and $\sum_i z_i = 1$.

Let us now define $\tilde{\Omega}_\gamma$ the ensemble of matrices $\gamma$ such that both $\gamma$ and $1-\gamma$ are positive definite (in contrast to Eq. (21) where they are positive semi-definite), i.e.

$$ \tilde{\Omega}_{\gamma} = \{\gamma \in \mathcal{M}_N(\mathbb{R}) \mid \eta_k > 0 \ \& \ 1 - \eta_k > 0, \forall k \}. \quad (23) $$

The convexity of $\tilde{\Omega}_\gamma$ follows the same proof as in Appendix A for $\Omega_\gamma$. Besides, we introduce

$$ \delta\Omega_{\gamma} = \Omega_{\gamma} \setminus \tilde{\Omega}_{\gamma}, \quad (24) $$

such that $\delta\Omega_\gamma$ contains the density matrices $\gamma$ that exhibit at least one integer eigenvalue $\eta_k = \{0, 1\}$. Contrary to $\Omega_\gamma$ and $\tilde{\Omega}_\gamma$, $\delta\Omega_\gamma$ is not convex. We show in Appendix B that $\tilde{\Omega}_\gamma$ corresponds to the topological interior of $\Omega_\gamma$, $\delta\Omega_\gamma$ is the topological boundary of $\Omega_\gamma$ and $\Omega_\gamma = \tilde{\Omega}_\gamma \cup \delta\Omega_\gamma$ is a topological closed set.

Finally, let us define $\delta\Omega_\gamma^* \subset \delta\Omega_\gamma$ the set of idempotent matrices, i.e. $\gamma = \gamma^2$ if $\gamma \in \delta\Omega_\gamma^*$ and

$$ \delta\Omega_{\gamma}^{*} = \{\gamma \in \Omega_{\gamma} \mid \eta_{k} = \{0, 1\}, \forall k\}, \quad (25) $$

which is also non convex. $\delta\Omega_\gamma^*$ corresponds to the set of extreme points of $\Omega_\gamma$, i.e. given $\gamma_* \in \delta\Omega_\gamma^*$ there is no $\gamma^A$ and $\gamma^B$ in $\Omega_\gamma$ such that $\gamma_* = z\gamma^A + (1-z)\gamma^B$ with $0 < z < 1$. Following the Krein-Milman theorem, there exists a decomposition in terms of idempotent matrices, i.e. for any $\gamma \in \Omega_\gamma$,

$$ \gamma(\mathbf{z}) = \sum_i z_i \gamma_*^{(i)}, \quad (26) $$

with $\mathbf{z} = \{z_i\}$, $\sum_i z_i = 1$ and $\gamma_*^{(i)} \in \delta\Omega_\gamma^*$. Given Eqs. (8) and (26), it follows the variational ansatz:

$$ E_0 = \min_{\mathbf{z}} E[\gamma(\mathbf{z})]. \quad (27) $$

Hence, the energy is minimized with respect to the weights $\mathbf{z}$ defining the interpolation of $\gamma(\mathbf{z})$ in terms of idempotent 1-RDM. We show in Appendix C that for each $\gamma \in \tilde{\Omega}_\gamma$, the sum in Eq. (26) runs up to $2^N$ idempotent matrices. It follows that the number of variational

parameters rapidly becomes prohibitive when the size of the system increases. To reduce the number of parameters, one could define the trial set as every Slater determinants corresponding to the singly- and doubly-excited determinants with respect to the Hartree-Fock state. The efficiency and accuracy of such ansatz strongly depends on the size and relevance of the trial set of $\{\gamma_*^{(i)}\}$.

Rather than interpolating $\gamma$ in terms of idempotent 1-RDM, we can use the property that for each $\gamma \in \Omega_\gamma$ there exists a decomposition in terms of boundary matrices belonging to $\delta\Omega_\gamma$ [see Eq. (24)]. Then, Eq. (26) turns to

$$ \gamma(\mathbf{z}) = \sum_i z_i \gamma_B^{(i)}. \quad (28) $$

As $\gamma_B^{(i)}$ can be seen as a linear combination of $\gamma_*^{(i)}$, one can expect the amount of variational parameters $\{z_i\}$ to be drastically reduced when too much $\{\gamma_*^{(i)}\}$ would be needed to reproduce $\gamma$. Indeed, any non-idempotent $\gamma \in \tilde{\Omega}_\gamma$ can in principle be obtained by interpolating only two boundary matrices. Hence, while the set of idempotent matrices $\{\gamma_*^{(i)}\}$ in Eq. (26) is very large but straightforward to construct, only a small relevant set of boundary matrices $\{\gamma_B^{(i)}\}$ would suffice but is non-trivial to find. An efficient algorithm would consist in constructing this non-trivial set of relevant $\{\gamma_B^{(i)}\}$ on which the minimization in Eq. (27) is performed. Eqs. (27) and (28) constitute the density-matrix interpolation variational ansatz (DIVA). Note that by construction, all density matrices of the form given by Eq. (28) are ensemble N-representable.

Interestingly, the topological properties presented in $\mathcal{M}_N(\mathbb{R})$ also hold if we restrain the considered metric space to $\mathcal{M}_N^{N_e}(\mathbb{R})$ and $\mathcal{M}_N^{\mathbf{n}}(\mathbb{R})$, corresponding to the spaces of $N \times N$ symmetric and real matrices with a fixed trace (fixed number of electron $N_e$) and fixed diagonal part $\mathbf{n}$, respectively. Considering these restrictions should automatically reduce the number of trial 1-RDMs in Eq. (28).

## B. Iterative gradient-based implementations of DIVA

Let us now turn to the numerical implementations of DIVA. Starting from a guess $\gamma^{(0)}$ that already belongs to $\Omega_\gamma$, we want to construct iteratively a minimal trial set of boundary density matrices to minimize the 1-RDMFT energy [Eq. (27)] according to Eq. (28),

$$ \gamma(\mathbf{z}) = z_0 \gamma^{(0)} + \sum_i z_i \gamma_B^{(i)}. \quad (29) $$

In the following, we highlight two ways of minimizing the energy functional within DIVA. The first one consists in minimizing the energy functional with respect to only one variational parameter $z_i$ at each iteration $i$ and is
<figure><img src="image_2.png" /><figcaption>FIG. 1: Illustration of the mono-parameter DIVA. The set $\Omega_{\gamma}$ is illustrated by the union of $\delta\Omega_{\gamma}$ (blue circle) and the topological interior. Red arrows represent the direction of the energy gradient, used to find the $\gamma_B^{(i)}$ belonging to the border of the representability domain $(\delta\Omega_{\gamma})$. Each $\gamma^{(i)}$ is obtained as an interpolation between $\gamma^{(i-1)}$ and $\gamma_B^{(i-1)}$. The star represents the ground-state 1-RDM $\gamma_0$ (at the level of the functional approximation) that is recovered after convergence.</figcaption></figure>

called the mono-parameter DIVA. For the second one, the minimization is done over the whole set of parameters $\{z_j\}_{j\le i}$ at iteration $i$, and is called the multi-parameter DIVA. These approaches are pictured in Figs. 1 and 2 and detailed below.

### 1. Mono-parameter DIVA

The mono-parameter DIVA is implemented according to the following steps.

1. Initialization (iteration $i = 0$): set $\gamma^{(0)}$ to be the (idempotent) density matrix coming from an Hartree-Fock or KS-DFT calculation, thus ensuring that $\gamma^{(0)}$ belongs to $\Omega_{\gamma}$, and even to $\delta\Omega_{\gamma}^* \subset \delta\Omega_{\gamma} \subset \Omega_{\gamma}$. To avoid divergence of the energy gradient with respect to $\gamma^{(0)}$ having occupation numbers equal to 0, which happens for many NOFs, the first gradient is set to $\nabla_{ij}[\gamma] = (1 - \delta_{ij})\gamma_{ij}^{(0)}$.

2. Search the boundary density matrix $\gamma_B^{(i)} \in \delta\Omega_{\gamma}$ following the opposite direction of the energy gradient. This is done in two steps:

(a) search for $\gamma_{out} = \gamma^{(i)} - \theta\nabla[\gamma^{(i)}]$ by incrementing $\theta$ until $\gamma_{out} \notin \Omega_{\gamma}$

(b) Once $\gamma_{out}$ is found, vary $\theta$ using a bracketing method until $\gamma_{out} \in \delta\Omega_{\gamma}$ and set $\gamma_B^{(i)} = \gamma_{out}$

<figure><img src="image_3.png" /><figcaption>FIG. 2: Illustration of the multi-parameter DIVA. The set $\Omega_{\gamma}$ is illustrated by the union of $\delta\Omega_{\gamma}$ (blue circle) and the topological interior. Red arrows represent the direction of the energy gradient, used to find the $\gamma_B^{(i)}$ belonging to the border of the representability domain $(\delta\Omega_{\gamma})$. Each $\gamma^{(i)>2}$ is a linear combination of all $\gamma_B^{(j)}$ with $j < i$. The star represents the ground-state 1-RDM $\gamma_0$ (at the level of the functional approximation) that is recovered after convergence.</figcaption></figure>

3. start a new iteration, $i = i + 1$

4. Interpolate between $\gamma^{(i-1)}$ and $\gamma_B^{(i-1)}$ to find $\gamma^{(i)} = (1 - z_i)\gamma^{(i-1)} + z_i\gamma_B^{(i-1)}$ such that

$$ \frac{\partial E[\gamma^{(i)}]}{\partial z_i} = 0. \qquad (30) $$

5. Compute the energy gradient $\nabla[\gamma^{(i)}]$. Note that the number of electrons is conserved by setting the diagonal part $\nabla_{ii}[\gamma^{(i)}] = \nabla_{ii}[\gamma^{(i)}] - \mu^{(i)}$ with $\mu^{(i)} = (1/N) \sum_i \nabla_{ii}$. If the orbital occupations should be kept unchanged, set $\nabla_{ii}[\gamma^{(i)}] = 0, \forall i$.

6. If $\delta E = |E^{(i-1)} - E^{(i)}| < \nu$ and $\delta \gamma = |\gamma^{(i-1)} - \gamma^{(i)}| < \varepsilon$ where $\nu$ and $\varepsilon$ are the convergence criteria, set $\gamma^{\text{opt}} = \gamma^{(i)}$, which should be equal to the ground-state 1-RDM $\gamma_0$ corresponding to the level of approximation of the chosen 1-RDM functional (the star in Fig. 1). Else, go to step 2.

At convergence, the density matrix fulfills Eq. (28), $\gamma(\mathbf{z}) = \sum_i Z_i \gamma_B^{(i)}$, with $Z_i = z_{i+1} \prod_{i+2}^M (1 - z_i)$ and $M$ the number of iterations needed to achieve convergence. Following Ref. [53], step 2 of the algorithm could be made more efficient by projecting $\gamma^{(i)}$ to the boundary domain $\delta\Omega_{\gamma}$. However, it would lead to a $\gamma_B^{(i)}$ that deviates from
the direction pointed by the gradient. Also the opposite energy gradient direction can be replaced by the conjugate gradient. This is what we implemented in this work, as we observed a much faster convergence using it.

## 2. Multi-parameter DIVA

Apart from step 4, the multi-parameter approach remains the same as the mono-parameter one. Rather than doing the interpolation between $\gamma^{(i-1)}$ and $\gamma_B^{(i-1)}$ only, the interpolation is done over the whole set of previously determined $\gamma_B^{(j)}$ with $j < i$:

$$ \gamma^{(i)} = z^{(0)}\gamma^{(0)} + \sum_{j=1}^{i-1} z_j \gamma_B^{(j)} \quad (31) $$

The minimization is done over the whole set of parameters $\mathbf{z} = (z_0, \dots, z_{(i-1)})$,

$$ \frac{\partial E[\gamma^{(i)}]}{\partial \mathbf{z}} = 0. \quad (32) $$

In contrast to the previous implementation, we expect to achieve a much faster convergence at the expense of more expensive optimizations of the $\mathbf{z}$ parameters.

## 3. SOFT-DIVA

According to Sec. III, by splitting the variational minimization of the energy into successive minimizations over the orbital occupations and the off-diagonal elements of the 1-RDM, one can extract the Hxc potential, functional of the orbital occupations used in SOFT [see Eq. (20)]. Everytime the second saddle point equation in Eq. (18) is fulfilled, the Hxc potential of Eq. (17) is computed. In practice Eq. (18) is solved via the DIVA algorithm while setting the diagonal part of the energy gradient to zero, such that the orbital occupations are fixed (see step 5 in Sec. IV B 1). At convergence of the DIVA algorithm, the Hxc potential can be extracted by aligning the chemical potentials of the interacting and non-interacting (KS) systems as follows

$$ \begin{aligned} V_{\text{Hxc}}[\mathbf{n}] &= \left. \frac{\delta F[\mathbf{n}, \gamma^*]}{\delta \mathbf{n}} \right|_{\gamma^*} \\ &= \left. \frac{\delta E[\mathbf{n}, \gamma^*]}{\delta \mathbf{n}} \right|_{\gamma^*} - \frac{\delta T_s[\mathbf{n}]}{\delta \mathbf{n}}, \\ &= \left. \frac{\delta E[\mathbf{n}, \gamma^*]}{\delta \mathbf{n}} \right|_{\gamma^*} - \varepsilon_H[\mathbf{n}], \end{aligned} \quad (33) $$

where $\varepsilon_H[\mathbf{n}]$ is the HOMO orbital energy of the non-interacting system with fixed orbital occupations $\mathbf{n}$. Hence, the Hxc potential is updated every time the minimizing off-diagonal elements of the 1-RDM are found.

<figure><img src="image_4.png" /><figcaption>FIG. 3: Ground-state energy of the 1D Hubbard model as a function of the electron filling $n$ at different values of the Coulomb integral $U/t$. Results are given for the exact Bethe–Ansatz (gray dots), and using DIVA with the Töws–Pastor functional (full orange lines) and the Müller functional (dashed blue lines).</figcaption></figure>

This Hxc potential is then plugged into the single-particle Hamiltonian (or Fock matrix), thus generating a new idempotent 1-RDM with a new set of orbital occupations. The off-diagonal elements will then be optimized again using DIVA, and so on. Upon convergence, one obtains a numerical Hxc potential, functional of the orbital occupations, with a level of theory corresponding to the 1-RDM functional approximation used for the calculation of the energy.

# V. RESULTS AND DISCUSSIONS

We propose to study the performance of DIVA on the one-dimensional (1D) Hubbard model and on a simple molecular system. For the 1D Hubbard model, in the spirit of Mitxelenal et al. [42], we compare a NOF (the Müller functional) with a functional derived in lattice representation such as proposed by Töws and Pastor to approximate the correlation energy [61–63]. For molecular systems, only NOFs are available and we also employ the Müller functional.

## A. Hubbard model

The homogeneous Hubbard model can be defined by the following Hamiltonian,

$$ \hat{H} = -t \sum_{\langle ij \rangle, \sigma} \hat{c}_{i\sigma}^{\dagger} \hat{c}_{j\sigma} + \frac{U}{2} \sum_{i} \sum_{\sigma\sigma'} \hat{c}_{i\sigma}^{\dagger} \hat{c}_{i\sigma'}^{\dagger} \hat{c}_{i\sigma'} \hat{c}_{i\sigma} \quad (34) $$

where the notation $\langle \cdot \rangle$ refers to first neighbor lattice sites. In contrast to Eq. (1), the indices correspond to spatial orbitals (or sites) and $\sigma = \{\uparrow, \downarrow\}$ denotes the spin.
Focusing on the homogeneous model allows us to bypass Eq. (17) and to solely focus on the solution of Eq. (18) within the DIVA framework. In practice, chains with 202 sites have been considered with periodic boundary conditions, and the bisection algorithm has been employed for bracketing and line minimization in steps 2b and 4, respectively (see Sec. IV B 1). We have considered the Müller functional, defined as a function of the natural orbital occupation number (Blöch function for periodic systems), and the Töws-Pastor functional that is defined in the lattice site representation. Note that the Müller functional has been investigated in the context of finite Hubbard chains [42, 64], and the Töws-Pastor functional has been applied in the context of the Anderson model [61–63]. In terms of performance, we found that fast convergence is achieved with DIVA (5 iteration max to reach an energy difference of $\delta E/t < 10^{-7}$) for both functionals. However, computational time is smaller for the Töws-Pastor functional compared to the Müller one. Indeed, for the Müller functional depending only on natural orbital occupation numbers, a standard minimization algorithm in the natural orbital representation appears well-suited [33], while a minimization scheme directly in the lattice representation such as DIVA appears more adapted for the Töws-Pastor functional. For the latter, the computation of the energy and energy gradients (needed several time at each minimization step) do not require to switch between lattice and natural orbital representations.

In Fig. 3, the ground-state energy is shown as a function of the electron filling $n$. It shows that the Töws-Pastor functional is in fair agreement with the exact Bethe-Ansatz (BA) solution [65, 66] regardless of the electron filling and the strength of the Coulomb repulsion $U/t$. On the contrary, the Müller functional leads to inaccurate energies away from the half-band filling and especially for large $U/t$. Note that this wrong behaviour goes together with an unphysical negative value of $\partial E/\partial U$ corresponding to the number of double occupations (not shown). These unphysical results for the Müller functional arises from equivalent spin on-site exchange contributions to the energy, associated with the operator $\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}\hat{c}_{i\sigma}$ that is supposed to be nil.

In Fig. 4, the occupation numbers of the natural orbitals are represented as a function of the plane-wave number $k$, at different electron fillings $n$, different strengths of the Coulomb repulsion $U/t$, and for both the Töws-Pastor and Müller functionals. At half-band filling ($n=1$), the Töws-Pastor and the Müller functionals lead to very similar occupation numbers with a $k$-dependency that is close to a step function, which is due to the local nature of the Coulomb interaction. The step-function-like behaviour also explains the fast convergence of the algorithm as it corresponds to a quasi homogeneous renormalization of all off-diagonal terms of the density matrix in the lattice site representation, and thus to a renormalization of the kinetic energy. Turning to the three-quarter band filling, this behavior remains, however with large

<figure><img src="image_5.png" /><figcaption>FIG. 4: Occupation numbers as a function of the Blöch wave number $k$ for the half-band filled 1D Hubbard model at $U/t = 1, 4$ and $8$ (in ascending order following the black arrow). The top panel corresponds to the half-band filling ($n=1$) case while the bottom panel corresponds to a three-quarter filling case ($n=3/2$). Results are shown using the Töws-Pastor functional (full orange lines) and the Müller functional (dashed blue lines).</figcaption></figure>

discrepancies between the Töws-Pastor and the Müller functional. As already mentioned away from half-band filling, the Müller functional leads to negative double occupations in particular when $U/t$ increases, which comes along with a drastic reduction of the kinetic energy and thus drastic modification of the occupation numbers. In contrast, occupation numbers obtained with the Töws-Pastor functional are less dependent on $U/t$, which is consistent with the kinetic energy remaining finite even in the strongly interacting regime away from half-band filling.

Let us now investigate the SOFT exchange and correlation potential $V_{xc}$ obtained by Eq. (33). In Fig. 5, $V_{xc}$ obtained using the Töws-Pastor functional is compared with the exact BA one and with the Bethe-Ansatz Local Density Approximation (BALDA) [56, 57], with respect to the electron filling $n$ and for different values of $U/t$. $V_{xc}$ is not shown for the Müller functional as it leads to unphysical results away from half-band filling. A major positive result for the Töws-Pastor functional is the ability to fairly reproduce the $V_{xc}$ discontinuity at half-band filling, enabling the description of the Mott-Hubbard transition from metal to insulating regimes. This highlights the ability to design sound exchange and correlation potentials for SOFT from 1-RDMFT approximations. Away from half-filling ($0 \le n < 1$ and $1 < n \le 2$), the Töws-Pastor functional leads to a monotonous behavior of $V_{xc}$, in contrast to BALDA that shows unphysical inflection points and even positive values for $V_{xc}$, particularly in
<figure><img src="image_6.png" /><figcaption>FIG. 5: SOFT exchange and correlation potential $V_{xc}$ for the 1D Hubbard model as a function the electron filling $n$ and at different values of the Coulomb integral $U/t$. Results for the Töws-Pastor functional are compared with BALDA and with BA. The embedded subplot shows the curve associated to $U/t = 1$.</figcaption></figure>

the weakly correlated regime away from the half-filled case (see the zoomed panel in Fig. 5). This unphysical behaviour was already observed and shown in Eqs. (31) and (32) in Ref. [67]. Finally, in contrast to BALDA that has been specifically designed for the 1D Hubbard model, the Töws-Pastor functional can *a priori* be applied to any lattice dimension.

## B. Molecular systems

In contrast to the Hubbard model, both Eqs. (17) and (18) have to be solved for molecular systems. As a proof of concept, we apply our three different DIVA implementations described in Sec. IV B to the hydrogen molecule $H_2$ in the 6-31G basis, using the Müller functional. As DIVA is general and not tight to any particular orbital representation, we decided to work in the orthogonal atomic orbital (OAO) basis whose first orbital occupations are given by a DFT calculation using a local spin density approximation functional (SVWN) from the PSI4 program [68]. Numerical energy gradients with respect to $\gamma^*$ and $n$ have been computed with finite difference.

### 1. Convergence of the mono-parameter, multi-parameter, and SOFT-DIVA

In Fig. 6, we compare the convergence of the energy with respect to the number of iterations, using either the mono- or multi-parameter DIVA implementations, for an interatomic distance of $R = 1.0$ Å. A convergence criteria of $10^{-8}$ Hartree for the energy and $10^{-5}$ for the Frobenius norm of the 1-RDM have been used.

<figure><img src="image_7.png" /><figcaption>FIG. 6: 1-RDMFT energy with respect to the number of iterations, using SOFT-DIVA, mono-parameter DIVA, and multi-parameter DIVA, for $R = 1.0$ Å. The 6-31G basis and the Müller functional were used.</figcaption></figure>

As expected, the multi-parameter implementation converges in much fewer iterations than the mono-parameter one, although it does not mean that it is computationally faster than the mono-parameter DIVA. Indeed, each iteration of mono-parameter DIVA has the same cost, i.e. the cost of the optimization of a single variational parameter $z$. On the contrary, the number of parameters to optimize in multi-parameter DIVA increases at each additional iteration. Hence, there is a trade-off between the number of iterations and the cost of optimizing many parameters at once. Note that in this work, multi-parameter DIVA was always computationally much faster than mono-parameter DIVA. Besides, the former always converges in around 20 iterations while thousands of iterations are sometimes required using mono-parameter DIVA for some interatomic distances. This very slow convergence has also been observed in other gradient-based algorithms and is attributed to the fact that NOFs are not proper functionals of the 1RDM, since these functionals are not invariant with respect to the rotation of two orbitals having degenerate occupation numbers [53]. We believe that the convergence of mono-parameter DIVA could be significantly improved by using analytical gradients. Alternatively, considering a simultaneous perturbation stochastic approximation can reduce significantly the number of energy evaluations needed to compute the numerical gradient, compared to standard finite difference. This is left for future work.

Let us now consider the convergence of SOFT-DIVA, where the multi-parameter implementation has been used for the optimization of the off-diagonal elements of the 1-RDM [Eq. (18)]. As readily seen in Fig. 6, there is no numerical advantage in splitting the variational principle into two parts. Indeed, although constraining the orbital occupations fixed leads to faster convergence when solving Eq. (18), one has to do it for each new orbital occupations obtained when solving Eq. (17) with the pre-
<figure><img src="image_8.png" /></figure>

FIG. 7: 1-RDMFT energy with respect to the interatomic distance in the 6-31G basis using the Müller functional. SOFT-DIVA$_n$ corresponds to the *n*-th iteration of SOFT-DIVA, i.e. when the off-diagonal elements have been optimized for a fixed orbital occupation.

<figure><img src="image_9.png" /></figure>

FIG. 8: Occupation numbers of each natural orbital *i* denoted NO$_i$ with respect to the interatomic distance. Dots correspond to the reference results obtained with MOLGW.

viously obtained Hxc potential, which is manifested by each peak of SOFT-DIVA in Fig. 6. However, in practice it allows for the estimation of the Hxc potential, functional of the orbital occupation, at the Müller functional level of approximation.

## 2. Energies and occupation numbers

In order to check the accuracy of our implementation, we compare the energy and occupation numbers obtained by multi-parameter DIVA with the 1-RDMFT implemented in MOLGW [69, 70]. As readily seen in Figs. 7 and 8, the multi-parameter DIVA energy and occupation numbers are on top of the reference ones obtained with

<figure><img src="image_10.png" /></figure>

FIG. 9: OAO occupations with respect to the interatomic distance. For symmetry reasons, there are only two different OAO occupation values (corresponding to the 1s and 2s of the hydrogens), shown in blue and orange, respectively. Converged DIVA occupations and DFT ones are depicted with circles and squares, respectively. SOFT-DIVA$_n$ corresponds to the *n*-th iteration of SOFT-DIVA.

MOLGW. One can also see how the fractional occupation numbers of the natural orbitals in Fig. 8 are slowly getting degenerate when stretching the bond, which is a very well known result showing that 1-RDMFT can easily account for static correlation, in contrast to KS-DFT for which the 1-RDM is idempotent by construction.

Turning to SOFT-DIVA, by stopping the algorithm at each new update of the Hxc potential, one can see how much the change in orbital occupations impacts the energy. In other words, one can look at the best 1-RDM functional energy corresponding to a fixed orbital occupation (in the OAO basis), where only the off-diagonal elements of the 1-RDM have been optimized. In Fig. 7, the energy of the first iteration deviates from the reference one when stretching the hydrogen bond, meaning that the orbital occupations are not optimal in this case. This deviation is already corrected at the second iteration, corresponding to the first update of the Hxc potential that gives sufficiently accurate orbital occupations for the energy to be on top of the converged one.

Turning to the orbital occupations, the occupations of the OAO are given for the first three iterations of SOFT-DIVA in Fig. 9. Note that the initial orbital occupations of SOFT-DIVA are obtained by taking the diagonal (in the OAO basis) of the Hxc potential from the converged DFT calculation, which explains why they are not equal but close to the DFT occupations. According to Fig. 9, the OAO occupations at the second iteration of SOFT-DIVA are indeed much closer to the converged one, thus explaining why the energy in Fig. 7 matches the converged one. However, another iteration is required to get very accurate OAO occupations.
# VI. CONCLUSIONS AND PERSPECTIVES

We have introduced a new variational minimization scheme for 1-RDMFT called density-matrix interpolation variational ansatz (DIVA) that works in any orbital basis representation while satisfying the ensemble N-representability conditions. We have derived different implementations and successfully tested them on the uniform Hubbard model with the Müller and the Töws-Pastor functionals, as well as the hydrogen molecule with the Müller functional. In summary, every implementations give the correct 1-RDM and energy. The mono-parameter DIVA converges in a few iterations for the Hubbard model, but exhibits very slow convergence for the hydrogen molecule. The multi-parameter DIVA led to a much faster convergence, at the expense of more involved optimizations of the variational parameters. Finally, by allowing the estimation of the exchange and correlation potential for a given orbital occupation vector, the SOFT-DIVA implementation contributes to the development of orbital occupation functionals for chemistry. Besides, we stress that DIVA is more adapted to the computational basis representation, as it would require less swapping between the computational and the NO basis for the calculation of energy and energy gradients, thus motivating the future development of 1-RDM functionals beyond NOFs. Work is in progress in these directions.

# VII. ACKNOWLEDGMENTS

We thank Mauricio Rodriguez-Mayorga for fruitful discussions about 1-RDMFT in general. QM and MS would like to thank the ANR (Grant No. ANR-19-CE29-0002 DESCARTES project) for funding.

## Appendix A: Convexity of the ensemble N-representable 1-RDM domain.

The ensemble N-representable domain $\Omega_\gamma$ is convex if, for any ensemble N-representable 1-RDM $\gamma^A$ and $\gamma^B$, $\gamma = z\gamma^A + (1-z)\gamma^B$ with $z \in [0, 1]$ is also ensemble N-representable.
Let us consider a N dimensional vector $\mathbf{x}$. Thus we have

$$ \mathbf{x}^T \gamma \mathbf{x} = z \mathbf{x}^T \gamma^A \mathbf{x} + (1-z) \mathbf{x}^T \gamma^B \mathbf{x} \geq 0 \quad (A1) $$

since $\mathbf{x}^T \gamma^A \mathbf{x} \geq 0$ and $\mathbf{x}^T \gamma^B \mathbf{x} \geq 0$ ($\gamma^A$ and $\gamma^B$ are positive semi-definite) and $0 \leq z \leq 1$. Consequently, $\gamma$ is also positive semi-definite. Similarly,

$$ \mathbf{x}^T (1 - \gamma) \mathbf{x} = z \mathbf{x}^T (1 - \gamma^A) \mathbf{x} + (1 - z) \mathbf{x}^T (1 - \gamma^B) \mathbf{x} \geq 0 \quad (A2) $$

since $0 \leq z \leq 1$ and $\mathbf{x}^T (1 - \gamma^A) \mathbf{x} \geq 0$ and $\mathbf{x}^T (1 - \gamma^B) \mathbf{x} \geq 0$) since $1 - \gamma^A$ and $1 - \gamma^B$ are positive semi-definite. Consequently, $1 - \gamma$ is also semi positive definite. Finally,

$\gamma$ and $1 - \gamma$ being positive semi-definite, $\gamma$ is ensemble N-representable and $\Omega_\gamma$ is convex.

## Appendix B: Interior and boundary of $\Omega_\gamma$

We define the neighborhood $\mathcal{N}(\gamma)$ of a matrix $\gamma \in \mathcal{M}_N(\mathbb{R})$ as the open ball

$$ \mathcal{N}(\gamma) = \{\gamma' \in \mathcal{M}_N(\mathbb{R}) \mid \|\gamma - \gamma'\| < \epsilon, \epsilon \to 0\} \quad (B1) $$

where $\|\cdot\|$ refers to the Frobenius norm. We consider $\gamma' \in \mathcal{N}(\gamma)$ and $\delta\gamma = (\gamma' - \gamma)/\lambda$ with $\lambda = \|\gamma' - \gamma\| < \epsilon$, such that

$$ \gamma' = \gamma + \lambda \delta \gamma. \quad (B2) $$

Following a perturbation approach with $\gamma' = \gamma + \lambda\delta\gamma$ and $\lambda \to 0$, we propose to develop the eigenvalues $\eta' = \eta + \lambda\delta\eta$ and eigenvectors $\mathbf{u}' = \mathbf{u} + \lambda\delta\mathbf{u}$ as a power of infinitesimal $\lambda$ where the $\eta$ and $\mathbf{u}$ correspond to the eigenvalues diagonal matrix and eigenvector matrix of $\gamma$, respectively,

$$ \gamma \mathbf{u}^{\dagger} = \mathbf{u}^{\dagger} \eta. \quad (B3) $$

The unitary transformation $\mathbf{u}'$ that diagonalizes $\gamma'$ corresponds to

$$ \gamma' \mathbf{u}^{\dagger} = \mathbf{u}^{\dagger} \eta', \quad (B4) $$

and Eq. (B3) is recovered at the zero order development. At first order it reads

$$ \delta\gamma\mathbf{u}^{\dagger} + \gamma\delta\mathbf{u}^{\dagger} - \delta\mathbf{u}^{\dagger}\eta - \mathbf{u}^{\dagger}\delta\eta = 0. \quad (B5) $$

Multiplying by $\mathbf{u}$ on the left side leads to

$$ \begin{aligned} & \mathbf{u}\delta\gamma\mathbf{u}^{\dagger} + \mathbf{u}\gamma\delta\mathbf{u}^{\dagger} - \mathbf{u}\delta\mathbf{u}^{\dagger}\eta - \mathbf{u}\mathbf{u}^{\dagger}\delta\eta = 0 \\ & \mathbf{u}\delta\gamma\mathbf{u}^{\dagger} - \delta\eta + [\eta, \mathbf{u}\delta\mathbf{u}^{\dagger}] = 0. \end{aligned} \quad (B6) $$

On the diagonal part of Eq. (B6), the last term on the right hand side cancels, such that

$$ \delta\eta = \mathbf{u}\delta\gamma\mathbf{u}^{\dagger}. \quad (B7) $$

Let us define

$$ d(\gamma) = \min(\eta_1, 1 - \eta_N), \quad (B8) $$

where the eigenvalues are sorted in ascending order. $d(\gamma)$ appears as a pseudo-distance from $\delta\Omega_\gamma$ as $d(\gamma) \geq 0$ for $\gamma \in \Omega_\gamma$, $d(\gamma) > 0$ for $\gamma \in \hat{\Omega}_\gamma$, $d(\gamma) = 0$ for $\gamma \in \delta\Omega_\gamma$ and finally, $d(\gamma) < 0$ for non-representable matrices $\gamma \notin \Omega_\gamma$. It follows that

$$ d(\gamma') = \min(\eta_1 + \delta\eta_1, 1 - (\eta_N + \delta\eta_N)), \quad (B9) $$

and $d(\gamma')$ is contained in the interval $d(\gamma') \in ]d(\gamma) - \lambda, d(\gamma) + \lambda[$ since $|\delta\eta| \leq \lambda|\mathbf{u}\delta\gamma\mathbf{u}^{\dagger}| \leq \lambda$.

Consequently, for each $\gamma \in \hat{\Omega}_\gamma$ there is a neighborhood with $\epsilon < d(\gamma)$ such that $\mathcal{N}(\gamma) \subseteq \Omega_\gamma$, thus making $\hat{\Omega}_\gamma$ the topological interior of $\Omega_\gamma$. On the contrary, for $\gamma \in \delta\Omega_\gamma$, since $d(\gamma) = 0$ it becomes trivial to find matrices in $\mathcal{N}(\gamma)$ that do not belong to $\Omega_\gamma$. Hence, $\delta\Omega_\gamma$ is the topological boundary of $\Omega_\gamma$.
Appendix C: Density-matrix decomposition

We show that for any non-idempotent $\gamma$, there exists at least one decomposition in terms of (up to $2^N$) idempotent matrices such that $\gamma = \sum_i z_i \gamma_*^{(i)}$ with $\sum_i z_i = 1$. To that aim, we propose the following iterative algorithm. Consider $\gamma \in \Omega_\gamma \setminus \delta\Omega_\gamma^*$ that has $F(\gamma)$ fractional eigenvalues.

1. Choose one among the $F(\gamma)$ fractional eigenvalues $\eta_k$

2. Construct, in the diagonal representation of $\gamma$,

$$ \gamma'(\lambda) = \gamma + \lambda\gamma^{(k)} \quad (C1) $$

$$ \gamma_{ij}^{(k)} = 0 \forall (i, j) \neq (k, k), \text{ and } \gamma_{kk}^{(k)} = 1 \quad (C2) $$

with $\lambda$ a scalar and $\gamma'(\lambda) \in \Omega_\gamma$ for $\lambda \in [-\lambda^A, \lambda^B]$ with $\lambda^A = \eta_k$ and $\lambda^B = 1 - \eta_k$. By construction,

$$ F(\gamma'(-\lambda^A)) = F(\gamma'(\lambda^B)) = F(\gamma) - 1. $$

3. Decompose $\gamma$ in terms of matrices with higher number of integer eigenvalues,

$$ \gamma = (1 - \eta_k)\gamma'(-\lambda^A) + \eta_k\gamma'(\lambda^B) \quad (C3) $$

Repeat steps 1 and 2 for $\gamma'(-\lambda^A)$ and $\gamma'(\lambda^B)$ and insert their decomposition into Eq. (C3). Continue the algorithm until

$$ \gamma = \sum_i z_i \gamma_*^{(i)}, \quad (C4) $$

with $\gamma_*^{(i)} \in \delta\Omega_\gamma^*$ and $\sum_i z_i = 1$. Each iteration of the algorithm multiplies by two the number of matrices in the decomposition, making the sum in Eq. (C4) running ultimately up to $2^F(\gamma)$ elements.

1. P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).

2. W. Kohn and L. Sham, Phys. Rev. 140, A1133 (1965).

3. M. A. L. Marques, M. J. T. Oliveira, and T. Burnus, Comput. Phys. Commun. 183, 2272 (2012).

4. N. Mardirossian and M. Head-Gordon, Mol. Phys. 115, 2315 (2017).

5. C. J. Cramer and D. G. Truhlar, Phys. Chem. Chem. Phys. 11, 10757 (2009).

6. A. J. Cohen, P. Mori-Sánchez, and W. Yang, Science 321, 792 (2008).

7. P.-O. Löwdin, Phys. Rev. 97, 1474 (1955).

8. T. L. Gilbert, Phys. Rev. B 12, 2111 (1975).

9. A. J. Coleman, Rev. Mod. Phys. 35, 668 (1963).

10. D. T. Haar, Rep. Prog. Phys. 24, 304 (1961).

11. M. Piris and K. Pernal, Phys. Rev. A 96, 046501 (2017).

12. M. Altunbulak and A. Klyachko, Commun. Math. Phys. 282, 287 (2008).

13. C. Schilling, M. Altunbulak, S. Knecht, A. Lopes, J. D. Whitfield, M. Christandl, D. Gross, and M. Reiher, Phys. Rev. A 97, 052503 (2018).

14. F. Castillo, J.-P. Labbé, J. Liebert, A. Padrol, E. Philippe, and C. Schilling, in Annales Henri Poincaré, Vol. 24 (Springer, 2023) pp. 2241–2321.

15. I. Theophilou, N. N. Lathiotakis, M. A. L. Marques, and N. Helbig, J. Chem. Phys. 142 (2015).

16. A. J. Coleman, J. Math. Phys. 13, 214 (1972).

17. S. M. Valone, J. Chem. Phys. 73, 1344 (1980).

18. A. M. K. Müller, Phys. Lett. A 105, 446 (1984).

19. G. Csányi and T. A. Arias, Phys. Rev. B 61, 7348 (2000).

20. G. Csányi, S. Goedecker, and T. A. Arias, Phys. Rev. A 65, 032510 (2002).

21. M. A. Buijse and E. J. Baerends, Mol. Phys. 100, 401 (2002).

22. S. Sharma, J. K. Dewhurst, N. N. Lathiotakis, and E. K. U. Gross, Phys. Rev. B 78, 201103 (2008).

23. N. N. Lathiotakis, N. Helbig, A. Zacarias, and E. K. U. Gross, J. Chem. Phys. 130, 064109 (2009).

24. S. Goedecker and C. J. Umrigar, Phys. Rev. Lett. 81, 866 (1998).

25. A. Holas, Phys. Rev. A 59, 3454 (1999).

26. J. Cioslowski and K. Pernal, J. Chem. Phys. 111, 3396 (1999).

27. J. Cioslowski, C. Schilling, and R. Schilling, J. Chem. Phys. 158 (2023).

28. M. Piris and P. Otto, Int. J. Quantum Chem. 94, 317 (2003).

29. M. Piris, Recent Res. Devel. Quantum Chem. 4, 1 (2004).

30. P. Leiva and M. Piris, J. Mol. Struct. THEOCHEM 719, 63 (2005).

31. M. Piris, Int. J. Quantum Chem. 106, 1093 (2006).

32. M. Piris, X. Lopez, and J. M. Ugalde, J. Chem. Phys. 126, 214103 (2007).

33. M. Piris and J. M. Ugalde, J. Comput. Chem. 30, 2078 (2009).

34. M. Piris, J. M. Matxain, X. Lopez, and J. M. Ugalde, J. Chem. Phys. 132 (2010).

35. M. Piris, X. Lopez, F. Ruipérez, J. M. Matxain, and J. M. Ugalde, J. Chem. Phys. 134, 164102 (2011).

36. M. Piris, Int. J. Quantum Chem. 113, 620 (2013).

37. M. Piris and J. M. Ugalde, Int. J. Quantum Chem. 114, 1169 (2014).

38. X. Lopez, M. Piris, F. Ruipeérez, and J. M. Ugalde, J. Phys. Chem. A 119, 6981 (2015).

39. W. Kutzelnigg and D. Mukherjee, J. Chem. Phys. 110, 2800 (1999).

40. E. Ramos-Cordoba, X. Lopez, M. Piris, and E. Matito, J. Chem. Phys. 143, 164112 (2015).

41. J. Cioslowski, M. Piris, and E. Matito, J. Chem. Phys. 143 (2015).

42. I. Mitxelena, M. Piris, and M. Rodriguez-Mayorga, J. Phys. Condens. Matter 29, 425602 (2017).

43. M. Rodriguez-Mayorga, E. Ramos-Cordoba, M. Vianadal, M. Piris, and E. Matito, Phys. Chem. Chem. Phys. 19, 24029 (2017).

44. I. Mitxelena, M. Rodriguez-Mayorga, and M. Piris, Eur. Phys. J. B 91, 109 (2018).
45 M. Piris, Phys. Rev. Lett. **127**, 233001 (2021).

46 K. Pernal, Phys. Rev. Lett. **94**, 233002 (2005).

47 N. N. Lathiotakis, N. Helbig, A. Rubio, and N. I. Gidopoulos, Phys. Rev. A **90**, 032511 (2014).

48 I. Theophilou, N. N. Lathiotakis, N. I. Gidopoulos, A. Rubio, and N. Helbig, J. Chem. Phys. **143**, 054106 (2015).

49 S. Goedecker and L. Colombo, Phys. Rev. Lett. **73**, 122 (1994).

50 T. Baldsiefen and E. K. U. Gross, Comput. Theor. Chem. **1003**, 114 (2013).

51 M. Saubanère, M. B. Lepetit, and G. M. Pastor, Phys. Rev. B **94**, 045102 (2016).

52 J. Wang and E. J. Baerends, Phys. Rev. Lett. **128**, 013001 (2022).

53 E. Cancès and K. Pernal, J. Chem. Phys. **128**, 134108 (2008).

54 O. Gunnarsson and K. Schönhammer, Phys. Rev. Lett. **56**, 1968 (1986).

55 K. Schönhammer, O. Gunnarsson, and R. Noack, Phys. Rev. B **52**, 2504 (1995).

56 N. Lima, L. Oliveira, and K. Capelle, Europhys. Lett. **60**, 601 (2002).

57 N. A. Lima, M. F. Silva, L. N. Oliveira, and K. Capelle, Phys. Rev. Lett. **90**, 146402 (2003).

58 E. Fromager, Mol. Phys. **113**, 419 (2015).

59 J. P. Coe, Phys. Rev. B **99**, 165118 (2019).

60 N. G. Cartier and K. J. H. Giesbertz, J. Chem. Theory Comput. (2024).

61 W. Töws and G. M. Pastor, Phys. Rev. B **83**, 235101 (2011).

62 W. Töws and G. M. Pastor, Phys. Rev. B **86**, 245123 (2012).

63 W. Töws, M. Saubanère, and G. M. Pastor, Theor. Chem. Acc. **133**, 1 (2014).

64 E. Kamil, R. Schade, T. Pruschke, and P. E. Blöchl, Phys. Rev. B **93**, 085141 (2016).

65 E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. **20**, 1445 (1968).

66 H. Shiba, Phys. Rev. B **6**, 930 (1972).

67 B. Senjean, N. Nakatani, M. Tsuchiizu, and E. Fromager, Phys. Rev. B **97**, 235105 (2018).

68 D. G. A. Smith, L. A. Burns, A. C. Simonett, R. M. Parrish, M. C. Schieber, R. Galvelis, P. Kraus, H. Kruse, R. Di Remigio, A. Alenaizan, A. M. James, S. Lehtola, J. P. Misiewicz, M. Scheurer, R. A. Shaw, J. B. Schriber, Y. Xie, Z. L. Glick, D. A. Sirianni, J. S. O'Brien, J. M. Waldrop, A. Kumar, E. G. Hohenstein, B. P. Pritchard, B. R. Brooks, H. F. Schaefer, A. Y. Sokolov, K. Patkowski, A. E. DePrince, U. Bozkaya, R. A. King, F. A. Evangelista, J. M. Turney, T. D. Crawford, and C. D. Sherrill, J. Chem. Phys. **152**, 184108 (2020).

69 F. Bruneval, T. Rangel, S. M. Hamed, M. Shao, C. Yang, and J. B. Neaton, Comput. Phys. Commun. **208**, 149 (2016).

70 M. Rodriguez-Mayorga, Standalone NOFT module (1.0). Zenodo (2022), 10.5281/zenodo.7352733.