Variational determination of ground and excited-state two-electron reduced density matrices in the doubly occupied configuration space: A dispersion operator approach

This work implements a variational determination of the elements of two-electron reduced density matrices corresponding to the ground and excited states of N -electron interacting systems based on the dispersion operator technique. The procedure extends the previously reported proposal [Nakata et al. , J. Chem. Phys. 125 , 244109 (2006)] to two-particle interaction Hamiltonians and N -representability conditions for the two-, three-, and four-particle reduced density matrices in the doubly occupied configuration interaction space. The treatment has been applied to describe electronic spectra using two benchmark exactly solvable pairing models: reduced Bardeen–Cooper–Schrieffer and Richardson–Gaudin–Kitaev Hamiltonians. The dispersion operator combined with N -representability conditions up to the four-particle reduced density matrices provides excellent results.


I. INTRODUCTION
Although the full configuration interaction (FCI) method constitutes the exact treatment to describe N-electron systems in a determined Hilbert subspace, its practical use is limited to studies of systems using small or medium basis sets due to its high computational cost.This shortcoming can be partially overcome by means of the determination of the corresponding two-particle reduced density matrix (2-RDM), which provides the information for the calculation of most relevant physical observables, including the energy.In the variational scheme, the 2-RDM elements are treated as parameters in energy minimization processes.The procedure, denominated variational determination of the two-particle reduced density matrix (v2RDM), requires the introduction of constraints in the The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpminimization algorithms so that the resulting 2-RDMs are limited to those arising from N-electron wave functions (or density matrices).The so-called P, Q, and G two-positivity (2-POS) conditions impose that the 2-RDM and its linearly associated twohole and particle-hole matrices turn out to be positive semidefinite.These three requirements, which are known as ensemble N-representability conditions of the 2-RDM, 1,2 constitute necessary although not sufficient conditions to guarantee the N-representability.In quantum chemistry, these conditions have provided results that could be considered as feasible when describing ground states of molecular systems at equilibrium geometries but present serious deficiencies otherwise.4][5] Although these conditions have allowed us to reach more accurate results, this treatment is limited to the description of the electronic ground state corresponding to the Hamiltonian of each system symmetry.
A different scenario appears when the v2RDM approach is applied to N-electron Hamiltonians in model systems used in nuclear and condensed matter physics. 6,7In fact, the description of strongly correlated systems or entangled states 8,9 requires the mandatory use of N-representability conditions involving the 3-and 4-RDMs, 4 otherwise no precise results are attained.Furthermore, the description of electronic spectra in these systems requires to identify discrete excited states, apart from the continuous ones.To this end, other variational approaches have been proposed, 10,11 which also require using higher order N-representability conditions than the 2-POS P, Q, and G ones. 12 Nakata et al. 10 reported suitable electronic spectra for one-particle Hamiltonians in model systems obtained by means of the so-called dispersion operator, whose expectation values are minimized according to several N-representability condition sets.The aim of this work is to go beyond, extending the dispersion operator technique to two-particle interaction Hamiltonian models and using a complete set of necessary N-representability conditions of the 2-, 3-, and 4-RDMs.This treatment has been accomplished with doubly occupied configuration interaction (DOCI) wave functions.The variational determination of the 2-RDM elements corresponding to those wave functions is known as v2RDM-DOCI.
This article is organized as follows: In Sec.II, we present the notation, theoretical aspects, and the features of the dispersion operator, which is the main tool used in this work.Using the formalism of hard-core bosons, 13 we formulate the dispersion operator as well as the elements of the 2-, 3-, and 4-RDMs for an N-electron Hamiltonian model of two-particle interactions in the DOCI space.This section also reports two exactly solvable models that have been chosen to perform numerical comparisons: reduced Bardeen-Cooper-Schrieffer (BCS) 14,15 and Richardson-Gaudin-Kitaev (RGK) 16,17 Hamiltonians.Section III reports the results obtained by these models along with the corresponding discussion, and Sec.IV gathers the concluding remarks.Finally, in the Appendix, we summarize the complete N-representability condition set used for the calculation of RDM elements.

