Phase dilemma in natural orbital functional theory from the N-representability perspective

Any rigorous approach to first-order reduced density matrix (Γ) functional theory faces the phase dilemma, that is, having to deal with a large number of possible combinations of signs in terms of the electron–electron interaction energy. This problem was discovered by reducing a ground-state energy generated from an approximate N-particle wavefunction into a functional of Γ, known as the top-down method. Here, we show that the phase dilemma also appears in the bottom-up method, in which the functional E[Γ] is generated by progressive inclusion of N-representability conditions on the reconstructed two-particle reduced density matrix. It is shown that an adequate choice of signs is essential to accurately describe model systems with strong non-dynamic (static) electron correlation, specifically, the one-dimensional Hubbard model with periodic boundary conditions and hydrogen rings. For the latter, the Piris natural orbital functional 7 (PNOF7), with phases equal to −1 for the inter-pair energy terms containing the exchange-time-inversion integrals, agrees with exact diagonalization results.


Introduction
The first-order reduced density matrix (Γ ) functional theory, that is, the theory where the ground-state energy (E) is represented in terms of Γ , has emerged in recent years as a promising method to study strongly correlated materials. The seminal article of Gilbert [1] on the existence of the functional along with the works of Donnelly and Parr [2], Levy [3] and Valone [4] laid the foundations, however, the computational schemes based on these exact formulations are several times more expensive than solving directly the Schrödinger equation, so practical applications require a different approach for E [Γ ].
In 1967 [5], Rosina demonstrated that there is a one-to-one mapping from the ground-state two-particle reduced density matrix (D) to the N-particle wavefunction in the case of a Hamiltonian with at most two-body interactions. Taking advantage of the Rosina's theorem, the existence theorem of Gilbert implicitly establishes a one-to-one correspondence between the ground-state D and Γ , therefore, the functional E [Γ ] must match the exact well-known functional E [D]. It should be noted that the unknown functional in a Γ -based theory is the electron-electron potential energy (V ee ) since the rest of the Hamiltonian is actually a single-particle operator. Unfortunately, the exact reconstruction V ee [Γ ] has been an unattainable goal so far, and we have to settle for making approximations. The typical approach is to employ the exact energy expression E [D] but using solely a reconstruction functional D [Γ ]. Approximating the energy functional in that way implies that theorems obtained for the exact functional E [Γ ] are no longer valid, since an approximate functional still depends on D [6]. An undesired implication is that the functional N-representability problem arises [7][8][9], that is, we have to comply the requirement that D reconstructed in terms of Γ must satisfy the same N-representability conditions as those imposed on unreconstructed two-matrices to ensure a physical value of the approximate ground-state energy; otherwise, there will not be an N-electron system with an energy value E [D].
In general, Γ -based functionals have been proposed using reasonable heuristic or physical arguments [10], so that most of the approximate functionals currently in use are not N-representable, and that is why energy is often obtained far below true energy. It has been assumed that there is no N-representability problem of the functional, since it is believed that only N-representable conditions [11] on Γ are sufficient, but the latter is only true for the exact reconstruction of V ee [Γ ]. The ensemble N-representability constraints for acceptable Γ are easy to implement, but are insufficient to guarantee that the reconstructed D is N-representable, and thereby the approximate functional either. To date, only the functionals proposed by Piris and coworkers [12][13][14][15] relies on the reconstruction of D subject to ensemble N-representability conditions.
One issue related to the functional N-representability is that approximate N-representable functionals are not invariant with respect to a unitary transformation of the orbitals. The fact is that apart from the simple Hartree-Fock (HF) approximation, none of the known approximate functionals are explicitly given in terms of Γ , including the familiar functional which accurately describes twoelectron closed-shell systems [16,17]. It is worth noting that there are functionals [18][19][20][21] that seem to depend properly on Γ . However, these functionals violate the antisymmetric requirement for D [22], consequently none of these functionals affords an N-representable two-matrix, nor can they reproduce the simplest two-electron case. Extensive N-representability violations have been recently reported [23] for this kind of functionals.
Approximations for E [Γ ] can be obtained essentially using two methods, namely, the top-down and bottomup methods [8,24]. The top-down method consists in the reduction of an N-particle ground-state energy generated from an approximate wavefunction into a functional of Γ , whereas, in the bottom-up method E [Γ ] is generated by progressive inclusion of N-representability conditions [25] on the reconstructed D [Γ ].
The use of the top-down method with a parametrization of coefficients in a configuration interaction (CI) expansion reveals a very serious bottleneck affecting any rigorous approach to E [Γ ], namely the phase dilemma that stems from the necessity to carry out minimization over a large number of possible combinations of CI coefficient signs [26]. As expected, the phase dilemma also appears when the bottom-up method is used, i.e., we have to deal with a large number of possible combinations of signs in terms of the electron-electron interaction energy.
In the next section, we analyze how the phase dilemma arises when applying N-representabilty conditions to the reconstructed D [Γ ]. In Sections 3 and 4, we demonstrate that a suitable choice of signs is essential to describe accurately model systems with strong non-dynamic (static) electron correlation. This leads us to the formulation of the Piris natural orbital functional 7 (PNOF7) with phases equal to −1 for the inter-pair energy terms containing the exchange-time-inversion integrals, which captures the electron correlation energy close to the exact diagonalization values.

