Diverging exchange force and form of the exact density matrix functional

TL;DR

This paper reveals how the exact density matrix functional in lattice models depends on the geometry of the representability polytope, introduces the concept of an 'exchange force' from divergence of the gradient, and derives the functional for the Hubbard square.

Contribution

It uncovers the geometric structure of the exact functional in RDMFT and introduces the concept of an exchange force related to the boundary of the representability polytope.

Findings

Functional depends only on natural occupation numbers within symmetry sectors.

The polytope of representable matrices is described by linear constraints.

Gradient of the functional diverges at the boundary, indicating an 'exchange force'.

Abstract

For translationally invariant one-band lattice models, we exploit the ab initio knowledge of the natural orbitals to simplify reduced density matrix functional theory (RDMFT). Striking underlying features are discovered: First, within each symmetry sector, the interaction functional $\mathcal{F}$ depends only on the natural occupation numbers $\bf{n}$. The respective sets $\mathcal{P}^1_N$ and $\mathcal{E}^1_N$ of pure and ensemble $N$-representable one-matrices coincide. Second, and most importantly, the exact functional is strongly shaped by the geometry of the polytope $\mathcal{E}^1_N \equiv \mathcal{P}^1_N $, described by linear constraints $D^{(j)}(\bf{n})\geq 0$. For smaller systems, it follows as $\mathcal{F}[\bf{n}]=\sum_{i,i'} \overline{V}_{i,i'} \sqrt{D^{(i)}(\bf{n})D^{(i')}(\bf{n})}$. This generalizes to systems of arbitrary size by replacing each $D^{(i)}$ by a linear combination of $\{D^{(j)}(\bf{n})\}$ and adding a non-analytical term involving the interaction $\hat{V}$. Third, the gradient $\mathrm{d}\mathcal{F}/\mathrm{d}\bf{n}$ is shown to diverge on the boundary $\partial\mathcal{E}^1_N$, suggesting that the fermionic exchange symmetry manifests itself within RDMFT in the form of an "exchange force". All findings hold for systems with non-fixed particle number as well and $\hat{V}$ can be any $p$-particle interaction. As an illustration, we derive the exact functional for the Hubbard square.

Figures

Equations

$$\bm{n}=\sum_{\bm{q}\in\mathcal{I}^{(Q)}}\Gamma_{\bm{q}\bm{q}}\,\bm{v}_{\bm{q}}\,\stackrel{{\scriptstyle\small{\hat{\Gamma}\equiv\mbox{$|\Psi\rangle$}\!\mbox{$\langle\Psi|$}}}}{{=}}\,\sum_{\bm{q}\in\mathcal{I}^{(Q)}}|\alpha_{\bm{q}}|^{2}\,\bm{v}_{\bm{q}}\,,$$

$${\mathcal{P}}_N^1(Q)={\mathcal{E}}_N^1(Q).$$

$$D^{(1)}(\text{bbf}n)$$

$$D^{(2)}(\text{bbf}n)$$

$$D^{(3)}(\text{bbf}n)$$

$$D^{(4)}(\text{bbf}n)$$

$${\mathcal{F}}_p[\text{bbf}n]=\underset{{\mathcal{H}}^{(Q)}\ni \Psi \mapsto \text{bbf}n}{\min }⟨\Psi |\hat{V}|\Psi ⟩.$$

$$D^{(r)}(\text{bbf}n)=|{\alpha }_r{|}^2.$$

$${\mathcal{F}}_p[\text{bbf}n]=\sum _{r,r^{\prime }}{V}_{r{r}^{\prime }}{\overline{\eta }}_r^{\ast }{\overline{\eta }}_{r^{\prime }}\sqrt{D^{(r)}(\text{bbf}n)D^{(r^{\prime })}(\text{bbf}n)}.$$

$$\left|\frac{\mathrm{d}{\mathcal{F}}_p}{\mathrm{d}\text{bbf}n}[\text{bbf}n]\right|\sim {\mathcal{G}}^{(r)}\frac{1}{\sqrt{D^{(r)}(\text{bbf}n)}},\text{mbox}as{D}^{(r)}(\text{bbf}n)\to 0.$$

$$D^{(j)}(\text{bbf}n)=\sum _{r=1}^R{D}^{(j)}({\mathbf{v}}_r)|{\alpha }_r{|}^2,$$

$${\mathcal{F}}_p[\text{bbf}n]=\sum _{r,r^{\prime }=1}^R{V}_{r{r}^{\prime }}{\overline{\eta }}_r^{\ast }{\overline{\eta }}_{r^{\prime }}\sqrt{{\tilde{D}}_r(\text{bbf}n,\hat{V})}\sqrt{{\tilde{D}}_{r^{\prime }}(\text{bbf}n,\hat{V})},$$

$$\displaystyle\tilde{D}_{r}({\bbf{n}},\hat{V})\equiv\sum^{J}_{j=1}b^{(j)}_{r}D^{(j)}({\bbf{n}})+\overline{a}_{r}\big{(}\{D^{(i)}({\bbf{n}})\},\hat{V}\big{)}\,.$$

$${\mathcal{F}}_p[\text{bbf}n]={\mathcal{F}}_p^{(j)}+{\mathcal{G}}_p^{(j)}\sqrt{D^{(j)}(\text{bbf}n)}+\mathcal{O}(D^{(j)}(\text{bbf}n)).$$

$$\mathcal{F}[\text{bbf}n]$$

$$E_0(u)/t$$

$$n_2(u)$$

$$E_0(u)/t$$

$$n_2(u)$$

$$\begin{gathered} D^{(j)}(\text{bbf}n)=\sum _{r^{\prime }=1}^R{A}_{j{r}^{\prime }}|{\alpha }_{r^{\prime }}{|}^2 \\ . \end{gathered}$$

$$\begin{gathered} \tilde{\mathcal{F}}[\{|{\alpha }_r|\}]=\sum _r{V}_{rr}|{\alpha }_r{|}^2-\sum _{r\ne r^{\prime }}|V_{r{r}^{\prime }}||{\alpha }_r||{\alpha }_{r^{\prime }}| \\ . \end{gathered}$$

$$(|{\alpha }_1{|}^2,\cdots ,|{\alpha }_R{|}^2)=\sum _{l=1}^R{a}^{(l)}{\mathbf{w}}^{(l)}$$

$$|{\alpha }_r|(\mathbf{n},\mathbf{a})$$

$$\mathcal{F}_{p}[{\bbf{n}}]=\sum^{R}_{r,r^{\prime}=1}V_{rr^{\prime}}\sqrt{\sum^{J}_{j=1}\Big{[}b^{(j)}_{r}D^{(j)}({\bbf{n}})+\overline{a}_{r}(\{D^{(j)}({\bbf{n}})\},\hat{V})\Big{]}}\,\,\sqrt{\sum^{J}_{j=1}\Big{[}b^{(j)}_{r^{\prime}}D^{(j)}({\bbf{n}})+\overline{a}_{r^{\prime}}(\{D^{(j)}({\bbf{n}})\},\hat{V})\Big{]}}\$$

$$\tilde{\mathcal{F}}[\{|{\alpha }_r|\}]$$

$$\tilde{\mathcal{F}}[\{|{\alpha }_r^{(j)}|\}]=\sum _{r,r^{\prime }\in I_j}{V}_{rr}|{\alpha }_r^{(j)}||{\alpha }_{r^{\prime }}^{(j)}|$$

$$\begin{gathered} \delta \tilde{\mathcal{F}}[\{|{\alpha }_r^{(j)}|\},\mathbf{\beta }]=2\sum _{r\in I_j,r^{\prime }\notin I_j}{V}_{r{r}^{\prime }}|{\alpha }_r^{(j)}|{\beta }_{r^{\prime }} \\ . \end{gathered}$$

$${\mathcal{F}}_p[\text{bbf}n]={\mathcal{F}}_p^{(j)}[\text{bbf}{n}^{(j)}]-2\sum _{r\in I_j,r^{\prime }\notin I_j}|V_{r{r}^{\prime }}||{\alpha }_r^{(j)}|(\text{bbf}{n}^{(j)}){\beta }_{r^{\prime }}(\text{bbf}{n}^{(j)})\sqrt{D^{(j)}(\text{bbf}n)}+\mathcal{O}(D^{(j)}(\text{bbf}n))$$

$$\tilde{\mathcal{F}}[\{|{\alpha }_r|\},\mathbf{\eta }]=\sum _{r,r^{\prime }}{V}_{r{r}^{\prime }}{\eta }_r^{\ast }{\eta }_{r^{\prime }}|{\alpha }_r||{\alpha }_{r^{\prime }}|,$$

$$D^{(r)}(\text{bbf}n)=|{\alpha }_r{|}^2.$$

$${\mathcal{F}}_p[\text{bbf}n]$$

$$\tilde{\mathcal{F}}[\{|{\alpha }_r|\}]$$

$$n_3=1-n_0,n_4=1-n_1,n_5=1-n_2.$$

Full text

Introduction.—

Reduced density matrix functional theory (RDMFT) Gilbert (1975); Cioslowski (2000); Piris (2007); Pernal and Giesbertz (2016); Schade et al. (2017) has the potential of overcoming the shortcomings and fundamental limitations of the widely used density functional theory (DFT) Hohenberg and Kohn (1964); Parr and Yang (1995); Gross and Dreizler (2013); Jones (2015). Involving the full one-particle reduced density matrix (1RDM) $\gamma$ facilitates not only an exact description of the single particle potential energy, $\mathcal{U}[\gamma]\equiv\mbox{Tr}[\hat{U}\gamma]$, but also of the kinetic energy, $\mathcal{T}[\gamma]\equiv\mbox{Tr}[\hat{T}\gamma]$. It remains to derive accurate approximations to the interaction term $\mathcal{F}[\gamma]$. Moreover, RDMFT allows explicitly for fractional occupation numbers as it is required in the description of strongly correlated systems Pernal and Giesbertz (2016). At the same time, involving the full 1RDM lies, however, also at the heart of possible disadvantages of RDMFT relative to DFT: While both methods avoid the use of exponentially complex $N$-electron wave functions, the 1RDM involves $d^{2}$ degrees of freedom compared to $d$ for the spatial density used in DFT, where $d$ is the basis set size. To be more specific, one often uses the spectral representation $\gamma\equiv\sum_{j}n_{j}|\varphi_{j}\rangle\!\langle\varphi_{j}|$ and then minimizes the total energy functional $\mathcal{E}[\gamma]=\mathcal{T}[\gamma]+\mathcal{U}[\gamma]+\mathcal{F}[\gamma]$ with respect to the natural occupation numbers (NONs) $n_{j}$ and natural orbitals $|\varphi_{j}\rangle$, separately. The dependence on the latter makes the minimization of $\mathcal{E}$ particularly difficult and one often encounters slow convergence (see, e.g., Pernal (2005)).

The general situation drastically changes in favour of RDMFT for the important class of periodic one-band lattice systems as studied in solid state physics. The 1RDM inherits the translational symmetry of the ground state Davidson (2012) and the natural orbitals are known from the very beginning. They are given for all systems by plane waves (multiplied by some spin state). Thus, various possible disadvantages of RDMFT compared to DFT disappear and RDMFT simplifies de facto to a NON-functional theory.

Based on this observation and the fact that in general the significance of symmetries in physics can hardly be overestimated, we will explore in this letter the role of the translational symmetry within RDMFT and reveal universal and far-reaching consequences. In that sense, our work complements previous studies of the homogeneous electron gas Cioslowski and Pernal (1999); Lathiotakis et al. (2007); Sharma et al. (2008); Lathiotakis et al. (2009a, b), periodic polymers Piris and Otto (2000, 2005) and of lattice systems Schindlmayr and Godby (1995); Carlsson (1997); López-Sandoval and Pastor (2000); Hennig and Carlsson (2001); López-Sandoval and Pastor (2002, 2004); Silva et al. (2005); Töws and Pastor (2011); Saubanère and Pastor (2011); Töws and Pastor (2012); Capelle and Campo Jr (2013); Töws et al. (2013); Saubanère and Pastor (2014); Carrascal et al. (2015); Di Sabatino et al. (2015); Kamil et al. (2016); Saubanère et al. (2016); Cohen and Mori-Sánchez (2016); Mitxelena et al. (2017, 2018a); Müller et al. (2018); Mitxelena et al. (2018b) in which the crucial role of symmetries was not further explored. In particular, we determine the sets $\mathcal{P}^{1}_{N}$ and $\mathcal{E}^{1}_{N}$ of pure and ensemble $N$-representable 1RDMs and show that they coincide. Then, in the form of an analytic derivation, we discover the general form of the exact functional $\mathcal{F}$ which will illustrate the fundamental role of one-body $N$-representability constraints. Finally, we show that the fermionic exchange symmetry manifests itself within RDMFT in the form of an “exchange force” which diverges on the boundary $\partial\mathcal{E}^{1}_{N}$ of the polytope $\mathcal{E}^{1}_{N}=\mathcal{P}^{1}_{N}$. All those universal features will be illustrated in two lattice cluster systems.

One-body $N$-representability constraints.—

We consider translationally invariant systems of $N$ electrons on a one-band lattice in $D$ dimensions with periodic boundary conditions and $L$ sites in each direction. Due to the translational invariance, the symmetry-adapted “orbital” part of the one-electron states are plane waves with momenta $\vec{k}=(2\pi/L)(\nu_{1},\ldots,\nu_{D})^{t}\equiv(2\pi/L)\vec{\nu}$, where $\nu_{i}=0,1,\ldots,L-1$. The spin-orbitals follow as $|\vec{\nu}m\rangle$ ($m=\pm\frac{1}{2}$) and we introduce for the following the collective quantum number $q\equiv(\vec{\nu}m)$. On the $N$-fermion level, a symmetry-adapted basis is then given by the Slater determinants \mbox{|\bm{q}\rangle}\equiv\mbox{|q_{1},\ldots,q_{N}\rangle}. The translational and spin symmetries decompose the $N$-fermion Hilbert space $\mathcal{H}$ into irreducible sectors $\mathcal{H}^{(Q)}$, $Q\equiv(\vec{K},M_{z})$, each of which is spanned by the Slater determinants \{\mbox{|\bm{q}\rangle}\}_{\bm{q}\in\mathcal{I}^{(Q)}} with total momentum $\vec{K}=\sum^{N}_{n=1}\vec{k}_{n}$ and magnetization $M_{z}=\sum^{N}_{n=1}m_{n}$. The respective set of configurations $\bm{q}$ is denoted by $\mathcal{I}^{(Q)}$.

The crucial observation is now that any two Slater determinants belonging to the same symmetry sector $Q$ differ in at least two entries $q_{n}$. As a consequence the 1RDM \mbox{\langle q|}\gamma\mbox{|q^{\prime}\rangle}=\mbox{Tr}[c_{q^{\prime}}^{\dagger}c_{q}\hat{\Gamma}] for an $N$-fermion density operator \hat{\Gamma}=\sum_{\bm{q},\bm{q^{\prime}}\in\mathcal{I}^{(Q)}}\Gamma_{\bm{q}\bm{q^{\prime}}}\mbox{|\bm{q}\rangle}\!\mbox{\langle\bm{q^{\prime}}|} (including pure states \hat{\Gamma}\equiv\mbox{|\Psi\rangle}\!\mbox{\langle\Psi|}, \mbox{|\Psi\rangle}=\sum_{\bm{q}\in\mathcal{I}^{(Q)}}\alpha_{\bm{q}}\mbox{|\bm{q}\rangle}) is diagonal. Its diagonal elements, the NONs $\bm{n}=(n_{q})$, are given by

[TABLE]

where \bm{v}_{\bm{q}}\equiv(\mbox{\langle\bm{q}|}c_{q}^{\dagger}c_{q}\mbox{|\bm{q}\rangle}) is the vector of spin-momentum occupation numbers of the Slater determinant state |\bm{q}\rangle$$\langle\bm{q}|. Its entries are one whenever $q$ is contained in $\bm{q}$ and zero otherwise. Since any $\bm{n}$ is given as the convex combination of the vectors $\{\bm{v}_{\bm{q}}\}_{\bm{q}\in\mathcal{I}^{(Q)}}$, the respective sets $\mathcal{E}_{N}^{1}(Q)$ and $\mathcal{P}_{N}^{1}(Q)$ of ensemble and pure $N$-representable 1RDMs are given as the polytope with vertices $\{\bm{v}_{\bm{q}}\}_{\bm{q}\in\mathcal{I}^{(Q)}}$ and in particular they do coincide (cf. Eq. (1)),

[TABLE]

Since not all vertices of the hypercube $[0,1]^{d}$ with particle number $N$ contribute to those sets, the $N$-representability constraints for each sector $Q\equiv(\vec{K},M_{z})$ are more restrictive than Pauli’s exclusion principle $0\leq n_{q}\leq 1$. Yet, it is important to notice that the calculation of those symmetry-adapted generalized Pauli constraints is considerably simpler than the calculation of the one-body pure $N$-representability constraints for systems without symmetries.

As an illustration, we consider three fully polarized electrons on a ring of six lattice sites with $K=0$ (for details, see supporting information 111See the Supplemental Material at url for technical details on the derivation of the exact functional, the diverging exchange force and the solution of lattice cluster systems, which includes Refs. Klyachko (2009); Schilling et al. (2017); Borland and Dennis (1972); Ruskai (2007); Fradkin (2013)). It is an elementary exercise to determine all $(\nu_{1},\nu_{2},\nu_{3})$ with $\sum_{n=1}^{3}\nu_{n}\,(\mbox{mod}\,6)=0$. One gets $(0,1,5)$, $(0,2,4),(1,2,3),(3,4,5)$ and the respective polytope (2) is then given by the convex hull of the four vertices $(1,1,0,0,0,1),(1,0,1,0,1,0),(0,1,1,1,0,0)$ and $(0,0,0,1,1,1)$. By solving linear equations this vertex representation of $\mathcal{P}_{N}^{1}=\mathcal{E}_{N}^{1}$ can be transformed into a half space representation, $\{D^{(j)}(\bm{n})\geq 0\}$, with the following four $N$-representability constraints:

[TABLE]

with the linearly dependent variables $n_{3}=1-n_{0}$, $n_{4}=1-n_{1}$ and $n_{5}=1-n_{2}$. For larger settings, the easy-to-determine vertex representation of (2) can be transformed into a half space representation by resorting to standard softwares.

Interaction functional $\mathcal{F}$ and exchange force.—

To elaborate on the structure of the exact interaction functional $\mathcal{F}$, we resort to Levy’s construction Levy (1979) (see also Ref. Lieb (1983)). For general systems (and by ignoring possible symmetries), the exact $\mathcal{F}[\gamma]$ follows as the minimization of the interaction energy over the set of all $N$-fermion pure states $|\Psi\rangle$ with 1RDM $\gamma\in\mathcal{P}_{N}^{1}$, i.e. $\mathcal{F}_{p}[\gamma]=\min_{\Psi\mapsto\gamma}\langle\Psi|\hat{V}|\Psi\rangle$. This leads to a “pure RDMFT” on $\mathcal{P}_{N}^{1}$. In practice, one tries, however, to avoid the highly intricate generalized Pauli constraints Borland and Dennis (1972); Klyachko (2006); Altunbulak and Klyachko (2008) by relaxing the minimization to $N$-fermion ensemble states $\hat{\Gamma}$ Valone (1980). This then leads to an “ensemble RDMFT” with an interaction functional $\mathcal{F}_{e}$ on the set $\mathcal{E}^{1}_{N}$ which is described by the simple Pauli exclusion principle constraints only Coleman (1963). Yet, this cannot allow one to “circumvent” the mathematically proven complexity of the ground state problem Liu et al. (2007); Schuch and Verstraete (2009) and the complexity is just shifted from the set of underlying 1RDMs to the derivation of the functional $\mathcal{F}_{e}$ and/or its minimization Schilling (2018). In that context, with regard to approximated functionals such as Cioslowski and Pernal (1999); Lathiotakis et al. (2007); Sharma et al. (2008); Lathiotakis et al. (2009a, b); Müller (1984); Goedecker and Umrigar (1998); Csányi and Arias (2000); Baerends (2001); Yasuda (2001, 2002); Buijse and Baerends (2002); Csányi et al. (2002); Herbert and Harriman (2003); Cioslowski et al. (2003); Kollmar and Heß (2003, 2004); Pernal and Cioslowski (2004); Kollmar (2004); Gritsenko et al. (2005); Piris (2005); Kollmar (2006); Rohr et al. (2008); Lathiotakis and Marques (2008); Marques and Lathiotakis (2008); Piris et al. (2011); Benavides-Riveros and Várilly (2012); Piris (2012); Pernal (2013); Piris et al. (2013); Piris and Ugalde (2014); Piris (2017); Schade et al. (2017); Piris (2018); Mitxelena et al. (2018b); Benavides-Riveros and Marques (2018), it is unclear why those based on pure state ansatzes with fixed $N$ are treated within “ensemble RDMFT”, as well. For more details the reader is referred to the reviews Pernal and Giesbertz (2016); Schade et al. (2017) and references therein.

As already stressed above, for periodic one-band lattice systems the interaction functionals simplify drastically to functionals (or more precisely to functions) of the spin-momentum occupation numbers $\bm{n}$. For each $Q\equiv(\vec{K},M_{z})$, Levy’s construction Levy (1979) is restricted to $|\Psi\rangle$ in the respective symmetry-sector (see also Refs. Gritsenko (2018); Wang and Knowles (2018))

[TABLE]

In the following, we simplify the notation by enumerating all configurations $\bm{q}\in\mathcal{I}^{(Q)}$, denote the respective Slater determinants by $|r\rangle$, $r=1,\ldots,R\equiv\dim{(\mathcal{H}^{(Q)})}$ and introduce V_{rr^{\prime}}\equiv\mbox{\langle r|}\hat{V}\mbox{|r^{\prime}\rangle}. Moreover, we will focus on $\mathcal{F}_{p}$. As it is proven in the supporting information Note (1), the equivalence $\mathcal{F}_{e}\equiv\mathcal{F}_{p}$ holds, at least whenever there exists phase factors $\eta_{r}$ such that $V_{rr^{\prime}}\equiv-\eta_{r}\eta_{r^{\prime}}|V_{rr^{\prime}}|$.

It is instructive to derive in a first step our main results for systems in which $\mathcal{P}_{N}^{1}$ takes the form of a simplex, i.e., each of its facets contains all vertices except one. Equivalently, it means that the number of independent coefficients, $\{\alpha_{\bm{q}}\}$, equals the number of independent NONs, ${\bf{n}}$. This condition is valid for several smaller systems, but also for systems of arbitrary size in case their underlying Hilbert space is restricted within (4) to a subspace involving only $\mathcal{O}(d)$ CI coefficients (yielding an approximate functional). A prime example is the one of three fully polarized electrons on six sites as already discussed above (for details see Note (1)). We thus label the one-body $N$-representability constraints $D^{(r)}(\bbf{n})\geq 0$ such that the respective facet does not contain the vertex $\bm{v}_{r}$, i.e. we have $D^{(r)}(\bm{v}_{r^{\prime}})=0$ whenever $r\neq r^{\prime}$. Moreover, we “normalize” each $D^{(r)}\geq 0$ such that $D^{(r)}(\bm{v}_{r})=1$. Using Eq. (1) and the linearity of $D^{(r)}$, we find

[TABLE]

It is exactly the simplicial structure of $\mathcal{P}_{N}^{1}$ which implies this crucial one-to-one relation between $\{D^{(r)}(\bbf{n})\}$ and $\{|\alpha_{r}|^{2}\}$. Consequently, Levy’s construction (4) with the ansatz \mbox{|\Psi\rangle}=\sum_{r}\eta_{r}|\alpha_{r}|\,\mbox{|r\rangle} is trivial to carry out up to the phase factors $\eta_{r}$ of $\alpha_{r}$. Their minimization leads to some $\overline{\eta}_{r}\equiv\overline{\eta}_{r}(\bbf{n},\hat{V})$ and eventually we obtain

[TABLE]

The result (6) for the exact interaction functional valid for any symmetry-respecting interaction $\hat{V}$ could hardly be more striking: $\mathcal{F}_{p}$ is fully determined (up to phase factors $\overline{\eta}_{r}(\bbf{n},\hat{V})$) by the geometry of the simplex $\mathcal{P}_{N}^{1}$. Moreover, the presence of an exchange force, as we shall call it, follows immediately which diverges on the boundary of $\mathcal{P}_{N}^{1}$,

[TABLE]

Remarkably, the exchange force is always repulsive in the sense that it is repelling n from the polytope boundary (see supporting information Note (1)).

Generalizing the results (6) and (7) to systems with arbitrary underlying polytope $\mathcal{P}_{N}^{1}\equiv\mathcal{E}_{N}^{1}$ is quite intricate: Relation (S11) takes the form (see supporting information Note (1))

[TABLE]

for all $j=1,\ldots,J$, where typically $D^{(j)}(\bm{v}_{r})>0$ for more than one $r$. We also introduced $J$, the number of $N$-representability constraints. As a consequence, n does not uniquely determine $\{|\alpha_{r}|\}$ anymore and instead a set of $d$ linear equations with $R>d$ variables has to be solved. The constrained search in (4) then amounts to a non-trivial minimization over the $R-d$ remaining variables. This purely technical and less informative derivation (see supporting information Note (1)) leads to the general final form

[TABLE]

The coefficients $b^{(j)}_{r}$ are solely determined by the geometry of the polytope $\mathcal{P}_{1}^{N}$ and $\overline{a}_{r}(\{D^{(j)}({\bbf{n}})\},\hat{V})$ follow from the minimization of the degrees of freedom not fixed by $\bm{n}$. This highly involved minimization, as discussed in the supporting information Note (1), leads to an implicit additional dependence of $\mathcal{F}_{p}$ on n and the interaction $\hat{V}$.

At the same time, the general form (9) offers excellent prospects for a perturbation theoretical approach by expanding $\overline{a}_{r}(\{D^{(j)}({\bbf{n}})\},\hat{V})$ (see Hubbard square below).

Whenever n approaches the facet described by $D^{(j)}\equiv 0$, it follows from Eq. (8) that $|\alpha_{r}|\to 0$ for all $r$ whose vertices ${\bbf{v}}^{(r)}$ do not belong to that facet. This fact must reflect itself in the n-dependence of $\mathcal{F}_{p}$. Indeed, one obtains for each $j$ the singular n-dependence Note (1)

[TABLE]

This result presents in a particularly striking form the crucial role of the $N$-representability constraints $D^{(j)}({\bbf{n}})\geq 0$. In particular, as an extension of (7), it confirms that the fermionic exchange symmetry manifests itself within RDMFT in the form of an exchange force diverging on the boundary of the polytope $\mathcal{P}_{N}^{1}=\mathcal{E}_{N}^{1}$.

Hubbard square.—

Now, as an illustration, we apply the general framework from above to the one-dimensional one-band Hubbard model with $N=4$ electrons, $L=4$ sites (half filling) and nearest neighbor hopping with hopping rate $t>0$. This will emphasize from a different perspective the drastic simplification of RDMFT in case all symmetries are fully exploited: The boundaries of exact functional calculation are extended from the commonly studied Hubbard dimer López-Sandoval and Pastor (2000, 2002, 2004); Carrascal et al. (2015); Di Sabatino et al. (2015); Kamil et al. (2016); Cohen and Mori-Sánchez (2016) with an underlying six-dimensional Hilbert space to the Hubbard square with a Hilbert space of dimension $70=\binom{8}{4}$.

The kinetic energy functional for the Hubbard square reads $\mathcal{T}[{\bbf{n}}]=-4t\sum_{\nu=0}^{3}\cos{(2\pi{\nu}/4)}(n_{\nu\uparrow}+n_{\nu\downarrow})$ and the Hubbard on-site interaction has strength $U\geq 0$ (Coulombic repulsion). We will present only the essential steps and refer to the supporting information Note (1), where all details of the following discussion are presented.

The ground state for $U\geq 0$ is a singlet state with total momentum $K=\frac{2\pi}{4}2=\pi$ and parity $p=-1$. Taking all these symmetries into account leads to a rather simple polytope $\mathcal{P}_{N}^{1}=\mathcal{E}_{N}^{1}\cong[0,1]$ of $N$-representable 1RDMs: It is $n_{\nu\uparrow}=n_{\nu\downarrow}\equiv n_{\nu}$, $n_{1}=n_{3}=1/2$ and $n_{0}=1-n_{2}$. Hence there is only one independent variable ($n_{2}$) (which is identified with n) constrained by Pauli’s exclusion principle $0\leq n_{2}\leq 1$, only. This is a particular incidence and in larger systems in a singlet state, the translational symmetry implies constraints which are more restrictive than Pauli’s exclusion principle.

For given n, Levy’s construction (4) cannot be fully carried out by analytical means since it involves the root of a polynomial of degree six. The exact functional $\mathcal{F}\equiv\mathcal{F}_{p}=\mathcal{F}_{e}$ Note (1) as function of $n_{2}$ is determined numerically instead and we depict it in Figure 1. Its graph demonstrates the divergence of the slope on the “facets” $n_{2}=0,1$, as predicted by (10). Also the particle-hole duality $\mathcal{F}[n_{2}]=\mathcal{F}[1-n_{2}]$ Yasuda (2001) is obvious and the convexity of $\mathcal{F}$ is consistent with the fact that “ensemble functionals” $\mathcal{F}_{e}$ are always convex Lieb (1983); Zumbach and Maschke (1985).

Using a perturbative approach for (9), the functional $\mathcal{F}$ simplifies in the asymptotic regimes of weak ($0\leq U\ll t$) and strong ($U\gg t$) coupling Note (1),

[TABLE]

Using $\mathcal{T}[{\bbf{n}}]=-8t(\frac{1}{2}-n_{2})$ and the results from Eq. (11), one obtains from the minimization of $\mathcal{E}[{\bbf{n}}]$ the ground state energy $E_{0}$ and the corresponding NON $n_{2}$ in the weak coupling regime as a function of $u=U/t$

[TABLE]

and for strong coupling

[TABLE]

The asymptotically exact results (Hubbard square.—),(Hubbard square.—) are shown in Figure 2 (left). This figure also contains the exact result and those of PNOF5 Piris et al. (2011, 2013) and PNOF7(-) Mitxelena et al. (2018b), the best approximate functionals among all used in Ref. Mitxelena et al. (2017, 2018a). Result (Hubbard square.—) fits perfectly the exact result for all $u>10$. The convergence to zero for $u\to\infty$ (a general property of the Hubbard model at half filling in any dimension Fradkin (2013)) is reproduced also by PNOF5 and PNOF7(-). In order to check the quality of the approximate functionals more, we have also plotted the relative error $\Delta E/E_{0}$ in Figure 2 (right). We observe that this error is about $60\%$ and $10\%$ for PNOF5 and PNOF7(-), respectively, and practically zero for our approximate result (Hubbard square.—) for all $u>10$.

Summary and conclusions.—

We have demonstrated how the ab initio knowledge of the natural orbitals for translationally invariant one-band lattice models significantly simplifies reduced density matrix functional theory (RDMFT). For each symmetry sector, the sets $\mathcal{P}^{1}_{N}$ and $\mathcal{E}^{1}_{N}$ of pure and ensemble $N$-representable one-matrices coincide, the interaction functionals $\mathcal{F}_{p/e}$ depend only on the natural occupation numbers n and RDMFT therefore reduces de facto to a natural occupation number “functional” theory.

Those insights have tremendous consequences. Based on Levy’s construction Levy (1979) they allowed us, to discover the form of the exact functional $\mathcal{F}_{p}[{\bm{n}}]$ (cf. (9)) which differs considerably from the approximate functionals proposed so far Pernal and Giesbertz (2016); Schade et al. (2017). Intriguingly, $\mathcal{F}_{p}[{\bm{n}}]$ is given by a bilinear form of square roots (generalizing the two-electron result Löwdin and Shull (1956)), whose radicants contain two terms. The first one is linear in the one-body $N$-representability constraints $\{D^{(j)}({\bm{n}})\}$, while the second summand depends nonlinearly on $\{D^{(j)}({\bm{n}})\}$ and on the interaction $\hat{V}$ (cf. Eq. (9)). This summand deserves particular attention: First, it arises in the constrained-search (4) from those degrees of freedom of $\Psi$ which are not determined by the one-matrix. Therefore, it represents within RDMFT irreducible correlations, a crucial concept recently established in quantum information theory Linden et al. (2002); Linden and Wootters (2002). Second, its dependence on $\hat{V}$ emphasizes that the construction of highly accurate functionals based, e.g., on tensor properties Csányi and Arias (2000); Csányi et al. (2002) or $N$-representability conditions for the 2RDM Piris (2005, 2017) would necessitate information on the interaction $\hat{V}$, as well. Third a finite series expansion of that term, \overline{a}_{r}\big{(}\{D^{(i)}({\bbf{n}})\},\hat{V}\big{)}, with respect to $\{D^{(i)}({\bbf{n}})\}$ in conjunction with a fitting scheme would allow one to establish a hierarchy of approximate functionals similar to Jacob’s ladder in DFT Perdew and Schmidt (2001).

Another potentially transformative key result of our work is the discovery of an “exchange force” emerging from the fermionic exchange symmetry: The gradient of the exact functional diverges, $\left|\mathrm{d}\mathcal{F}_{p}/\mathrm{d}\bbf{n}\right|\sim c_{i}/\sqrt{D^{(i)}(\bbf{n})}$, as n approaches a facet of the polytope $\mathcal{P}^{1}_{N}=\mathcal{E}^{1}_{N}$, defined by $D^{(i)}\equiv 0$. This repulsive divergence on the boundary of $\mathcal{E}^{1}_{N}$ also explains why fermionic occupation numbers $n_{k}$ typically cannot take the extremal values [math] or $1$. In turn, studying the equation $c_{i}(\hat{V})=0$ would allow one to systematically identify all (highly non-generic) systems (such as Cioslowski (2018)) for which occupation numbers can be pinned to [math] or $1$. It will be one of the crucial future challenges to generalize those new concepts to systems without translational symmetry, with particular focus on ensemble RDMFT (i.e., $\mathcal{F}_{e}$ on $\mathcal{E}^{1}_{N}$).

Finally, we would like to stress that all our findings hold for systems with non-fixed particle number, as well and $\hat{V}$ can be any (spin-dependent) $p$-particle interaction obeying translational symmetry.

I Derivation of the general form of $\mathcal{F}_{p}[\emph{n}]$

In a first step, we recall the one-body $N$-representability conditions. We remind the reader that we enumerated the $N$-particle configurations $\bm{q}$ by $r=1,.\ldots,R$. Then, an $N$-particle state has the form \mbox{|\Psi\rangle}=\sum^{R}_{r^{\prime}=1}\alpha_{r^{\prime}}\mbox{|r^{\prime}\rangle}, where $|r^{\prime}\rangle$ represents a Slater determinant formed from one-particle states $|q\rangle$, $q=1,\ldots,d$. $d$ is the dimension of the one-particle Hilbert space. Let us consider a single Slater determinant $|r\rangle$, i.e., $\alpha_{r^{\prime}}=1$ for $r^{\prime}=r$ and [math] otherwise. The corresponding natural occupation numbers (NONs) are denoted by the vector ${\bbf{v}^{(r)}}$. Its $q$-th component is $v^{(r)}_{q}=1$ if $|r\rangle$ contains the one-particle state $|q\rangle$, and otherwise zero. ${\bbf{v}^{(r)}}$, $r=1,\cdots,R$ are the extremal points of the set $\mathcal{P}^{1}_{N}$ of pure-state $N$-representable $\bbf{n}=(n_{q})$ and build the vertices of a polytope $P\equiv\mathcal{P}^{1}_{N}$ in the $d$-dimensional space of the NONs. Having determined the vertices one has to find the polytope’s facets. This in general is a nontrivial task. Each facet, $F_{j}$, $j=1,\cdots,J$, of $P$ is part of a $(d-1)$-dimensional hyperplane defined by $D^{(j)}({\bbf{n}})=0$ where $D^{(j)}({\bbf{n}})=\kappa^{(j)}_{0}+\sum\limits_{k=1}^{d}\kappa^{(j)}_{k}n_{k}$. The coefficients, $\kappa^{(j)}_{i}$, $i=0,1,\cdots,d$, are integers. The polytope is the intersection of the hyperplanes defined by $D^{(j)}({\bbf{n}})\geq 0$, for all $j$. Therefore, necessary and sufficient conditions for the pure-state $N$-representability are the constraints $D^{(j)}({\bbf{n}})\geq 0$ , $j=1,\cdots,J$.

The functions $D^{(j)}({\bbf{n}})$ also allow us to decompose the vertices into two sets. For given $j$ we decompose the index set $\{1,\cdots,R\}$ into a set $I_{j}=\{r_{1},\cdots,r_{j}\}$ and its complement such that ${\bbf{v}}^{(r)}\in F_{j}$ for $r\in I_{j}$ and ${\bbf{v}}^{(r)}\notin F_{j}$ otherwise. It is $D^{(j)}({\bbf{v}}^{(r)})=0$ for $r\in I_{j}$ and $D^{(j)}({\bbf{v}}^{(r)})>0$ for $r\notin I_{j}$ Klyachko (2009); Schilling et al. (2017).

What remains is the derivation of the relation between $\{D^{(j)}({\bbf{n}})\}$ and $\{|\alpha_{r}|^{2}\}$. Each function $D^{(j)}({\bbf{n}})$ determines an operator $D^{(j)}(\hat{{\bbf{n}}})$ with $\hat{n}_{q}=c^{\dagger}_{q}c_{q}$. Since $\langle\Psi|\hat{n}_{q}|\Psi\rangle=n_{q}$ we have $\langle\Psi|D^{(j)}(\hat{{\bbf{n}}})|\Psi\rangle=D^{(j)}({\bbf{n}})$. On the other hand we can substitute $|\Psi\rangle=\sum^{R}_{r^{\prime}=1}\alpha_{r^{\prime}}|r^{\prime}\rangle$ on its l.h.s.. With $({\bf{A}})_{jr}\equiv A_{jr}:=D^{(j)}({\bbf{v}}^{(r)})$ this leads to

[TABLE]

This equation establishes a relation between the NONs and $\{|\alpha_{r}|^{2}\}$ involving the functions $\{D^{(j)}({\bbf{n}})\}$ which define the domain of pure-state representability.

For the constrained minimization of $\langle\Psi|\hat{V}|\Psi\rangle=\sum_{r,r^{\prime}}V_{rr^{\prime}}\alpha^{*}_{r}\alpha_{r^{\prime}}$ for fixed n we assume for a moment that $V_{rr^{\prime}}=-\tilde{\eta}_{r}\tilde{\eta}_{r^{\prime}}|V_{rr^{\prime}}|$ for $r\neq r^{\prime}$, as well as, real coefficients $\alpha_{r}=\eta_{r}|\alpha_{r}|$) with $\tilde{\eta}_{r}=\pm 1$ and $\eta_{r}=\pm 1$. In that case the minimization with respect to the phase factors $\{\eta_{r}\}$ is accomplished by the choice $\tilde{\eta}_{r}\equiv\eta_{r}$. Then the expectation value $\langle\Psi|\hat{V}|\Psi\rangle$ takes the form

[TABLE]

To derive $\mathcal{F}[\bm{n}]$ we have to determine $\{|\alpha_{r}|\}$ as a function of n. This can be done as follows. Introducing the symmetric and semi-definite matrix ${\bf{C}}={\bf{A}}^{t}\bf{A}$ and operating with ${\bf{A}}^{t}$ on Eq. (S1) one obtains $\sum^{J}_{j=1}({\bf{A}}^{t})_{rj}D^{(j)}({\bbf{n}})=\sum^{R}_{r^{\prime}=1}C_{rr^{\prime}}|\alpha_{r^{\prime}}|^{2}$. If $d=R$, there are as many NONs as coefficients $\{|\alpha_{r}|\}$. In that case $\{|\alpha_{r}|\}$ is uniquely determined by $\{D^{(j)}({\bbf{n}})\}$, i.e. by the NONs, n. This always holds if the poytope $P$ is a simplex, which was the case for the example of three fully polarized electrons on a ring of six lattice sites. Note that there are $J$ functions $\{D^{(j)}({\bbf{n}})\}$ and only $d\leq J$ NONs. In case $d<R$ (occurs for $L$ and $N$ large enough), however, $\{|\alpha_{r}|\}$ are not uniquely determined by the NONs. In that case $\bf{C}$ has zero-eigenvalues, i.e., its rank, $d$, is smaller than $R$. Let $\{{\bf{w}}^{(l)}\}$ be the eigenvectors of $\bf{C}$ and $\{c_{l}\}$ its corresponding eigenvalues. $c_{l}>0$ for $l=1,\cdots,d$ and $c_{l}=0$ for $l=d+1,\cdots,R$. Substituting the expansion

[TABLE]

into the equation $\sum^{J}_{j=1}({\bf{A}}^{t})_{rj}D^{(j)}({\bbf{n}})=\sum^{R}_{r^{\prime}=1}C_{rr^{\prime}}|\alpha_{r^{\prime}}|^{2}$ from above and taking the orthonormality of $\{{\bf{w}}^{(l)}\}$ into account allows us to determine $a^{(l)}$ for $l=1,\ldots,d$. This yields a^{(l)}=\sum^{d}_{l=1}\big{(}c^{-1}_{l}{\bf{w}}^{(l)t}{\bf{A}}^{t}{\bf{D}}({\bbf{n}})\big{)}w_{r}^{(l)} for $l=1,\ldots,d$. Substituting these $a^{(l)}$ into the r.h.s. of Eq. (S3) we arrive at

[TABLE]

The coefficients $\{b^{(j)}_{r}\}$ follow from \sum^{d}_{l=1}\big{(}c^{-1}_{l}{\bf{w}}^{(l)t}{\bf{A}}^{t}{\bf{D}}({\bbf{n}})\big{)}w_{r}^{(l))}=\sum^{J}_{j=1}b^{(j)}_{r}D^{(j)}({\bbf{n}}) where ${\bf{D}}({\bbf{n}})=(D^{(1)}({\bbf{n}}),\cdots,D^{(J)}({\bbf{n}}))^{t}$. The absolute values $\{|\alpha_{r}|\}$ are fixed by the NONs through $\{D^{(j)}({\bbf{n}})\}$ and by the independent real variables ${\bf{a}}=(a^{(d+1)},\ldots,a^{(R)})$.

To get $\mathcal{F}_{p}$ the result Eq. (S4) has to be substituted into Eq. (S2) with a subsequent minimization with respect to ${\bf{a}}$. This is a nontrivial problem which in general can not be performed analytically. Substituting its solution ${\bf{a}}(\{D^{(j)}({\bbf{n}})\},\hat{V})$ into Eq. (S4) and this expression into Eq. (S2) yields the functional

[TABLE]

where $V_{rr^{\prime}}=-|V_{rr^{\prime}}|$ for all $r\neq r^{\prime}$ and $\overline{a}_{r}(\{D^{(j)}(\bbf{n})\},\hat{V})=\sum^{R}_{l=d+1}\,a^{(l)}(\{D^{(j)}(\bbf{n})\},\hat{V})\,w^{(l)}_{r}$. Note, the dependence on $\hat{V}$ occurs through the matrix elements $\{V_{rr^{\prime}}\}$ of $\hat{V}$.

The result (S5) simplifies for n close to a facet $F_{j}$. Remember that $D^{(j)}({\bbf{n}})\to 0$ for ${\bbf{n}}\to F_{j}$. As described above we can decompose the set $r=1,\ldots,R$ of the vertex-indices into two subsets, $I_{j}$ and its complement. Then it follows from Eq. (S1) that $|\alpha_{r}|=\sqrt{D^{(j)}({\bbf{n}})}\,{\beta}_{r}$ for all $r\notin I_{j}$ and $|\alpha_{r}|=|\alpha^{(j)}_{r}|+\mathcal{O}(D^{(j)}({\bbf{n}}))$ for all $r\in I_{j}$. The real and non-negative variables $\bm{\beta}=(\beta_{r})$ have to fulfil $\sum_{r^{\prime}\notin I_{j}}A_{jr^{\prime}}({\beta}_{r^{\prime}})^{2}=1$, which follows from Eq. (S1). $\{\alpha_{r}^{(j)}\}$ are the coefficients of the normalized $N$ particle state $|\Psi^{(j)}\rangle=\sum_{r^{\prime}\in I_{j}}\alpha_{r^{\prime}}^{(j)}|r^{\prime}\rangle$ build from Slater determinants $|r\rangle$ corresponding to the vertices of the facet $F_{j}$, only. Substituting these quantities into Eq. (S2) yields

[TABLE]

with

[TABLE]

and

[TABLE]

Again, it is $V_{rr^{\prime}}=-|V_{rr^{\prime}}|$ for all $r\neq r^{\prime}$.

$\tilde{\mathcal{F}}[\{|\alpha^{(j)}_{r}|\}]$ is like $\tilde{\mathcal{F}}[\{|\alpha_{r}|\}]$ but restricted to a subspace spanned by all basis states $|r\rangle$ with $r\in I_{j}$. Its minimization with respect to $\{|\alpha^{(j)}_{r}|\}$ has to be performed in analogy to that of $\mathcal{F}_{p}[\{|\alpha_{r}|\}]$, but now under the constraint ${\bbf{n}^{(j)}}$ fixed. $\bbf{n}^{(j)}$ is a chosen reference point in $F_{j}$ which is the limiting point of $\bm{n}\to F_{j}$. This minimization process yields $\{|\alpha^{(j)}_{r}|({\bbf{n}^{(j)}})\}$ and finally $\mathcal{F}_{p}[\{|\alpha^{(j)}_{r}|({\bbf{n}^{(j)}})\}]=\mathcal{F}^{(j)}_{p}[{\bbf{n}^{(j)}}]$. Note, the dependence of $\{|\alpha^{(j)}_{r}|\}$ on $\{V_{r,r^{\prime}}\}$, $r,r^{\prime}\in I_{j}$ is suppressed. Furthermore, $\{|\alpha^{(j)}_{r}|\}$ in the second line of Eq. (S6) has to be replaced by $\{|\alpha^{(j)}_{r}|({\bbf{n}^{(j)}})\}$. It remains the minimization of $\delta\tilde{\mathcal{F}}$ with respect to $\bm{\beta}$ which yields $\bm{\beta}({\bbf{n}^{(j)}})$ where the dependence on $\{V_{rr^{\prime}}\}$ is suppressed, as well. This completes the minimization of $\tilde{\mathcal{F}}[\{|\alpha_{r}|\}]$ for $\bf{n}$ approaching ${\bf{n}}^{(j)}$ in $F_{j}$. The functional takes the final form

[TABLE]

In case that the interaction matrix does not have the form $V_{rr^{\prime}}=-\tilde{\eta}_{r}\tilde{\eta}_{r^{\prime}}|V_{rr^{\prime}}|$ we have to minimize the functional

[TABLE]

where $\bm{\eta}=\{\eta_{r}\}$ again are the phase factors of $\{\alpha_{r}\}$. Similar as above , for fixed $(\bm{n},\bm{\eta})$ we require additional parameters ${\bf{a}}=(a^{(l)})$ in order to fix $\{|\alpha_{r}|\}$. Performing the minimization on the r.h.s. of Eq. (S10) with respect to $\bf{a}$ yields $a^{(l)}(\{D^{(j)}({\bbf{n}})\},\hat{V},\bm{\eta})$ and $\{|\alpha_{r}|\}(\{D^{(j)}({\bbf{n}})\},\hat{V},\bm{\eta})$ follows from Eq. (S4) by substituting $a^{(l)}(\{D^{(j)}({\bbf{n}})\},\hat{V},\bm{\eta})$. Then, substitution of $\{|\alpha_{r}|\}(\{D^{(j)}({\bbf{n}})\},\hat{V},\bm{\eta})$ into Eq. (S10) yields $\tilde{\mathcal{F}}[\bm{n},\bm{\eta}]$. The final step concerns the minimization with respect to the Ising-like variables $\bm{\eta}$. Let $\overline{\bm{\eta}}(\{D^{(j)}({\bbf{n}})\},\hat{V})$ denote the minimizing phase factors. The substitution of those into $\tilde{\mathcal{F}}[\bm{n},\bm{\eta}]$ yields the final result for $\mathcal{F}_{p}$ given in Eq. (9) of the main text with \overline{a}_{r}\big{(}\{D^{(j)}({\bbf{n}})\},\hat{V}\big{)}=\sum^{R}_{l=d+1}\,a^{(l)}\big{(}\{D^{(j)}({\bbf{n}})\},\hat{V},\overline{\bm{\eta}}(\{D^{(j)}({\bbf{n}})\},\hat{V})\big{)}w^{(l)}_{r}.

I.1 Derivation of the exchange force for the case of $\mathcal{P}_{N}^{1}$ being a simplex

In the case where $\mathcal{P}_{N}^{1}$ takes the form of a simplex, the derivation of the exchange force is apparently much easier (cf. Eq. (6)) than for the general case of an arbitrary polytope $\mathcal{P}_{N}^{1}=\mathcal{E}_{N}^{1}$. We prove in the following that this exchange force is repulsive in the sense that it repels n from the boundary of $\mathcal{P}_{N}^{1}$. For this, we revisit Levy’s construction where we use again the ansatz \mbox{|\Psi\rangle}=\sum_{r=1}^{R}\eta_{r}|\alpha_{r}| and assume that the interaction matrix elements $V_{rr^{\prime}}$ and therefore also the phases factors $\eta_{r}$ are real-valued, i.e., $\eta_{r}=\pm 1$. As in the main text, we label the one-body $N$-representability constraints $D^{(r)}(\bbf{n})\geq 0$ such that the respective facet does not contain the vertex $\bm{v}_{r}$, i.e. we have $D^{(r)}(\bm{v}_{r^{\prime}})=0$ whenever $r\neq r^{\prime}$. For simplicity, we “normalize” each $D^{(r)}\geq 0$ such that $D^{(r)}(\bm{v}_{r})=1$. Moreover, we recall Eq. (5), i.e.

[TABLE]

Let us now consider n very close, in a distance $\varepsilon$ to the facet described by $D^{(s)}\equiv 0$ and assume that the distances $D^{(r)}(\bbf{n})$, to all other facets are much larger, i.e., $D^{(r)}(\bbf{n})\gg D^{(s)}(\bbf{n})\equiv\varepsilon$ for all $r\neq s$. W.l.o.g. we assume $s=1$. Resorting to Levy’s construction and the general ansatz for $|\Psi\rangle$, we find

[TABLE]

We remind the reader that the one-body $N$-representability constraints read $D^{(r)}({\bbf{n}})=\kappa^{(r)}_{0}+\sum^{d}_{q=1}\kappa^{(r)}_{q}n_{q}\geq 0$. The gradient $\nabla_{{\bm{n}}}\mathcal{F}_{p}[\bbf{n}]$ contains products of $\partial\mathcal{F}_{p}/\partial D^{(r)}$ and $\nabla_{{\bm{n}}}D^{(r)}({\bbf{n}})$. The latter equals the vector ${\bm{\kappa}^{(r)}}=(\kappa^{(r)}_{1},\ldots,\kappa^{(r)}_{d})^{t}$, which is anti-parallel to the normal vector of the corresponding facet. Taking now the gradient of $\mathcal{F}_{p}[\bbf{n}]$, only the term in the middle yields a contribution which diverges in the limit $\varepsilon\rightarrow 0^{+}$ (since we assumed $D^{(r)}(\bbf{n})\gg D^{(1)}(\bbf{n})\equiv\varepsilon$ for all $r>1$). $\partial\mathcal{F}_{p}/\partial D^{(1)}$ is proportional to $1/\sqrt{D^{(1)}(\bbf{n})}$ which is positive, and its prefactor $-\,\left|\sum_{r>1}\eta_{r}V_{r1}\sqrt{D^{(r)}(\bbf{n})}\right|$ is apparently negative. Consequently, the exchange force ${\bm{f}}_{ex}(\bm{n})=-\nabla_{{\bm{n}}}\mathcal{F}_{p}[\bbf{n}]$ is parallel to ${\bm{\kappa}^{(1)}}$, i.e., it points towards the interior of the polytope. Hence, the exchange force is repulsive in the sense that it repels n from the polytope’s boundary.

II Proof of $\mathcal{F}_{p}=\mathcal{F}_{e}$

We assume $\alpha_{r}=\eta_{r}|\alpha_{r}|$ in $|\Psi\rangle=\sum^{R}_{r=1}\alpha_{r}|r\rangle$ to be real and that the interaction matrix elements are of the form $V_{rr^{\prime}}\equiv\langle r|\hat{V}|r^{\prime}\rangle=-\tilde{\eta}_{r}\tilde{\eta}_{r^{\prime}}|V_{rr^{\prime}}|$ for all $r\neq r^{\prime}$. $\eta_{r}=\pm 1$ and $\tilde{\eta}_{r}=\pm 1$. Then the minimization of the expectation value $\langle\Psi|\hat{V}|\Psi\rangle$ with respect to the phase factors $\{\eta_{r}\}$ is done for $\eta_{r}\equiv\tilde{\eta}_{r}$ leading to

[TABLE]

Choose an $N$-particle ensemble $\hat{\Gamma}=\sum^{R}_{r,r^{\prime}=1}\Gamma_{rr^{\prime}}|r\rangle\langle r^{\prime}|$ . Then it follows $\langle\hat{V}\rangle_{\hat{\Gamma}}=Tr_{N}(\hat{V}\hat{\Gamma})=\sum^{R}_{r=1}V_{rr}\Gamma_{rr}-\sum^{R}_{r\neq r^{\prime}=1}\tilde{\eta}_{r}\tilde{\eta}_{r^{\prime}}|V_{rr^{\prime}}|\Gamma_{rr^{\prime}}$. A necessary condition for $\hat{\Gamma}\geq 0$ is $|\Gamma_{rr^{\prime}}|^{2}\leq\Gamma_{rr}\Gamma_{r^{\prime}r^{\prime}}$ for all $r\neq r^{\prime}$. The choice $\Gamma_{rr^{\prime}}=\tilde{\eta}_{r}\tilde{\eta}_{r^{\prime}}\sqrt{\Gamma_{rr}\Gamma_{r^{\prime}r^{\prime}}}$ minimizes $\langle\hat{V}\rangle_{\hat{\Gamma}}$ for fixed diagonal elements $\{\Gamma_{rr}\}$ and leads to $\min_{\{\Gamma_{rr^{\prime}}\}_{r\neq r^{\prime}}}\langle\hat{V}\rangle_{\hat{\Gamma}}\equiv\tilde{\mathcal{F}}[\{\sqrt{\Gamma_{rr}}\}]$. This choice also implies $\hat{\Gamma}=|\Phi\rangle\langle\Phi|$ with $|\Phi\rangle=\sum^{R}_{r=1}\tilde{\eta}_{r}\sqrt{\Gamma_{rr}}\,|r\rangle$, i.e. the corresponding $N$-particle density operator, $\hat{\Gamma}$, is positive semi-definite. Final minimization of $\tilde{\mathcal{F}}\{|\alpha|_{r}\}$ and $\tilde{\mathcal{F}}\{\sqrt{\Gamma_{rr}}\}$ with respect to $\{|\alpha|_{r}\}$ and $\{\sqrt{\Gamma_{rr}}\}$ under constraints $\bm{n}=\{n_{q}\}$ fixed, leads to $\mathcal{F}_{p}=\mathcal{F}_{e}\equiv\mathcal{F}$.

III Derivation of $\mathcal{F}_{p}[\emph{n}]$ for $N$=3 fully polarized electrons in one dimension and $L$=6

The one-particle momenta $k=(2\pi/6)\,\nu$ from the first Brillouin zone are chosen as $\nu=0,1,2,3,4,5$. Taking only nearest neighbor hopping into account this leads to the one-particle energies $\varepsilon_{\nu}=-2t\cos(2\pi\nu/6)$. $t>0$ is the nearest neighbor hopping parameter. The ground state for noninteracting spinless electrons ( i.e., $\hat{V}\equiv 0$) is \mbox{|\nu_{1},\nu_{2},\nu_{3}\rangle}^{(0)}=|0,1,5\rangle for which the total momentum is $K=(2\pi/6))(1+5)$(mod$6$)$=0$. If $|\langle{\bm{q}}|\hat{V}|{\bm{q}^{\prime}}\rangle|$ for the $3$-particle states $|{\bm{q}}\rangle=c^{\dagger}_{q_{1}\uparrow}c^{\dagger}_{q_{2}\uparrow}c^{\dagger}_{q_{3}\uparrow}|0\rangle$ is below a critical value for all ${\bm{q}}$, ${\bm{q}^{\prime}}$ with $K=0$ the ground state of the interacting system will stay in this symmetry sector. It is easy to show that the zero-momentum space is spanned by four states \mbox{|\nu_{1},\nu_{2},\nu_{3}\rangle}=|0,1,5\rangle, $|0,2,4\rangle$, $|1,2,3\rangle$ and $|3,4,5\rangle$, denoted by $|r\rangle$, $r=1,\cdots,4$. Then a general three-particle state in this symmetry sector is represented as $|\Psi\rangle=\sum\limits_{r=1}^{4}\alpha_{r}|r\rangle$. Note that the number, $d\equiv L=6$, of one-particles states $|q\rangle$ is larger than, $R=4$, the dimension of the three-particle subspace. This implies besides the normalization $\sum^{5}_{\nu=0}n_{\nu}=3$ additional identities for the NONs ($n_{\nu}$) independent on $\{\alpha_{r}\}$. Since fully polarized electrons correspond to spinless fermions, the spin variables are suppressed.