II. THEORETICAL ASPECTS
The dispersion operator D( Ĥ, λ) corresponding to a general N-electron Hamiltonian Ĥ and a real parameter λ is defined as 10,18,19 where Î is the N-electron unit operator.The minimization of the expectation values of the dispersion operator D( Ĥ, λ), as a function of the parameter λ, leads to the determination of the zero points of that operator.These, in turn, yield the eigenvectors (or an arbitrary linear combination when degenerated) and eigenvalues of the Hamiltonian Ĥ.Consequently, finding out the zero points of the dispersion operator provides the solutions to the Schrödinger equation.
In this work, we will refer to N-electron Hamiltonians defined in the subspace of wave functions of seniority zero [20][21][22] or DOCI subspace, 23,24 which will be formulated by means of the SU(2) pair algebra.Within this approach, we will use the creation/annihilation hard-core boson operators b † i /bi, defined as b † i = a † i a † ī and bi = a¯iai, respectively, where a † i /ai are the standard fermionic particle creation/annihilation operators corresponding to an orthonormal single-particle basis set {i, j, k, l, . ..}.The notation (i, ī) defines a pairing scheme involving two particles with either opposite spins (i α , i β ), momenta (i, −i), or in general any pairing of conjugate quantum numbers in doubly degenerate single-particle levels and provides the formulation of the particle number operator as ni = b † i bi.As is well known, the hard-core boson operators satisfy the relations 13 [bi, b These operators allow us to express an N-electron Hamiltonian in the seniority zero subspace as 25 in which εi are the energies of the single-particle levels and Bij and Zij stand for the pairing and monopole interactions, respectively.Equation (3) allows us to relate N-electron Hamiltonians with hardcore boson operators and particle number operators.The commutation rules written down in formulas (2) enable us to express the dispersion operator of a seniority-zero N-electron Hamiltonian formulated by Eq. (3) as where The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp and M = N 2 is the number of hard-core bosons of an N-electron system.Expressions (4)-( 10) constitute a generalization of the dispersion operator formulation, which was reported in Ref. 10 for one-particle Hamiltonians.
7][28][29] This methodology provides optimized values for the elements of these RDMs arising from the wave function Ψ(N), 10 which FIG. 1.The v2RDM-DOCI and FCI energies for the constant pairing Hamiltonian (L = 8, M = 4) at G/Gc = 0.5 as functions of λ.The average energy is evaluated with the minimizer of the dispersion operator D( Ĥ; λ).The vertical dashed lines correspond to the exact excitation energies, and the v2RDM-DOCI results are computed under up to 4-POS conditions.

FIG. 2.
The expectation value of the dispersion operator D( Ĥ, λ) for the constant pairing Hamiltonian (L = 8, M = 4) as a function of λ.The different panels correspond to interaction strengths G/Gc = 0.5, 1.0, and 1.5, as indicated.The vertical dashed lines correspond to the exact excitation energies, and the v2RDM-DOCI results are computed under up to 4-POS conditions.

ARTICLE
scitation.org/journal/jcpthen can also be used to evaluate the energy and other observables of interest.The 2-, 3-, and 4-RDM matrix elements can be block-classified according to the seniority numbers of their creation or annihilation strings. 26In the case of the 2-RDM elements, these blocks are where Π and D elements correspond to the 0-and 2-seniority numbers, respectively.When the i and j indices are coincident, in which ρ i are the elements of the 1-RDM.The 3-RDM elements present two blocks of seniority numbers 1 and 3, which will be denoted respectively as and The 4-RDM presents the following three blocks: for the seniority numbers 0, 2, and 4, respectively.
In order to get physically meaningful results, these matrix elements must satisfy certain N-representability conditions.These conditions require that all those RDMs and their associated hole and particle-hole ones must be Hermitian, positive semidefinite, properly normalized, related by contraction mappings to lower-order RDMs and to fulfill determined consistency relationships between their elements.These requirements are imposed as variational constraints in the algorithms of the minimization processes.In the TABLE I. Absolute energy differences of the v2RDM-DOCI results with respect to the FCI ones for the ground (0) and excited (1, 2, and 3) states of the constant pairing Hamiltonian model at different interaction strengths, for a system with L = 8 and M = 4.The v2RDM-DOCI results are computed using up to 4-POS conditions.FCI energies are given in square brackets.The Journal of Chemical Physics