Natural orbital functional theory and N-representability
The present-day functionals are constructed in the basis where Γ is diagonal, which is the definition of a natural orbital functional (NOF) [27,28]. In this context, the natural orbitals (NOs) are the orbitals that diagonalize the one-matrix corresponding to an approximate ground-state energy, so it is more appropriate to speak of a NOF rather than a functional of Γ due to the explicit dependence on D mentioned above for approximate functionals.
Accordingly, the ground-state electronic energy is given in terms of the NOs and their occupation numbers (ONs), namely, where H ii denotes the diagonal elements of the core-Hamil tonian, kl|ij are the matrix elements of the two-particle interaction, and D[n i , n j , n k , n l ] represents the reconstructed D from the ONs. Restrictions on the ONs to the range 0 ≤ n i ≤ 1, also known as Pauli constraints, represent the necessary and sufficient conditions for ensemble N-representability of Γ [11] under the normalization condition i n i = N. On this respect, it is worth noting that we focus on the N-representability problem for statistical ensembles. Conditions named generalized Pauli constraints have been obtained [29,30] for pure-state N-representability of Γ . However, the number of these conditions increases dramatically with the number of NOs, so it becomes quite difficult to handle them in practical implementations [31]. Anyway, in order to guarantee the pure-state Nrepresentability conditions in the minimization of E [Γ ] only Pauli constraints are necessary if the functional is the appropriate one [4,32]. Indeed, if the approximate NOF is pure-state N-representable, i.e., it is obtained from the reconstruction of a pure-state N-representable two-matrix, then contraction of D will always lead to a pure-state N-representable one-matrix.
It is clear that the construction of an N-representable functional given by (1) is related to the N-representability problem of D[n i , n j , n k , n l ]. The use of ensemble Nrepresentability conditions [33] for generating a reconstruction functional was proposed in reference [12]. This particular reconstruction is based on the introduction of two auxiliary matrices ∆ and Π expressed in terms of the ONs to reconstruct the cumulant part of D [34]. In this work, we address only singlet states and adopt a restricted spin theory, so that energy (1) becomes where J pq , K pq , and L pq are the direct, exchange, and exchange-time-inversion integrals [35]. Appropriate forms of matrices ∆ and Π lead to different implementations known in the literature as PNOFi (i = 1-7) [12][13][14][15].
Remarkable is the case of PNOF5 [36,37] which turned out to be strictly pure N-representable [38,39]. The conservation of the total spin allows us to derive the diagonal elements ∆ pp = n 2 p and Π pp = n p [40]. The N-representability D and Q conditions of the two-matrix impose the following inequalities on the off-diagonal elements of ∆ [12], whereas to fulfill the G condition, the elements of the Πmatrix must satisfy the constraint [7] Here, h p denotes the hole 1 − n p . Notice that for singlets the total hole for a given spatial orbital p is 2h p . For a given approximation of ∆ qp that satisfies the inequalities of equation (3), it is evident that the modulus of Π matrix elements is determined from equation (4) assuming the equality, however, there is not any hint to determine the sign of Π qp . Consequently, a large number of possible combinations of these signs looms up for those terms containing exchange-time-inversion integrals L pq in equation (2). We need to solve this phase dilemma to propose a NOF, that is, make an adequate selection of the Π qp signs.
We now focus on the simplest case of two electrons. Fortunately, an accurate NOF is well-known for this system from the exact wavefunction [16] assuming that all amplitudes, with the exception of the first one, are negative if the first amplitude is chosen to be positive [17]. The two-electron singlet energy reads as It is worth mentioning that in some stretched twoelectron molecules, small contributions to energy may have opposite signs to those adopted in equation (5) [41,42]. A recent study [23] on the two-electron Harmonium atom reveals similar small deviations in the highcorrelation regime. Nevertheless, the convention of signs adopted in equation (5) provides very accurate results for almost all correlation regimes in two-electron systems, including those with strong non-dynamic correlation.
The requirement that for any two-electron singlet system the NOF (2) yields the accurate expression (5), together with the cumulant sum rules, and the Nrepresentability conditions (3) and (4), imply that ∆ qp = n q n p and |Π qp | = √ n q n p , respectively [7]. Furthermore, the phase factor of Π qp is +1 if q, p ∈ (1, ∞), and −1 otherwise.
For systems with N > 2, the generalization of these expressions for auxiliary matrices ∆ and Π leads to the independent-pair model (PNOF5) [36,37]: where we have divided the orbital space Ω into N/2 mutually disjoint subspaces Ω g , so each orbital belongs only to one subspace. Each subspace contains one orbital g below the Fermi level (N/2), and N g orbitals above it, which is reflected in additional sum rules for the ONs: p∈Ωg n p = 1; g = 1, 2, . . . , N/2.
Taking into account the spin, each subspace contains solely an electron pair, and the normalization condition for Γ (2 p n p = N) is automatically fulfilled. The energy (2) of PNOF5 can be then conveniently written as The first term of the energy draws the system as independent N/2 electron pairs, whereas the second term contains the contribution to the HF mean-field of the electrons belonging to different pairs. To go beyond the independent-pair approximation, let's maintain ∆ qp = 0 and consider nonzero the Π-elements between orbitals belonging to different subspaces [15]. From equation (4), note that provided the ∆ qp vanishes, |Π qp | ≤ Φ q Φ p with Φ q = n q h q . Assuming equality once again, the sign of Π qp remains undetermined, so there is a large number of possible combinations of signs that affect now the inter-pair interactions. Again, we need to solve this phase dilemma to propose a NOF, but this time we need to make a proper selection of the Π qp signs for orbitals q and p belonging to different pairs (subspaces). In contrast to the intra-pair interactions, there is no indication to determine the phase factor for the inter-pair Π qp .
Recently [15], the generalization of the sign convention adopted for Π g qp in equation (6), namely Π Φ qp = Φ q Φ p if q, p > N/2, and Π Φ qp = −Φ q Φ p otherwise, led to a new functional denoted as PNOF7. The resulting energy is It is obvious that a possible election that favors decreasing of the energy (9) is to consider all the phase factors negative, i.e., Π Φ qp = −Φ q Φ p . From now on we will denote by PNOF7(+) the functional that considers +1 the phase factors of Π Φ qp for q, p > N/2, whereas PNOF7(−) will be employed to denote the NOF (9) with these phases equal to −1.
Since we do not have an accurate functional like (5) that helps us to determine which is the best combination of signs for Π Φ qp , in the following sections, we analyze several examples with strong non-dynamic (static) electron correlation in order to make a proper phase choice after comparing with exact diagonalization calculations.

