Benchmarking the variational reduced density matrix theory in the doubly-occupied configuration interaction space with integrable pairing models

The variational reduced density matrix theory has been recently applied with great success to models within the truncated doubly-occupied configuration interaction space, which corresponds to the seniority zero subspace. Conservation of the seniority quantum number restricts the Hamiltonians to be based on the SU(2) algebra. Among them there is a whole family of exactly solvable Richardson-Gaudin pairing Hamiltonians. We benchmark the variational theory against two different exactly solvable models, the Richardson-Gaudin-Kitaev and the reduced BCS Hamiltonians. We obtain exact numerical results for the so-called PQGT N-representability conditions in both cases for systems that go from 10 to 100 particles. However, when random single-particle energies as appropriate for small superconducting grains are considered, the exactness is lost but still a high accuracy is obtained.


Introduction
One of the main problems in many-body quantum mechanics, which includes condensed matter, nuclear physics, and quantum chemistry, is the so-called exponential wall problem, 1 namely, the exponential growth of the dimension of the Hilbert space with the number of particles composing the studied system. A complete diagonalization of the corresponding Hamiltonian in the many-particle space provides the exact answer but at a prohibitively expensive computational cost. Therefore, research efforts have been focused on the development of approximate methods capturing the relevant degrees of freedom present in the wavefunction at a feasible computational cost, i.e., with a polynomial increase.
A great variety of such approximate methods that have been developed over the years can be broadly classified into approximations that improve over a reference state and variational theories. In the former case, a standard approach is to start from a mean-field reference state and improve on this by adding perturbative corrections 2 or excitations with increasing complexity within Coupled Cluster Theory. 3,4 However, these methods break down in strong correlation regimes where multi-reference approximations are needed. New variational methods overcoming this issue were developed in the last decades. For instance, variational algorithms like tensor-network-state approaches, [5][6][7][8] variational Monte Carlo methods, [9][10][11][12] or stochastic techniques [13][14][15][16] can be made, in principle, as accurate as the exact diagonalizations while extending its computational limits beyond. Some of these many-body methods were recently benchmarked in the hydrogen chain. 17 A very different approach to tackle the exponential wall problem that is applicable to any correlation regime concentrates on the second-order reduced density matrix (2RDM), 18,19 while dispensing with the wavefunction altogether. The 2RDM is a much more compact object than the wavefunction and it holds all the necessary information to evaluate the expectation values of one-and two-particle observables of physical interest. As the energy of any pairwise-interacting system can be written as an exact but simple linear function of the 2RDM, it can be used to variationally optimize this matrix at polynomial cost. 20 This optimization should be constrained to the class of 2RDMs that can be derived from a wavefunction (or an ensemble of wavefunctions), the so-called N -representable 2RDMs. 21,22 Since the complete characterization of this class of 2RDMs is known to be a quantum Merlin Arthur (QMA) complete problem, 23 one has to use an incomplete set of necessary but not (in general) sufficient constraints on the 2RDM. The optimization thus finds a lower bound to the exact ground-state energy and an approximation to the exact ground-state 2RDM. Such an approach, known as the variational second-order reduced density matrix (v2RDM) method has been applied with different degrees of success in quantum-chemistry problems, 24-27 nuclear-physics, 28,29 and condensed-matter. [30][31][32] Recently, the computational efficiency of the v2RDM method has been substantially improved for systems whose states can be accurately described in terms of doubly-occupied single-particle states only. This lies at the heart of the doubly-occupied configuration interaction (DOCI) method, widely used in quantum chemistry to reduce the dimension of the configuration interaction Hilbert space. DOCI corresponds to the subspace of the Hilbert space of seniority zero, where the seniority quantum number 33 counts the number of unpaired particles. It has been recognized that the DOCI subspace captures most of the static correlations, serving as the first rung on a seniority ladder leading to the exact full configuration interaction (CI) solution. [34][35][36][37] The assumptions in DOCI drastically simplify the structure of the 2RDM 38,39 and reduce the scaling of the v2RDM method [40][41][42] while, expectedly, retain most of the correlation. Several applications of the v2RDM for seniority nonconserving Hamiltonians were already implemented and their accuracy tested against exact diagonalizations for small systems. 40,42,43 Here we will take advantage of the seniority-zero nature of DOCI space that restricts the Hamiltonians to be seniority conserving and therefore, to be based on the SU(2) algebra. An important class of SU (2) Hamiltonians are the pairing Hamiltonians, where the fundamental physics lies in the specific form of the paired states.
The quantum integrable and exactly solvable Richardson-Gaudin pairing models [44][45][46] are ideal Hamiltonians to test the performance of the v2RDM method within the DOCI space.
In this paper, we will benchmark the method for two different integrable Richardson-Gaudin models: the Richardson-Gaudin-Kitaev model 47 describing a chain of spinless fermions with p-wave pairing, and the constant pairing or reduced BCS Hamiltonian with uniform 48 and random single-particle energies. 49 We will also explore the behavior of the method for increasingly large systems addressing its extensivity properties