ARTICLE scitation.org/journal/jcp
Appendix, we summarize the complete set of N-representability conditions, which must be satisfied by the matrix elements of Eqs. ( 11)-( 17) in the numerical determinations performed in this work.1][32] The advantage of using the dispersion operator based treatment is that it allows us to describe excited states, in contrast to the more conventional minimization of the expectation value of the Hamiltonian Ĥ, which is limited to the description of ground states.
In Ref. 10, it was demonstrated that for one-body Hamiltonians, the enforcement of the 2-POS N-representability conditions warrants the positivity of ⟨Ψ(N)| D( Ĥ, λ)|Ψ(N)⟩.However, for two-body Hamiltonians, these conditions are insufficient and imposing the 4-POS ones becomes necessary.Furthermore, due to the lower-bound variational nature of the v2RDM-DOCI technique, the exact eigenvalues of the Hamiltonian, corresponding to zeroes of the dispersion operator, are always encountered as zeroes of the variationally determined counterpart.The converse is not true.Consequently, in order to avoid spurious solutions in the v2RDM-DOCI calculation, we will consider such a λ value as an excited state energy when no other zero points of the expectation value of the D( Ĥ, λ) operator exist in the neighborhood of that λ. 10 To deal with a tractable mathematical framework, in this work, we will use Hamiltonian models that have been widely used in condensed matter and nuclear physics for testing the accuracy of new theoretical treatments and approximations. 33The pairing Hamiltonians formulated by means of SU(2) algebra are relevant since they show the specific form of the paired states, revealing their fundamental physics aspects.The quantum integrable and exactly solvable Richardson-Gaudin pairing models [34][35][36] have proven to be suitable Hamiltonians to check the performance of the variational determination method of reduced density matrices within the DOCI space.Two different integrable Richardson-Gaudin models have been selected: the constant pairing or reduced BCS Hamiltonian 15,37,38 and the RGK model describing a chain of spinless fermions with p-wave pairing. 17,39he constant pairing Hamiltonian is formulated as TABLE II.Root-mean-square (rms) deviations for the Π (on the first row) and D (on the second row) seniority blocks of the 2-RDM resulting from the v2RDM-DOCI method with respect to that of the FCI one, for the ground (0) and excited (1, 2, and 3) states.The results correspond to the constant pairing Hamiltonian model at different interaction strengths, for a system with L = 8 and M = 4.The v2RDM-DOCI results are computed using up to 4-POS conditions.rms

G/Gc
The Journal of Chemical Physics

ARTICLE scitation.org/journal/jcp
where G is the strength of the infinite-range pairing interaction and εi are single-particle energies.We will only consider half-filled systems having the number of pairs M = L 2 (where L is the total number of levels) and equispaced single-particle energies εi = i L , i = 1, 2, . . ., L. This model predicts a metallic phase with no gap and a superconducting phase with finite gap.The critical value of G separating these phases is obtained from the gap equation in the limit of the gap going to zero and the chemical potential equal to as In the limit G → ∞, the exact ground state becomes a pair condensate denominated number projected BCS (PBCS) wave function in nuclear physics and antisymmetrized geminal power (AGP) in quantum chemistry. 40Therefore, the PQG conditions of the 2-RDMs are sufficient to reproduce the exact ground-state energies in this limit. 40ternatively, the RGK Hamiltonian model is formulated as where εi and G are again the single-particle energies and interaction strength, respectively; ηi = sin( i 2 ) and εi = η 2 i .Note that this value of εi describes the underlying fermion hopping between near neighbor sites in a 1D chain.If it is assumed antiperiodic boundary conditions and half-filled systems with the number of pairs M = L 2 , the allowed values for the indices i, j in formula (21) run in the set }.This Hamiltonian is a particular case of the hyperbolic or XXZ family of the Richardson-Gaudin models, 36 in which η's are an arbitrary set of real parameters and εi = η 2 i .The Richardson-Gaudin hyperbolic Hamiltonians ( 21) have an exact AGP ground state at a particular value of the interaction strength known as the Moore-Read (MR) point, [41][42][43]

ARTICLE scitation.org/journal/jcp
Therefore, the PQG conditions are sufficient to obtain the exact ground-state energies in this point. 39