Hubbard model
The Hubbard model is an ideal candidate for the study of electron correlation due to its conceptual simplicity. The corresponding one-dimensional (1D) Hamiltonian reads as [43] where Greek indices µ and υ denote sites of the model, µ, υ indicates only near-neighbors interactions, t > 0 is the hopping parameter, σ = α, β; U is the on-site interelectron repulsion parameter, and n µ,σ = c † µ,σ c µ,σ where c † µ,σ (c µ,σ ) corresponds to fermionic creation(annihilation) operator. It is known that the HF approximation retrieves the exact solution for the 1D Hubbard model at halffilling if U = 0, where, in the site basis, all the possible states | − − , |− ↑ , |− ↓ , and | ↓↑ have the same weight. Conversely, in the U/t → ∞ limit the singly occupied states |− ↑ and |− ↓ appear uniquely, so the antiferromagnetic scheme is recovered and the model becomes equivalent to the spin-1/2 Heisenberg model. The performance of commonly used NOF approximations in the 1D Hubbard model with periodic boundary conditions has been recently studied [44,45], showing that the here presented PNOF7(+) is in good agreement with exact results for the Hubbard model at half-filling. Nevertheless, the amount of electron correlation recovered by PNOF7(+) for large systems is slightly less than for small systems. Since the NOF theory is a promising approach for large many-body systems, it is crucial to develop approximations that do not deteriorate as the size of the system increases. In the following, we show that a proper choice of inter-pair interaction signs prevents the accumulation of errors as the number of electron pairs in the system gets larger.
In Figure 1, we show the differences in E values with respect to the exact diagonalization (ED) results (∆E = E N OF − E ED ) obtained for the 8, 10, 12 and 14 sites systems for a range of U/t values in order to cover all correlation regimes (exact energy curves corresponding to these systems are included in the supplementary material). Exact results are computed using a modified version of the code developed by Knowles and Handy [46,47], whereas results for NOF approximations have been computed using DoNOF code developed by Piris and coworkers.
First, we observe that energies obtained by using the independent-pair model PNOF5, which is equivalent to a special case of an antisymmetrized product of strongly orthogonal geminals [24], systematically underestimate the correlation effects for all U/t values regardless of the number of sites of the system. Therefore, it is mandatory to consider the interactions between electron pairs to get an accurate description for the Hubbard model.
Considering nonzero inter-pair interactions by introducing the Π Φ qp term given in equation (9), the amount of electron correlation recovered in the region 0 < U/t < 10 is larger, but the behavior for large U/t values is rather similar to neglecting inter-pair interactions if PNOF7(+) is used. As illustrated in Figure 1, this issue can be properly solved simply by considering a different choice of the signs for Π Φ qp . Thus, PNOF7(−) significantly improves the performance of PNOF7(+) not only for any correlation regime, but also for any size of the system. While PNOF7(+) produces larger errors as the number of sites increases, the accuracy of PNOF7(−) is independent of the system size. In the region U/t 1, when the on-site electronic repulsion gets larger, the PNOF7(−) curve attaches to the exact curve giving an outstanding description of the asymptotic behavior. Hence, this approximation is able to reproduce the antiferromagnetic nature of the model in this region, in contrast to other approaches based on electron-pair states, such as AP1roG [48], which fail to describe weak orbital-pair correlations arising from singly occupied states in the strong correlation limit.