It is straightforward to determine $\{n_{\nu}\}$ as a function of $\{\alpha_{r}\}$. From this relation and the normalization condition for $\{\alpha_{r}\}$ one obtains

[TABLE]

Accordingly, there are three-independent NONs, only. We choose ${\bbf{n}}=(n_{0},n_{1},n_{2})$.

Now one could follow the general scheme described in the main text to determine the vertices of the polytope and then the facets which yields the functions $\{D^{(j)}({\bbf{n}})\}$. Since for the present case the set of linear equations relating $\{n_{\nu}\}$ and $\{\alpha_{r}\}$ is rather simple one can solve this set directly. One obtains for $r=1,\cdots,4$

[TABLE]

with

[TABLE]

The validity of $\sum\limits_{r=1}^{4}|\alpha_{r}|^{2}=1$ is obvious. Eq. (III) ist identical to Eq. (3) of the main text.

$|\alpha_{r}|\geq 0$ and Eq. (S15) yields the generalized constraints, $D^{(j)}({\bbf{n}})\geq 0$,  $j=1,\ldots,4$ on n. These guarantee that n is pure-state $N$-representable, in the sector $K=0$. They define four planes building a three-dimensional polytope (a tetrahedra, which is a simplex). This polytope is identical to that of the so-called Borland-Dennis setting for three spinless fermions in a six-dimensional one-fermion Hilbert space without any symmetry conditions Borland and Dennis (1972); Ruskai (2007). Substituting $\{|\alpha_{r}|\}$ from Eq. (S15) into \mbox{\langle\Psi|}\hat{V}\mbox{|\Psi\rangle}=\sum_{r,r^{\prime}}V_{rr^{\prime}}\eta_{r}\eta_{r^{\prime}}\,|\alpha_{r}|\,|\alpha_{r^{\prime}}| and minimizing with respect to $\{\eta_{r}\}$ one obtains the final result which is of the form of Eq. (6) with $\{D^{(j)}(\bm{n})\}$ from Eq. (III).