Theory
In second quantization, an N -particle Hamiltonian with pairwise interactions can be written where t and V are the one-body energy and the two-body interaction terms, respectively. c † i and c j are the standard fermion creation and annihilation operators in a given orthonormal single-particle basis {i, j, k, l, ...}.
According to Eq. (1), the ground-state energy can be expressed solely in terms of the second-order reduced density matrix, 2RDM, 2 Γ 18 and is the two-particle reduced Hamiltonian with |ψ the ground-state wavefunction and N the number of particles.
The idea behind the variational 2RDM methodology is to minimize the energy func- with H (p) ξ and E 0 (H ξ ) being the p-particle reduced Hamiltonian and the exact ground-state energy of H ξ , respectively. Unfortunately, this theorem cannot be used in practice since it would require knowledge of the ground-state energy of every possible p-body Hamiltonian H ξ . However, it can be relaxed using a set of Hamiltonians for which a lower bound for the ground-state energy is known. This is the case of the group of all semidefinite Hamiltonians, which are completely defined by its extreme elements yielding the well-known P, Q and G two-index N -representability conditions 22,52 on the 2 Γ matrix if B is restricted to the forms B = ij p ij c i c j , B = ij q ij c † i c † j , and B = ij g ij c † i c j , respectively. It has been shown that these conditions are the necessary and sufficient conditions to assure the N -representability for one-body Hamiltonians, 22,55 as well as for two-body Hamiltonians with an exact antisymmetric geminal power (AGP) ground state. 55, 56 We will demonstrate this last assertion analytically and numerically in Section 3.1 for the case of the Richardson-Gaudin-Kitaev Hamiltonian. 27,57,58 coming from the 3RDM on the 2 Γ matrix, respectively.

Hamiltonians of the class
As these conditions are in general necessary but not sufficient, the v2RDM will always find a lower bound to the exact ground-state energy and an approximation to the exact groundstate 2RDM.
In this work we will focus our attention on Hamiltonians with pairing interactions in the seniority zero subspace. Assuming time-reversal symmetry, the single-particle levels are doubly degenerate in the spin degree of freedom. The seniority quantum number classifies the Hilbert space into subspaces with a given number of singly occupied levels. The most general pairing Hamiltonian conserving seniority is where i are the energies of L doubly degenerate single-particle levels, N i = c † i c i + c † i cī is the number operator, and V ij is the pairing interaction. The (i,ī) pair defines the pairing scheme, which can involve two particles with either opposite spins (i ↑, i ↓), momenta (i, −i), or in general any classification of conjugate quantum numbers in doubly degenerate singleparticle levels. For these Hamiltonians the seniority number is an exact quantum number, as unpaired particles do not interact with the rest of the system and the pairing Hamiltonian does not allow for pair breaking. The Hamiltonian thus becomes block diagonal in sectors labeled by the seniority quantum number.
The pairing Hamiltonian (8) is based on the SU(2) pair algebra with generators and commutation relations We note that in the seniority zero subspace N i = 2K + i K − i and therefore, the Hamiltonian (8) can be written in terms of the ladder SU(2) operators as where J ij = δ ij i + V ij . The ground-state energy is thus given by where the P matrix is This matrix together with the D matrix define the seniority blocks of the 2 Γ matrix. Notice that the diagonal elements of both matrices are equal (D ii = P ii ). According to these definitions, it follows that the P and D matrices are hermitian and fulfill where M is the number of particle pairs in a system with L doubly degenerate single-particle levels and total ψ|K z |ψ = M − L 2 . The P, Q, G, T 1 and T 2 N -representability conditions can thus be written in terms of the seniority blocks of the 2RDM as, [38][39][40][41][42][43] • The P condition: • The Q condition: where • The G condition: g 0 where • The T 1 condition: • The T 2 condition: where where the symbol 0 denotes that a matrix is positive semidefinite.
The variational optimization of the 2RDM subject to conditions (15) for the PQGT conditions. This will allow us to treat without excessive computational efforts systems of sizes up to L = 100. In our numerical calculations we use the semidefinite programming algorithm (SDPA) code. 63,64 This code solves semidefinite problems at several precision levels by means of the Mehrotra-type predictor-corrector primal-dual interior-point method, providing ground-state energies and the corresponding 2RDM.
We programmed our v2RDM method as a dual problem in the SDPA code, which does not allow for the equality constraints (15)- (16). These are included by relaxing them into inequality constraints with a sufficiently small summation error δ. 27,65 In our work we have set δ = 10 −7 , which effectively fixes the precision of the ground-state energies.