Hydrogen rings
In accordance with our previous benchmarking [44,45] and results showed in Figure 1, PNOF7(−) is the best approximation within NOF theory to study systems described by the Hubbard model. Within the limitations of the Hubbard model, the lack of long-range inter-electron interactions may be one of the most important. Therefore, in the following we focus on model systems with strong static electron correlation in order to examine if the conclusion obtained from the previous section still holds in presence of long-range interaction effects.
Let us consider a ring of hydrogen atoms and vary the number of atoms as done in the Hubbard model with the sites. We consider the non-relativistic many-electron Hamiltonian to describe these systems, i.e.  where the first term accounts for the inter-nuclear repulsion, the second term includes both the kinetic energy and the nuclear attraction, and the last term introduces Coulombic repulsion between electrons. Note that indices p, q, r, and t run over spatial orbitals, whereas τ and σ run over spin functions. This model may be the simplest example of strong electronic correlation in low dimensions since a multi-reference method is required to get an accurate description for R H−H = 2.0Å or larger bond distances, due to the strong correlation near the equilibrium geometry and at the dissociation limit [49]. The employed systems are illustrated in Figure 2 for chains of 2, 4, and 16 hydrogen atoms.
In Figure 3, we show the relative energies obtained by using PNOF7(+) and PNOF7(−) with respect to exact diagonalization, increasing the chain of hydrogen atoms from 2 to 16 at an internuclear distance of R H−H = 2.0Å. We use a minimal basis set in all the calculations. According to Figure 3, the results obtained employing PNOF7(+) show the same drawbacks already displayed for the Hubbard model. The relative errors shown by this approximation get larger as the size of the chain increases, so PNOF7(+) is not expected to give an accurate description of the electron correlation in the presence of many inter-pair interactions. In contrast, when we choose negatives, all the electron correlation functions Π Φ qp in equation (9), the relative errors with respect to the results of exact diagonalization remain equally small when increasing the number of hydrogens, as shown in Figure 3. Note that the accurate energy (5) is recovered for the twoelectron system by using either PNOF7(+) or PNOF7(−). The largest error obtained by using PNOF7(−) is below 0.007 Hartree (exact and approximated energies are given in supplementary material), so the latter is notably superior to PNOF7(+), and does not present any issue with the size of the system.

Closing remarks
In this work, we have presented a novel approach to tackle the phase dilemma in the context of natural orbital functional theory. The bottom-up method employed by Piris to develop approximate functionals does not require the use of a N-particle wavefunction and makes use of ensemble Nrepresentability conditions to get an explicit form of the functional. Nevertheless, there is still an indeterminacy with respect to the phase of the interaction between electron pairs with opposite spins. We have shown that this indefiniteness must be studied carefully, as it dramatically affects the performance of our approach.
For this purpose, we selected model systems with a strong static electron correlation, such as the onedimensional Hubbard model with periodic boundary conditions and molecular hydrogen rings. Despite of their simplicity, the Hubbard model and the hydrogen atom chain (located at R H−H = 2.0Å) present strong non-dynamic correlation effects, and can be viewed as benchmarking systems for testing multi-reference electronic structure methods. It has been demonstrated that the PNOF7 approach presented here captures the physics that appears in strongly correlated systems. After an adequate choice of sign factors for the inter-pair interactions, the so-called PNOF7 approximation gives a robust description of nondynamic correlation effects appearing in these systems, even in the region of strong correlation.
According to the results shown throughout the paper, the proper selection of phases amends the behavior of the functional when applying to large systems. Thus, the performance of the here presented method, denoted as PNOF7(−), does not deteriorate with the size of the system, so the latter could be used to study strongly correlated systems beyond small molecules, e.g. periodic polymers or heavy-element-containing molecules.