IV Derivation of $\mathcal{F}[\emph{n}]$ for the Hubbard-square, $N=4$, $L=4$, $K=2(2\pi/4)$, $S=0$ and parity $p=-1$

To determine all Slater determinants $|k_{1}m_{1},k_{2}m_{2},k_{3}m_{3},k_{4}m_{4}\rangle$ with total momentum $K=\sum\limits_{n=1}^{4}k_{n}({\rm mod}\,2\pi)=2\pi/4\,\sum\limits_{n=1}^{4}\nu_{n}({\rm mod}\,4)=2\cdot 2\pi/4$ and total magnetization $M_{z}=\sum\limits_{n=1}^{4}m_{n}=0$ ($m_{n}=\pm 1/2$) is straightforward. With $\nu_{n}\in\{0,1,2,3\}$ skipping $2\pi/4$ and use of $+1/2=\uparrow,\,-1/2=\downarrow$ one obtains ten states

[TABLE]

Due to the isotropy in spin space and the reflection symmetry $P:i\to L-i+1$ implying $P:\nu\to-\nu(modL)$ all basis states can be chosen to be eigenstates of the operator of the total spin squared, $\vec{\hat{S}}^{2}$, and the parity operator $\hat{P}$ with eigenvalues $S(S+1)$ and $p=\pm 1$, respectively. The ground state for zero interactions is two-fold degenerate. The degeneracy is lifted in first order in $U$. The corresponding groundstate for $U=0^{+}$ is given by $\frac{1}{\sqrt{2}}[|0\uparrow,0\downarrow,3\uparrow,3\downarrow\rangle-|0\uparrow,0\downarrow,1\uparrow,1\downarrow\rangle]$ which is an eigenstate of $\vec{S}^{2}$ and $\hat{P}$ with eigenvalues [math] and $p=-1$, respectively.