Richardson-Gaudin integrable models
The Richardson-Gaudin (RG) models are based on a set of integrals of motion (IM) or quantum invariants that are linear and quadratic combinations of the generators of the SU(2) algebra. By requiring the IM to commute with the total spin operators K z = i K z i , the most general expression for the IM is where X ij and Z ij are antisymmetric matrices and G is the pairing strength. The operators R i must commute among themselves to constitute a set of IM. Imposing these conditions leads to two families of integrable models: 1. The hyperbolic or XXZ family 2. The rational or XXX family where the η s are an arbitrary set of real parameters.
The common eigenstates of IM (35) are determined by the solution of the set of M non-linear coupled RG equations with we can write the RG equations as a set of L coupled quadratic equations 66 in the Λ variables where C is a constant that depends on the Gaudin algebra, 0 for the rational family and −1 for the hyperbolic family. This new system of equations is free of the divergences that plague the original set of RG equations (38), and it can be solved numerically with the Levenberg-Marquardt algorithm. Once we have determined the set of Λ i for a particular eigenstate, the eigenvalues of the IM are If the Hamiltonian is an arbitrary linear combination of the IM,
The complete set of eigenstates in the seniority zero subspace is given by a product pair The PBCS or AGP wavefunction, being exact at G M R , will display important consequences for the v2RDM approach. As mentioned above, the PQG conditions are sufficient to produce the exact v2RDM result at this point. This statement can be independently proven starting from the set of killers of an AGP wavefunction such that from which the Moore-Read Hamiltonian derives as the positive semidefinite operator with 0 ground-state energy. In addition to the ground-state energies we will test another magnitude that characterizes the pair mixing across the Fermi level, the canonical gap defined as It turns out that ∆ c coincides with the BCS gap ∆ when it is evaluated with a number non-conserving BCS wavefunction. In this case the BCS gap equation reduces to As a function of G the system has a phase transition from a metallic state characterized by ∆ = 0 to a superconducting state with finite gap. The critical value of G is obtained from the gap equations as Even though BCS predicts a non-superconducting state for G < G c (∆ = 0), for correlated number conserving wavefunctions like PBCS or AGP the gap is always greater than zero. 73 We have now all the tools for testing the different variational approximations with the exact solution of the RGK model. We start with a system of L = 50 doubly degenerate levels at half filling corresponding to M = 25 fermion pairs. The size of the Hilbert space is  v2RDM. However, we should keep mind that this difference is positive for PBCS due to its Ritz variational character, while it is negative for v2RDM because it provides lower bounds.
The inset displays the behavior of the correlation energy, which stays flat for weak pairing, and starts to decrease linearly with G entering the superconducting region. The correlation energy is defined as where |ψ(0) is the ground state of the noninteracting Hamiltonian.  To ensure that the v2RDM method is extensible to systems of arbitrary sizes we show in Fig. 3 the comparison of the total ground-state energy under the PQGT conditions with the exact energy for systems with sizes ranging from L = 10 to 100 levels. To compute systems of such larger sizes we have relaxed the summation error to δ = 3 · 10 −7 , which is marginally lower that the previous computations. Our results show that the exact groundstate energies are numerically exact to the required precision independently of the system sizes. The relative energy errors are of the same order of magnitude taking into account that the correlation energy (inset of Fig. 1) increases by one order of magnitude along the horizontal axis.