III. RESULTS AND DISCUSSION
We have considered a dual semidefinite programming (SDP) problem formulation of the variational optimization of the 2-RDM in the DOCI space under the 2-, 3-, and 4-POS N-representability conditions.We have developed codes that formulate and solve the SDP problem adapted for the dispersion operator eigenvalue minimization satisfying the p-positivity conditions according to the sparse structure of the matrices arising from the seniority-zero wave functions.In our numerical calculations, we use the semidefinite programming algorithm (SDPA) code 44,45 in the SDPA-DD version. 46This code provides a highly accurate multiple-precision arithmetic procedure based on the Mehrotratype predictor-corrector primal-dual interior-point method, which allows us the evaluation of the ground-and excited-state energies, dispersion values, and their corresponding 2-RDMs.Furthermore, the SDPA code does not allow for the equality constraints such as those arising from the contraction mappings and consistency relations of the v2RDM method, as detailed in the Appendix.These are included by relaxing them into inequality constraints with a sufficiently small summation error of 10 −10 .The same precision limit has also been used for the ground and excited state energies, dispersion values, and RDMs in the convergence processes.This higher precision is a requirement arising from the squared-Hamiltonian dependence of the dispersion operator (1).
We first discuss the behavior of the average energy ⟨ Ĥ⟩ on the variationally obtained state at given λ. Figure 1 compares the values of ⟨ Ĥ⟩ as a function of λ obtained using the dispersor operator within the v2RDM-DOCI method with those from the FCI one.The comparison is performed for the constant pairing model at the interaction strength G/Gc = 0.5 for a system having L = 8 single-particle levels at half-filling.Analogously to that previously observed for one-body Hamiltonians by Nakata et al., 10 within the v2RDM-DOCI treatment, we may identify several states in both the low-lying and the high-lying part of the spectra of this two-body Hamiltonian.These states manifest themselves as plateaus of the average energy centered around specific values of λ corresponding to the Hamiltonian eigenvalues.Therefore, the dispersor operator within the v2RDM-DOCI approximation is effective at determining excitations.The values of the energy excitations are extracted from the location of the minima of the dispersor, as shown in Fig. 2 for G/Gc = 0.5, 1.0, and 1.5.Due to the symmetry of the constant TABLE III.Absolute energy differences of the v2RDM-DOCI results with respect to the FCI ones for the ground (0) and excited (1, 2, and 3) states of the RGK Hamiltonian model at different interaction strengths, for a system with L = 8 and M = 4.The v2RDM-DOCI results are computed using up to 4-POS conditions.FCI energies are given in square brackets.The Journal of Chemical Physics ARTICLE scitation.org/journal/jcppairing Hamiltonian with equispaced ϵi, at half-filling, the energy spectra for negative G = −|G| can be obtained by mirroring the results around the center of the spectra and shifting the results by a constant, (L + 1)M/L.As the interaction increases, the determination of the location of the minima around the center of the spectra becomes increasingly difficult as the expectation value of D( Ĥ, λ) develops continuous regions in λ where it attains a zero value.Therefore, the central part of the spectra in these cases cannot be accounted for under up to 4-POS conditions and stronger N-representability conditions may be needed.

State
To quantitatively assess the quality of our results for the lowlying part of the spectra, in Table I, we collect the absolute values of the differences found between the energies arising from the v2RDM-DOCI method using up to 4-POS N-representability conditions and those from the FCI one.These results correspond to the ground and first three excited states of a system having L = 8 single-particle levels at half-filling, M = 4.They have been obtained for the sequence of interaction strengths lying in the interval G/Gc ∈ [−2.0, 2.0].As can be observed in Table I, the energies obtained from the v2RDM-DOCI method, for both attractive (G > 0), and repulsive (G < 0) branches, turn out to be very close to their corresponding FCI ones, for the ground state as well as for the excited states.However, slightly larger energy differences appear for high |G/Gc| strengths.The attractive branch of the third excited state presents spurious solutions for G/Gc ≥ 1.25 values, which correspond to situations in which zero points of the expectation value of the D( Ĥ, λ) operator exist in the neighborhood of the exact excited state energy, as mentioned above.The quality of the results is corroborated by evaluating the root-mean-square (rms) deviations of the 2-RDM elements (for both Π and D seniority blocks) arising from the v2RDM-DOCI method with respect to the FCI one.The results are gathered in Table II and show that the largest rms values correspond to the highest |G/Gc| strengths, mainly in the attractive branch of the excited states.
We have also studied the system constituted by L = 8 singleparticle levels at half-filling M = 4 within the RGK Hamiltonian model [Eq.( 21)], using identical N-representability conditions to the former model.The results for the expectation value of D( Ĥ, λ) as function of λ are shown in Fig. 3.At variance to the constant pairing model studied above, the energy spectra of the RGK Hamiltonian model at half-filling is not symmetric around its center, and thus, we have chosen to display the spectra for both positive and negative values of G.In this case, the dispersion operator method based on the v2RDM-DOCI allows us to identify a higher number of excited TABLE IV.Root-mean-square (rms) deviations for the Π (on the first row) and D (on the second row) seniority blocks of the 2-RDM resulting from the v2RDM-DOCI method with respect to that of the FCI one, for the ground (0) and excited (1, 2, and 3) states.The results correspond to the RGK Hamiltonian model at different interaction strengths, for a system with L = 8 and M = 4.The v2RDM-DOCI results are computed using up to 4-POS conditions.rms    states.Furthermore, the method captures better a broader region in the low-energy part of the spectra for repulsive interactions and in the high-energy part for attractive ones.In Tables III and IV, we report the results corresponding to energies and 2-RDM elements, respectively, in the interaction strength interval G/GMR ∈ [−2.0, 2.0].Again, no significant differences have been found between the FCI results and those obtained from the v2RDM-DOCI method, neither in energies nor in 2-RDM element values.Furthermore, the RGK Hamiltonian model leads to slightly more precise results than the constant pairing model.Finally, in order to compare the ability of our procedure to describe systems of different sizes, we have computed the TABLE VI.Root-mean-square (rms) deviations for the Π (on the first row) and D (on the second row) seniority blocks of the 2-RDM resulting from the v2RDM-DOCI method with respect to that of the FCI one, for the ground (0) and excited (1, 2, and 3) states.The results correspond to the RGK Hamiltonian model at G/G MR = 1.75 interaction strength, for systems with L = 6, 8, 10, and 12 at half-filling.The v2RDM-DOCI results are computed using up to 4-POS conditions.lowest lying and ground-state energies using the dispersion operator method for systems with L = 6, 8, 10, and 12 at half-filling arising from the RGK Hamiltonian model at G/GMR = 1.75.Tables V and VI show the results corresponding to the energies of the systems and the rms deviations of their 2-RDM element values with respect to the FCI ones, respectively.As can be observed in Tables V and VI, the precision level slightly decreases as the system size increases.All these results confirm that our methodology is able to describe the low-lying energy states of strongly correlated systems in the seniority-zero subspace.The accuracy of the numerical values predicted by the v2RDM-DOCI method, for both constant paring and RGK Hamiltonian models, must be attributed to the suitability of the N-representability conditions used in this work (see the Appendix).