Then we get for $S=0$ and $p=-1$ the following three basis states

[TABLE]

which will be denoted by $|r\rangle$ , $r=1,\cdots,3$. With $|\Psi\rangle=\sum\limits^{3}_{r=1}\alpha_{r}|r\rangle$ it is straightforward to express the NONs by $\{|\alpha_{r}|^{2}\}$ :

[TABLE]

where the normalization $\sum\limits^{3}_{r=1}|\alpha_{r}|^{2}=1$ was used. Due to $S=0$ it is $n_{\mu\uparrow}=n_{\mu\downarrow}$. Furthermore Eq. (IV) implies $n_{0\uparrow}+n_{2\uparrow}=1$. With $n_{\mu}\equiv n_{\mu\uparrow}=n_{\mu\downarrow}$ it follows from Eq. (IV)

[TABLE]

Accordingly there is one independent NON, only. We choose $n_{2}$ and identify n with $n_{2}$ , being restricted to $0\leq n_{2}\leq 1$. Therefore the “facets” are defined be $D^{(1)}({\bbf{n}})=1-n_{2}=0$ and $D^{(2)}({\bbf{n}})=n_{2}=0$

The matrix $(A_{jr})\equiv(D^{(j)}({\bbf{v}}^{(r)}))$ in Eq. (S1) becomes