The reduced BCS Hamiltonian
The reduced BCS or constant pairing Hamiltonian has been widely employed in condensed matter and nuclear physics to study superconducting properties of extensive as well as finite size systems in the BCS approximation. Few years after the celebrated BCS paper, Richardson solved this Hamiltonian exactly. 74 More recently, the exact solution has been generalized to families of exactly solvable pairing models. 44 In this subsection we will resort to the constant pairing Hamiltonian in the form used to describe ultrasmall superconducting grains 75 Richardson proposed a product pair ansatz for the exact eigenstates of the BCS Hamiltonian As in the RGK case, the pair energies, In small grains it is customary to assume equidistant levels and to express all quantities in units of the mean level spacing d, which in turn is inversely proportional to the volume of the grain. However, due to presence of disorder, the level spacing in small metallic grains follows a Wigner-Dyson distribution obtained from random matrix theory. We will take advantage of the two standard descriptions of small grains to benchmark the v2RDM. First, we will test it with uniformly distributed equidistant levels, and then investigate how robust is the method in the presence of random disorder.
In order to quantify pairing fluctuations around the Fermi level we make use of the canonical gap ∆ c For finite systems the BCS approximation has a metallic phase with no gap, and a superconducting phase with finite gap. The critical value of G is Since G c is a sensible value to assess the degree of superconducting correlations, we will study the BCS Hamiltonian for different system sizes as a function of G in units of G c .    Fig. 4). As in the RGK example, we have relaxed the summation error to δ = 3 · 10 −7 .
It is known that the energy levels of small metallic grains follow a Gaussian orthogonal ensemble distribution. For simplicity, most of the studies have been carried out assuming a uniform level spacing. However, the exact solution of the BCS Hamiltonian (54) is valid for arbitrary single-particle levels ε i . This feature has been exploited to study in an exact manner the interplay between randomness and interaction in the crossover from metal to superconductor as a function of the grain size. 49 Here, we will use this ability of the exact solution to test the robustness of the PQGT conditions against disorder in the single-particle levels spectrum. For each value of G/G c in Fig. 7 we generate 70 symmetric random matrices of size 2L × 2L. Upon diagonalization, we select the central L eigenvalues to avoid edge effects. In order to assure an average constant level spacing we rescale them From the exact solution, at this point the M pair energies E α converge to zero transforming the product of geminals (44) into the AGP (45). From the other side, starting with the AGP and making use of the killers we derived the Moore-Read Hamiltonian (49) that is contained in the G condition, and therefore the v2RDM with the PQG conditions should provide the exact solution. Fig. 1 gives the numerical proof of this statement in a highly non-trivial problem. This figure also shows that the variational method with the PQGT conditions gives the exact numerical ground-state energy from weak to strong pairing. Additional confirmation of the exactness of the PQGT conditions comes from the canonical gaps in Fig. 2, which also shows an exact value for PQG at the Moore-Read point. Similar results for the ground-state energies and gaps were obtained for the reduced BCS Hamiltonian with equidistant single-particle levels. We then tested the robustness of the PQGT Nrepresentability conditions against disorder in the single-particle levels as in the case of small metallic grains (see Fig. 7). Surprisingly, and even though the systems are always quantum integrable, the exactness of the numerical results was lost in the superconducting region (G > G c ). This fact might be explained by the complexity of the ground-state wavefunctions in most of the random instances, as can be deduced from the distribution of pair energies E α in the complex plane when the system enters the superconducting phase. However, relative errors of 10 −4 are still competitive with DMRG calculations 77 for equidistant levels, and with more recent approaches tested in the Richardson model for small size systems. 78,79 The exact solvability of these models allowed us to test the v2RDM method for large systems in order to asses its extensive properties. Fig. 3 and 6 demonstrate that the high accuracy of the PQGT is independent of the system size in the studied range from L = 10 to L = 100.
Before closing, we would like to point out that SU(2) Hamiltonians encompass the area of quantum magnetism with Heisenberg type Hamiltonians. The formalism developed in 42 and tested in this work could be directly applied to the study of spin systems. Due to the non-perturbative nature of v2RDM, it might be possible to describe with high accuracy exotic phases and quantum phase transitions.