IV. CONCLUDING REMARKS AND PERSPECTIVES
In this work, we have reported algorithms that provide the description of electronic spectra of two-particle interaction N-electron Hamiltonian models by means of the dispersion operator technique.This study constitutes an extension of the previously reported procedures limited to the one-body Hamiltonian case.The search for zero points vanishing the expectation value of the dispersion operator leads to the variational determination of the 2-, 3-, and 4-RDM elements of the ground and excited states.This task has been performed in the DOCI space within the v2RDM framework by means of codes that formulate and solve the SDP problem adapted for the dispersion operator eigenvalue minimization, imposing a complete set of 4-POS N-representability conditions.The numerical results have been obtained with this procedure for systems at half-filling, possessing L = 8 and L = 6, 8, 10, and 12 electrons, using the constant pairing and RGK Hamiltonian models, respectively.The almost coincidence of those results with their counterparts from the FCI method warranties the suitability of the N-representability conditions proposed in this work.In order to know all abilities of this methodology, we are currently exploring the description of electronic spectra in closed-and open-shell molecules.

ACKNOWLEDGMENTS
The authors acknowledge the financial support from the Universidad de Buenos Aires (Grant No. 20020190100214BA); the Consejo Nacional de Investigaciones Científicas y Técnicas (Grant Nos.PIP 11220130100377CO, PIP 11220130100311CO, and 11220150100442CO); and the Agencia Nacional de Promoción Científica y Tecnológica, Argentina (Grant No. PICT-201-0381).E.R. acknowledges the support from the Consejo Nacional de Investigaciones Científicas y Técnicas.J.D. and A.R.-G.acknowledge the financial support from the Spanish Ministerio de Ciencia e Innovación and the European regional development fund (FEDER), Project No. PGC2018-094180-B-I00.

APPENDIX: N-REPRESENTABILITY CONDITIONS FOR 2-, 3-, AND 4-RDMs ARISING FROM DOCI WAVE FUNCTIONS
The Hermitian matrices ( 11)-( 17) must satisfy the following normalization, contraction, and consistency conditions: and the 3-F one is The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp

)FIG. 3 .
FIG. 3. The expectation value of the dispersion operator D( Ĥ, λ) for the RGK Hamiltonian model (L = 8, M = 4) as a function of λ.The panels correspond to interaction strengths G/G MR = −1.5, −0.5, 0.5, and 1.5, as indicated.The vertical dashed lines correspond to the exact excitation energies, and the v2RDM-DOCI results are computed under up to 4-POS conditions.

TABLE V .
Absolute energy differences of the v2RDM-DOCI results with respect to the FCI ones for the ground (0) and excited (1, 2, and 3) states of the RGK Hamiltonian model at G/G MR = 1.75 interaction strength, for systems with L = 6, 8, 10, and 12 at half-filling.The v2RDM-DOCI results are computed using up to 4-POS conditions.FCI energies are given in square brackets.
with X = Π or D.