[TABLE]

which leads to (see part I)

[TABLE]

The eigenvalues $\{c_{l}\}$and the corresponding orthonormalized eigenvectors $\{(w^{(l)}_{r})\}$ are $c_{1}=1\ ,c_{2}=3/2\ ,c_{3}=0$ and $(w^{(1)}_{r})=(1/\sqrt{2})(1,-1,0)^{t}\ ,(w^{(2)}_{r})=(1/\sqrt{3})(1,1,1)^{t}\ ,(w^{(3)}_{r})=(1/\sqrt{6})(1,1,-2)^{t}$, respectively. Substituting these expressions into Eq. (S3) yields

[TABLE]

It is straightforward to calculate the matrix elements $V_{rr^{\prime}}=\langle r|\hat{V}|r^{\prime}\rangle$ of the Hubbard interaction $\hat{V}=U\sum^{4}_{i=1}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}$. As a result one obtains

[TABLE]

Note, For $U>0$ all nondiagonal elements are negative. Therefore, the Hubbard interaction belongs to the class of pair interactions for which Eqs. (S5) and (S9) hold. Therefore it is $\mathcal{F}_{p}[\bbf{n}]=\mathcal{F}_{e}[\bbf{n}]\equiv\mathcal{F}[\bbf{n}]$

To get the functional $\mathcal{F}[\bbf{n}]$ we have to minimize $\sum^{3}_{r,r^{\prime}=1}V_{rr^{\prime}}\alpha_{r}\alpha_{r^{\prime}}$ with respect to the single independent degree of freedom, $a^{(3)}$, and have to follow the scheme described in part I . This leads to $\mathcal{F}[\bbf{n}]$ of the form of Eq. (9). Here we illustrate another route where first $(V_{rr^{\prime}})$ is diagonalized. The eigenvalues of $(V_{rr^{\prime}})$ are $U(1,2,0)$ with corresponding orthonormalized eigenvectors $\vec{v}_{1}=(1/\sqrt{2})(1,-1,0)^{t}$, $\vec{v}_{2}=1/(2\sqrt{2})(1,1,-\sqrt{6})^{t}$ and $\vec{v}_{3}=(\sqrt{3/8})(1,1,2/\sqrt{6})^{t}$. That one of its eigenvalues vanishes, will be crucial when we will discuss the strong coupling limit $U\to\infty$ below.

Using the eigenvectors $\{\vec{v}_{r}\}$ it follows $\vec{\alpha}=(\alpha_{1},\alpha_{2},\alpha_{3})^{t}=\sum^{3}_{r=1}\bar{\alpha}_{r}\vec{v}_{r}$ with $\sum^{3}_{r=1}\bar{\alpha}_{r}^{2}=1$, where $\bar{a}_{r}=(\vec{\alpha}\cdot\vec{v}_{r})$ are real. Then it is

[TABLE]

The second line of Eq. (IV) (the constraint for the independent NON, $n_{2}$) becomes

[TABLE]

from which it follows

[TABLE]

$\bar{\alpha}^{(\pm)}_{2}(n_{2};\bar{\alpha}_{1})$ put into the r.h.s. of Eq. (S33) yields the functional

[TABLE]

which has to be minimized with respect to $\bar{a}_{1}$. Note that $\bar{\alpha}^{(+)}_{2}(n_{2};\bar{\alpha}_{1})\equiv-\bar{\alpha}^{(-)}_{2}(1-n_{2};\bar{\alpha}_{1})$. Since $\tilde{\mathcal{F}}[{\bbf{n}};\bar{\alpha}_{1}]$ involves \big{(}\bar{\alpha}^{(\pm)}_{2}(n_{2};\bar{\alpha}_{1}\big{)}^{2} the minimization with respect to $\bar{\alpha}_{1}$ yields a functional $\mathcal{F}[n_{2}]$ exhibiting the particle-hole symmetry $\mathcal{F}[n_{2}]=\mathcal{F}[1-n_{2}]$ [55]. The resulting equation from that minimization is a polynomial in $\bar{\alpha}_{1}^{2}$ of degree six. Its roots can not be calculated analytically. Therefore, we use this situation to demonstrate the power of a perturbative approach leading in the weak and strong coupling limit $0\leq U\ll 1$ and $U\gg 1$, respectively, to asymtotically exact results, obtained analytically.

IV.1 weak coupling limit

The ground state for $U\to 0^{+}$ is twofold degenerate (both states in the first line of Eq. (IV)). In first order in $U$ the degeneracy is lifted and the ground state is given by the state in the first line of Eq. (IV). Consequently the coefficients $\alpha_{r}$ in $|\Psi\rangle$ must fulfil $\alpha_{1}\to 1$ and $\alpha_{r}\to 0$ for $r=2,3$ which implies $\bar{\alpha}_{1}=1/\sqrt{2}+x_{1}$and $\bar{\alpha}_{2}=1/(2\sqrt{2})+x_{2}$ with $|x_{r}|\to 0$. Taking this small-U dependence into account it follows from Eq. (IV) that $n_{2}\to 0$, i.e., $n_{2}$ is the smallness parameter for the pertubative calculation of $\mathcal{F}$ for weak coupling. This is intuitively clear, because it is the occupation number of the highest one-particle level $|\nu=2\rangle$. Eq. (S36) becomes

[TABLE]

where $h.o.t.$ stands for higher order terms. Expanding the r.h.s. of Eq. (S35) with respect to $x_{1}$ leads to

[TABLE]

where the minus-sign in Eq. (S35) has to be used. The plus-sign has to be chosen for $n_{2}\to 1$. Putting this result into Eq. (S37) and minimizing with respect to $x_{1}$ one obtains

[TABLE]

The fact that the leading order of $x^{(min)}_{1}$ is proportional to $\sqrt{n_{2}}$ justifies a postiori that we have not taken into account the higher order terms in Eq. (S37). Finally substituting $x^{(min)}_{1}$ into the r.h.s. of Eq. (S38) we get from Eq. (S37) the functional in the weak coupling limit

[TABLE]

IV.2 strong coupling limit

It is known that the ground state energy of the Hubbard model at half filling converges to zero for $U\to\infty$ Fradkin (2013). Therefore, it follows from Eq. (S33) that in the strong coupling limit it must be $\bar{a}_{1}\to 0$ and $\bar{a}_{2}\to 0$, i.e. only the eigenvector $\vec{v}_{3}$ of $(V_{r^{\prime}r})$ with eigenvalue [math] contributes on the l.h.s. of Eq. (S33) . From Eq. (S34) we obtain $n_{2}\to(1/2)^{-}$, i.e., $\delta=(\frac{1}{2}-n_{2})\geq 0$ is the smallness parameter for the perturbative construction of the functional $\mathcal{F}$ in the strong coupling limit. Since $\bar{\alpha}_{1}\to 0$, Eq. (S35) implies first that we have to choose the minus-sign, and second the square root must converge to $\delta$, in order that $\bar{\alpha}_{2}\to 0$. The latter condition becomes satisfied if

[TABLE]

Substituting this expression into Eq. (S35) (with the minus-sign) and expanding its r.h.s. with respect to $\delta$ and $y_{1}$ leads to

[TABLE]

Next, $\bar{\alpha}_{1}$ and $\bar{\alpha}_{2}$ from Eqs. (S41) and (S42) are put into Eq. (S36) which yields

[TABLE]

Its minimum is taken at

[TABLE]

In a final step $y_{1}(\delta)$ is substituted into the r.h.s. of Eq. (S43) leading to the exact functional

[TABLE]