Variational theory combining number-projected BCS and coupled-cluster doubles
VVariational theory combining number-projected BCS and coupled-cluster doubles
V.V. Baran , ∗ and J. Dukelsky † Faculty of Physics, University of Bucharest, 405 Atomi¸stilor,POB MG-11, Bucharest-M˘agurele, RO-077125, Romania ”Horia Hulubei” National Institute of Physics and Nuclear Engineering,30 Reactorului, RO-077125, Bucharest-M˘agurele, Romania Instituto de Estructura de la Materia, CSIC, Serrano 123, E-28006 Madrid, Spain
The ground state pairing correlations in finite fermionic systems are described with a high degreeof accuracy within a variational approach based on a combined coupled-cluster and particle-number-projected BCS ansatz. The flexibility of this symmetry-preserving wavefunction enables a unifiedpicture valid from weak to strong coupling, both in small and large systems. The present variationalapproach consistently yields an energy upper bound while operating at the same level of precisionof the non-variational particle-number projected Bogoliubov-coupled-cluster theory [Phys. Rev. C99, 044301 (2019)].
I. INTRODUCTION
Pairing Hamiltonians are ubiquitous in quantummany-body physics. Starting from their variationaltreatment in the microscopic theory of superconduc-tivity given by Bardeen, Cooper and Schrieffer (BCS)[1], they were soon exported to nuclear physics forthe descripton of the large gaps observed in even-evennuclei [2]. However, the violation of particle num-ber in the BCS theory, which is negligible for macro-scopic systems, represents a major drawback when ap-plied to finite systems. Therefore, techniques to im-plement number projection on top of the BCS wave-function (PBCS) were developed in nuclear structure[3, 4] and more recently in quantum chemistry [5, 6]where the PBCS wavefunction is known as the anti-symmetrized geminal power (AGP). PBCS improvesover the BCS theory in finite systems, specially in thestrong coupling limit where superconductivity is wellestablished, but it still fails in the weak coupling limitdominated by pairing fluctuations, and all along thetransitional region [7, 8]. This fact was made evidentin ultrasmall superconducting grains, where PBCSpredicted an abrupt metal-superconductor transitionas a function of the grain size [9] while the exact so-lution showed a smooth crossover dominated by largefluctuations [10]. It was precisely in the field of ultra-small superconducting grains that the exact solutionof the constant pairing Hamiltonian given by Richard-son in the sixties [11] was recovered [12] and inten-sively used as a natural benchmark model for super-conducting theories beyond BCS [13–21].In the extreme weak coupling limit pair coupled-cluster doubles (pCCD) describes correctly the pair-ing fluctuating regime but it quickly overbinds dueto the non-variational character of the theory basedon the left projection to a Hartree-Fock Slater deter-minant [8, 22, 23]. It seems, therefore, that a com- ∗ vvbaran@fizica.unibuc.ro † [email protected] bination of pCCD and quasiparticle BCS would beable to approach both limits correctly. Indeed, theextension of pCCD to BCS quasiparticles, the BCS-CCD method [15], gave the correct behavior in theweak coupling limit, but still suffers from large devia-tions ( ∼
10 % for sizes of ∼
100 particles) across thetransition region. There were attempts to interpolatebetween pCCD and PBCS [8, 16] or to diagonalizethe pairing Hamiltonian in a subspace defined by thereference PBCS state and different two and four num-ber projected quasiparticle states [24]. Other possibleways to add correlations to the PBCS state were ex-plored in [25, 26].Perhaps, the most successful theory beyond BCS-CCD is simply its number projected version coined asparticle-number projected BCS coupled-cluster dou-bles (PBCS-CCD) theory [27]. The theory is not Ritzvariational and therefore, it cannot assure an upperbound for the ground state energy. However it has anaffordable computational cost that scales polynomi-ally with the system size. Moreover, it gives excellentnumerical results both in the weak and strong cou-pling limits, as well as in the transitional region.The aim of our paper is to design a Ritz varia-tional method with a trial wavefunction that com-bines pCCD and PBCS and produces numerical re-sults with the same level of precision of PBCS-CCDand with a similar computational cost. We will bench-mark our variational theory with the exact solution ofthe Richardson model [28] and with PBCS-CCD re-sults [27] where available.
II. THEORETICAL BACKGROUNDA. Pairing Hamiltonian
We consider the generic pairing Hamiltonian H = L (cid:88) i =1 (cid:15) i ( c † i c i + c † ¯ i c ¯ i ) + L (cid:88) i,j =1 V i,j c † i c † ¯ i c ¯ j c j , (1) a r X i v : . [ nu c l - t h ] F e b where i and ¯ i indicate one of the L pairs of conjugateddegenerate single particle levels with energy (cid:15) i = (cid:15) ¯ i .The Hamiltonian (1) preserves seniority, and forsimplicity we restrict ourselves to the seniority zero( v = 0) subspace . Then, in the absence of the inter-action term the ground state is given by the Hartree-Fock product state | HF (cid:105) = M (cid:89) i =1 c † i c † ¯ i | (cid:105) , (2)where M is the number of pairs in the system.In preparation for the discussion that follows, wepass to the particle-hole (ph) representation and de-note the particle levels with p > M and the hole levelswith h ≤ M . We define the particle and hole pair andnumber operators P † p = c † p c † ¯ p , P p = c ¯ p c p , N p = c † p c p + c † ¯ p c ¯ p ,P † h = c ¯ h c h , P h = c † h c † ¯ h , N h = 2 − c † h c h − c † ¯ h c ¯ h , (3)such that the Hartree-Fock state (2) is the vacuumto the P and N operators, P p | HF (cid:105) = P h | HF (cid:105) = 0, N p | HF (cid:105) = N h | HF (cid:105) = 0. The pairing Hamiltonian (1)is then expressed as H = E HF + L (cid:88) p = M +1 (cid:15) p N p + M (cid:88) h =1 ( − (cid:15) h − V hh ) N h + L (cid:88) p,p (cid:48) = M +1 V pp (cid:48) P † p P p (cid:48) + M (cid:88) h,h (cid:48) =1 V hh (cid:48) P † h P h (cid:48) + L (cid:88) p = M +1 M (cid:88) h =1 V ph (cid:16) P † p P † h + P h P p (cid:17) , (4)with E HF = (cid:104) HF | H | HF (cid:105) = (cid:80) Mh =1 (2 (cid:15) h + V hh ) beingthe energy of the Hartree-Fock state (2). B. Mean-field theory and its symmetryrestoration
The standard description of the pairing correlationsinduced by the Hamiltonian (1) is given within theBCS approximation [1] in terms of the pair condensate | BCS (cid:105) = exp[Γ † ( x )] | (cid:105) , Γ † ( x ) ≡ L (cid:88) i =1 x i c † i c † ¯ i , (5)which explicitly breaks the U (1) gauge symmetryassociated with particle number conservation. Formacroscopic systems the symmetry broken picture isexact and the particle number fluctuations are negli-gible. However, for finite systems like atomic nucleior small superconducting grains one speaks only of obscured or emergent symmetry breaking [29–32], inwhich case the quantum fluctuations inevitably liftany degeneracy associated with the broken symmetry. Much effort is thus devoted to restore the symmetryof the mean-field ansatz with the help of projectiontechniques [33–39]. In the BCS case upon particlenumber restoration we obtain the so-called number-projected BCS (PBCS) [3] or antisymmetrized gemi-nal power (AGP) in the context of quantum chemistry[40] | PBCS( x ) (cid:105) = P M | BCS (cid:105) = 1 M ! [Γ † ( x )] M | (cid:105) = 12 π (cid:90) π d θ e − iθM exp[Γ † ( e iθ x )] , (6)where P M is the projector onto the state of M pairs.While PBCS describes well the properties of super-fluid nuclei with a small number of valence nucleons[36], it cannot account for the weak pairing correla-tions that develop within larger spaces, e.g. thoseconsidered in the large-scale energy density functionaltreatments of finite nuclei or in the study of small su-perconducting grains [41]. For a working descriptionof the weak pairing regime one usually turns to RPA[22] or coupled-cluster [8, 15, 23] approaches. Gener-alizations of the PBCS ansatz have also been consid-ered based on its structural similarity with a particu-lar coupled-cluster ansatz [8]. Specifically, the PBCSrepresentation in the ph-basis is obtained as | PBCS (cid:105) ∝ M (cid:88) (cid:96) =0 (cid:96) ! [Γ † P ( x )Γ † H (1 /x )] (cid:96) | HF (cid:105) = 12 π (cid:90) π d θ exp[Γ † P ( e iθ x )] exp[Γ † H ( e − iθ /x )] | HF (cid:105) , (7)in terms of the particle and hole components of thecollective pairsΓ † P ( x ) = L (cid:88) p = M +1 x p P † p , Γ † H ( x ) = M (cid:88) h =1 x h P † h . (8)The structure of the PBCS state is then defined bythe inverse squared factorials appearing as expansioncoefficients in the collective ph-pair basis. C. Coupled-cluster theory
Analogously, by using a slightly modified expansioninvolving plain factorials one obtains a separable paircoupled-cluster doubles variational ansatz (vCCD sep ) | vCCD sep (cid:105) = M (cid:88) (cid:96) =0 (cid:96) ! [Γ † P ( y )Γ † H ( y )] (cid:96) | HF (cid:105) = exp[Γ † P ( y )Γ † H ( y )] | HF (cid:105) = exp (cid:88) p,h y p y h P † p P † h | HF (cid:105) . (9)This is a particular case of a pair coupled-cluster dou-bles variational wavefunction (vCCD) | vCCD (cid:105) = exp (cid:88) p,h z ph P † p P † h | HF (cid:105) , (10)involving the most general double excitations that donot break pairs through a fully non-separable struc-ture matrix z ph [22]. The separable case vCCD sep isthus recovered for z ph = y p y h .On the one hand, the full freedom in the structurematrix of the CCD ansatz allows for an excellent de-scription of the weak pairing regime (see also Fig. 1below), due to the additional four-body correlationsaccounted for relative to the separable case. Note,for example, that the Richardson solution involvingcomplex-conjugated pairs may be expressed exactlyas a product of four-body quartet structures [42]. Onthe other hand, the computational complexity of eval-uating operator matrix elements exactly in the non-separable case grows exponentially with the size of thesystem. The usual approximation of coupled clusterinvolves a left projection onto the subspaces of zeroand two ph pairs, thus breaking the Ritz variationalprinciple. Furthermore, even with a fully nonsepara-ble structure matrix, the validity of the vCCD ansatzbreaks down around the critical value of the pairingstrenght (see also Fig. 1 below). The choice of a BCSmean-field reference state does improve on this aspectat the cost of effectively breaking the particle numbersymmetry [15].Given the success of vCCD in the weak pairingregime and that of PBCS in the strong pairing regime,it is then natural to combine them for a precise unifieddescription of all regimes. D. Combining coupled-cluster andsymmetry-restored mean-field theories
The symmetry restoration of broken-symmetrycoupled-cluster theories has been only recently consid-ered [27, 43–45]. For the schematic pairing Hamilto-nian with the breaking and restoration of the particlenumber symmetry, the so-called particle-number pro-jected Bogoliubov-coupled-cluster theory yields highlyaccurate results [18]. In this context, a set of differen-tial equations is set up for obtaining the gauge-angle-dependent excitation operator. The practical needfor truncating this set of ODEs implies however anapproximate action of the projection operator and aviolation of the Ritz variational principle.It is the purpose of this work to explore a physicallytransparent alternative in the form of a variationalCCD-PBCS combined approach. Ideally the groundstate would involve the CCD excitations built directly on top of the particle-number-projected BCS as | vCCD-PBCS (cid:105) = exp (cid:88) p,h z ph P † p P † h | PBCS( x ) (cid:105) , (11)where the L mixing amplitudes x i of the PBCS stateand the M ( L − M ) structure matrix elements z ph ofthe CCD excitations are to be treated as free varia-tional parameters. Indeed, the energy minimizationprocedure E gs = min ψ (cid:104) ψ | H | ψ (cid:105)(cid:104) ψ | ψ (cid:105) (12)yields extremely precise results for the few small sys-tems where this ansatz may actually be applied.Computational access to larger systems is enabledupon various simplifications. Let us consider first theseparability approximation z ph = y p y h for the struc-ture matrix of the CCD excitations in Eq. (11), yield-ing the vCCD sep PBCS ansatz | vCCD sep PBCS (cid:105) = exp[Γ † P ( y )Γ † H ( y )] | PBCS( x ) (cid:105) , (13)which involves 2 L independent variational parameters x i , y i . Remarkably, while the PBCS ansatz (6) andthe vCCD sep ansatz (9) each individually fails in theweak pairing regime due to their separable structurematrices, the combined vCCD sep PBCS wavefunctionwill turn out to be quite accurate. This may be easilyunderstood by noting that for weak pairing the aboveansatz effectively involves a non-separable structurematrix in the form | vCCD sep PBCS (cid:105) ≈ [1 + (cid:88) p,h ( x p x h + y p y h ) P † p P † h ] | HF (cid:105) , (14)obtained after taking into account the PBCS (6)and vCCD sep (9) expansions and redefining x h → /x h in (7). The accuracy of this doubly-separable2 L -dimensional parametrization of the fully non-separable M ( L − M )-dimensional structure matrixnaturally degrades with increasing system size, theactual rate being numerically determined in the nextsection.An alternative approximation scheme would theninvolve a fully non-separable CCD excitations limitedto a relatively small finite window around the Fermilevel. We thus consider the combination | vCCD ( w ) PBCS( x, z ) (cid:105) = vCCD ( w ) ( z ) | PBCS( x ) (cid:105) (15)withvCCD ( w ) ( z ) = exp M + w/ (cid:88) p = M +1 M (cid:88) h = M +1 − w/ z ph P † p P † h , (16)where w denotes the size of the truncation window. Aswill be detailed in the next section, the quality of theresults obtained within this approach will turn out tobe inferior to that of the above vCCD sep PBCS ansatzof Eq. (13) involving two sets of global parameters.Overall the optimal compromise between the com-putational complexity and the accuracy of the resultsis found by combining the above two approximationschemes (9) and (15) into the ansatz | vCCD ( w )sep PBCS( x, y, z ) (cid:105) = vCCD ( w )sep ( y, z ) | PBCS( x ) (cid:105) , vCCD ( w )sep = vCCD ( w ) ( z ) vCCD sep ( y ) , (17)involving a total of 2 L + ( w/ free variational pa-rameters and leading to highly accurate results com-parable to those of Ref. [18], to be discussed lateron. Next, we shortly review the actual computationalstrategy for the above mentioned wavefunctions. E. Computational aspects
In this section we propose a novel efficient algorithmfor the evaluation of expectation values on the com-bined vCCD ( w )sep PBCS wavefunction (17). It is basedon a representation of the PBCS (7), vCCD sep (19)and vCCD ( w ) (16) terms as disentangled particle andhole gauge-angle-rotated BCS states.We start from the discrete exact representation forthe particle-number projection operation and decom-pose the PBCS ansatz using the definitions of Eqs.(7) and (8) | PBCS (cid:105) = L (cid:88) n =0 exp[Γ † H ( e − iθ n x )] exp[Γ † P ( e iθ n x )] | HF (cid:105) , (18)where θ n = 2 πn/ ( L + 1). Note that we neglect theirrelevant constant normalization factors throughoutthis section.The collective pairs (for the particle and hole sub-spaces) appearing in the vCCD sep operatorvCCD sep = M (cid:88) (cid:96) =0 (cid:96) ! [Γ † P ( y )Γ † H ( y )] (cid:96) (19)are also expanded as superpositions of gauge-angle-rotated BCS operators[Γ † P ( y )] (cid:96) (cid:96) ! = L − M (cid:88) k =0 e − iφ k (cid:96) exp[Γ † P ( e iφ k y )] , [Γ † H ( y )] (cid:96) (cid:96) ! = M (cid:88) k =0 e − iϕ k (cid:96) exp[Γ † H ( e iϕ k y )] , (20)with φ k = 2 πk/ ( L − M + 1) and ϕ k = 2 πk/ ( M + 1).Finally we obtain the representation | vCCD sep PBCS (cid:105) = L − M (cid:88) k p =0 M (cid:88) k h =0 L (cid:88) n =0 g k p k h × BCS P ( k p , n ) BCS H ( k h , n ) | HF (cid:105) (21) in terms of the gauge-angle-rotated BCS operatorsBCS P ( k p , n ) = exp[Γ † P ( e iθ n x + e iφ kp y )] , BCS H ( k h , n ) = exp[Γ † H ( e − iθ n x + e iφ kp y )] , (22)with expansion coefficients g k p k h = M (cid:88) (cid:96) =0 (cid:96) ! exp[ − i ( φ k p + ϕ k h ) (cid:96) ] . (23)We are interested in the expectation values ofgeneric particle-number-conserving ph-factorized op-erators O = O ( P ) O ( H ) on the vCCD sep PBCS state(21). As the particle-number projection operationneeds only to be performed once on the mean-fieldBCS wavefunction, we obtain (cid:104)O(cid:105) = (cid:104) vCCD sep BCS |O ( P ) O ( H ) | vCCD sep PBCS (cid:105) = L − M (cid:88) k p ,k (cid:48) p =0 M (cid:88) k h ,k (cid:48) h =0 g ∗ k p k h g k (cid:48) p k (cid:48) h L (cid:88) n =0 O ( P ) k p ,k (cid:48) p ,n O ( H ) k h ,k (cid:48) h ,n , (24)involving the matrix elements between particle andhole BCS states (22) O ( P ) k p ,k (cid:48) p ,n = (cid:104) BCS P ( k p , | O ( P ) | BCS P ( k (cid:48) p , n ) (cid:105) , O ( H ) k h ,k (cid:48) h ,n = (cid:104) BCS H ( k h , | O ( H ) | BCS H ( k (cid:48) h , n ) (cid:105) . (25)By considering all terms in Eq. (4), the expression(24) may be directly employed to compute the energyexpectation value E = (cid:104) H (cid:105) / (cid:104) I (cid:105) on the vCCD sep PBCSwavefunction.To include the effects of the non-separable vCCD ( w ) excitations of Eq. (17) we use a Hubbard-Stratonovichtransformation [46] and pass to a disentangled BCSrepresentation of the vCCD operatorvCCD( z ) = exp (cid:88) p,h z ph P † p P † h = (cid:90) d L ξ exp (cid:18) − ξ T Z − ξ (cid:19) exp[Γ † P ( ξ ) + Γ † H ( ξ )] , (26)with Z = (cid:18) zz T (cid:19) . (27)This BCS representation allows us to generalize theabove Eq. (24) to the expectation values of ph-factorized operators on the vCCD ( w )sep PBCS wavefunc-tion (17). In practice the matrix elements on the re-sulting BCS wavefunctions which generalize Eq. (25)are to be treated as polynomials in the integrationvariables ξ . Only specific terms corresponding to thevarious nonzero Wick contractions are to be selectedaccording to (cid:90) d L ξ exp (cid:18) − ξ T Z − ξ (cid:19) ξ i ξ j . . . ξ k ξ l = (cid:88) Wick Z ab . . . Z cd , (28)the set of indices { a, b, ..., c, d } representing a permu-tation of { i, j, ..., k, l } [47].With the non-separable CCD excitations restrictedto a small window of w levels around the Fermi levelas in Eq. (15), the matrix elements of the relevant op-erators on generic vCCD ( w ) BCS states may be com-puted analytically (the notebook in Cadabra2 [48] isavailable upon request from the authors). Their sub-sequent coupling to the separable subspace (i.e. thecomplement of the w -level window) leads to the finalvCCD ( w )sep PBCS results presented in the next section.
III. NUMERICAL RESULTS
In this section we benchmark the variational calcu-lations for the various wavefunctions presented aboveagainst the exact solution [28] for a constant pair-ing Hamiltonian ( V i,j = − G in Eq. (1)) of equallyspaced single particle levels with energies (cid:15) k = k(cid:15), k =1 , . . . , L . In order to avoid the singularities in theequations that solve exactly the picket fence modelwe use the same algorithm proposed by Richardson in[28].The numerical code used to compute the expec-tation value of the Hamiltonian (4) on the variousderivatives of the vCCD-PBCS wavefunction is freelyavailable upon request from the authors. The min-imization procedure for the energy function (12) isperformed using the e04ucf routine of the NAG library[49].We present in Figure (1) the errors for the correla-tion energy E c = (cid:104) ψ | H | ψ (cid:105)(cid:104) ψ | ψ (cid:105) − E HF , (29)relative to its exact value for L = 12 at half fill-ing. On the one hand, notice how PBCS (6) andCCD sep (9) are limited in the weak pairing regimeby their common separable structure, both reducingto (1 + Γ † P Γ † H ) | HF (cid:105) at G/(cid:15) (cid:28)
1. On the other hand,the full generality of the structure matrix of the vCCDwavefunction (10) allows it to capture precisely all thecorrelations in this regime. However, beyond the crit-ical value G cr of the HF to BCS transition the vCCDansatz quickly loses its ability to describe the strongercorrelations, becoming indistinguishable from its sep-arable version vCCD sep for G/(cid:15) > sep
PBCS G cr L =
12, M = sep PBCS + vCCD sep vCCD sep + PBCSvCCD sep
PBCS G / ϵ - E c / E c , e x ac t ( % ) FIG. 1. Error in the correlation energy (29) relative to itsexact value (in percentages) versus the pairing strength G (in units of the level spacing (cid:15) ) for the fully variationalvCCD (10), PBCS (6), and vCCD sep PBCS (13) wavefunc-tions in the case of an L = 12 level system at half fill-ing. The “vCCD sep +PBCS” procedure indicated by a reddashed line involves an energy minimization with respectto the vCCD sep amplitudes in the presence of the fixedoptimal PBCS reference. Within the “PBCS+vCCD sep ”procedure indicated by a green dashed line the PBCS am-plitudes are varied in the presence of the fixed optimalvCCD sep structure. wavefunction (13) leads to very accurate energeticsacross all regimes with the relative errors in the cor-relation energy not exceeding 3 · − . Note how-ever that it is essential to enable all parameters tovary freely in the minimization process as to retainthe full flexibility of the wavefunction and thus re-cover all available dynamical correlations. Indeed, byconsidering a variational CCD sep on top of the frozenoptimal PBCS reference (or viceversa) we find signif-icant improvements only in the weak pairing regime,as indicated by the dashed lines in Fig. (1). In partic-ular, by varying in this way just the PBCS amplitudeswithin vCCD sep PBCS while keeping fixed the optimalCCD sep structure (computed beforehand) we find noimprovement over the full vCCD results for
G > G cr .At G/(cid:15) (cid:28)
1, the results differ due to the limitationsdiscussed around Eq. (14).The relative errors for the correlation energy (29)relative to its exact value are shown in logarithmicscale in Fig. (2) for a slightly larger system of L = 20levels at half-filling. Note that the vCCD sep PBCS er-rors have slighly increased with respect to the L = 12case, but are still respectably accurate when comparedto PBCS. The limitations of the vCCD sep PBCS wave- G cr L =
20, M = PBCS vCCD ( w = ) PBCSvCCD sep
PBCS vCCD sep ( w = ) PBCSvCCD sep ( w = ) PBCS - - G / ϵ - E c / E c , e x ac t FIG. 2. Error in the correlation energy (29) relative to itsexact value (in percentages) versus the pairing strength G (in units of the level spacing (cid:15) ) for the fully varia-tional PBCS (6), vCCD sep PBCS (13), vCCD ( w =6) PBCS(15) and vCCD ( w =4 , PBCS (17) wavefunctions in the caseof an L = 20 level system at half filling. function in the weak pairing regime originate in theapproximate form of its structure matrix, as remarkedaround Eq. (14).This weak pairing behaviour may be improved ata reasonable computational cost by including ad-ditional non-separable CCD excitations on top ofvCCD sep PBCS. The vCCD ( w )sep PBCS ansatz (17) al-ready provides more than an order of magnitude lowererrors with respect to vCCD sep
PBCS at weak pairingeven for the modest w = 4. A larger non-separableCCD excitation window w = 6 naturally accounts foran additional amount of correlations, further reduc-ing the errors at weak pairing. Beyond the criticalvalue G cr however there is no significant benefit ofthe supplementary non-separable excitations. This isalso the situation for the vCCD ( w ) PBCS ansatz (15)which displays the same strong pairing behaviour asthe plain PBCS. With the Fermi sea being washed outat strong pairing, only a global deformation of PBCSsuch as vCCD sep
PBCS is able to bring substantial im-provement to the energetics, as opposed to any localdeformation such as vCCD ( w ) PBCS (15).The numerical values for the optimal non-separablestructure matrix elements z ph are typically found tobe very small, of the order 10 − − − upon thefull vCCD ( w )sep PBCS energy minimization. In practice,within the vCCD ( w )sep PBCS approach we chose to limit the non-separable CCD excitations to linear ordervCCD ( w ) ( z ) = 1 + M + w/ (cid:88) p = M +1 M (cid:88) h = M +1 − w/ z ph P † p P † h , (30)which allows for a decrease in the computational com-plexity without any noticeable loss in precision withrespect with the full form of Eq. (15).An additional computational speed-up is enabledby limiting the action of the particle number projec-tion operations in Eqs. (18) and (20) to a reducedsubspace of particle and hole collective pairs. For themoderate values of the pairing strength G consideredhere (2 to 3 times the value of G cr ) this still allowsfor an exact particle-number conservation (within thenumerical accuracy). This key computational aspectis ensured by the negligible contributions of the highorder terms within the PBCS (6) and vCCD sep (19)expansions, due to the small (subunitary) numericalvalues of their corresponding collective pair ampli-tudes. All results presented below for L = 60 and L = 100 were obtained upon projecting the particle-number only within a C = 15 ph-pair subspace forPBCS in Eq. (18) and within a C = 10 ph-pairsubspace for CCD sep in Eq. (20). The gauge angleswere adjusted accordingly, i.e. θ n = 2 πn/ ( C + 1), φ k = ϕ k = 2 πk/ ( C +1). We also note that for all con-sidered systems it was sufficient to restrict the amountof CCD sep excitations by imposing an upper bound at (cid:96) max = 7 in the (cid:96) -sum of Eq. (23), as to avoid thenumerical errors originating from the combination oflarge factorials and rapidly oscillating phases and, atthe same time, but preserving the numerical precision.We present in Fig. (3) the relative errors in thecorrelation energy (29) relative to its exact values for L = 60 (left panel) and L = 100 (right panel) athalf filling. The vCCD sep PBCS wavefunction remainssurprisingly accurate given that its doubly-separablestructure matrix (14) at weak coupling involves only L independent parameters (down from a total of 2 L parameters due to the present ph-symmetry) as com-pared to the fully non-separable CCD structure ma-trix that requires L / sep PBCS errors at weak coupling are0.3% and 0.5% for L = 60 and L = 100 respectively,increasing slightly until G ∼ G cr and then rapidly de-creasing at stronger couplings.While the inclusion of non-separable local CCD ex-citations within the vCCD ( w )sep PBCS ansatz (17) sig-nificantly improves the error at weak coupling evenfor relatively very small values of w , perfectly ac-curate energetics would still require a set of globalph-excitations on top of vCCD sep PBCS. Nevertheless,the overall quality of the vCCD ( w =6)sep PBCS results isat the level of the more involved particle-number pro-jected Bogoliubov-coupled-cluster theory of Ref. [18],while consistently providing an energy upper bound G cr L =
60, M = vCCD sep PBCS vCCD sep ( w = ) PBCSvCCD sep ( w = ) PBCS PBCS - CCD - G / ϵ - E c / E c , e x ac t ( % ) G cr L = = - G / ϵ - E c / E c , e x ac t ( % ) FIG. 3. Error in the correlation energy (29) relative to its exact value (in percentages) versus the pairing strength G (inunits of the level spacing (cid:15) ) for the fully variational vCCD sep PBCS (13) and vCCD ( w )sep PBCS (17) wavefunctions, in thecase of an L = 60 level system (left panel) and for an L = 100 level system (right panel) at half filling. The red squaresindicate the results of the particle-number projected Bogoliubov-coupled-cluster theory of Ref. [18]. Note also that for L = 100, the errors for the pure PBCS (6) reach a maximum of about 22% at G/(cid:15) = 0 .
24 and decrease down to 6% atvery small G (see also Fig. 1 of Ref. [18]). across all regimes.The situation is similar away from half-filling, asshown in Fig. (4). To provide a comparison with theresults of Ref. [18], we consider the L = 100 case atan interaction strength of G = 1 . G cr . For this valuethe particle-number symmetry is broken for all fillingfractions. Both PBCS and vCCD sep PBCS are mostaccurate for small (or large) filling fractions, but theirerrors live on different scales (0.5% vs 20% at half fill-ing). As seen in Fig. (3), in the considered regimethe non-separable local CCD excitations bring no im-provement over the vCCD sep
PBCS results, which arewell matched against those of Ref. [18].Finally, we consider the behaviour of the canonicalgap ∆ = G L (cid:88) i =1 (cid:112) n i (1 − n i ) , (31)where n i = (cid:104) c † i c i + c † ¯ i c ¯ i (cid:105) / i . This quantity exhibits amore pronounced sensitivity to the structure of thewavefunction than the correlation energy due to itsdependence on the occupation probabilities.We show in Fig. (5) the error for the canonical gap∆ relative to its exact value for L = 60 at half filling.Improving on both individual PBCS and vCCD sep wavefunctions, their combination vCCD sep PBCS onlyshows visible ∼
1% errors in the
G < G cr region.These are further reduced within the vCCD ( w )sep PBCS approach which exhibits highly accurate occupationsacross all regimes.
IV. SUMMARY AND CONCLUDINGREMARKS
In this work, we considered a variational approachfor the ground state of finite paired systems basedon a combined coupled-cluster and particle-number-projected BCS wavefunction. We benchmarked ourresults against the exact solution for a picket-fencemodel involving a pure pairing force acting on a spaceof doubly degenerate, equally distanced levels. Theanalyzed systems systems range from small ( L = M =12) to relatively large ( L = M = 100).We confirmed within the variational context thatthe combination of symmetry-restored mean field the-ory and coupled-cluster theory leads to a wavefunc-tion that is significantly better than either of the twotaken separately, which is also the main conclusionof Ref. [18]. By incorporating pure four-body corre-lations, our vCCD ( w )sep PBCS (17) is able to reproducewell the physics at weak pairing while also offering or-ders of magnitude smaller errors relative to the stan-dard pair-only PBCS (6) across all other regimes. Ourresults match the high level of precision of the particle-number projected Bogoliubov-coupled-cluster theoryof Ref. [18], while consistently providing an energyupper bound. L = PBCSvCCD sep
PBCSPBCS - CCD
Filling Fraction - E c / E c , e x ac t ( % ) FIG. 4. Error in the correlation energy (29) relative to itsexact value (in percentages) versus the filling fraction
M/L of an L = 100 level system for the fully variational PBCS(6) and vCCD sep PBCS (13) wavefunctions. The pairingstrength is G = 1 . G cr , with G cr the critical strengthat half filling. The blue diamonds indicate the results ofthe particle-number projected Bogoliubov-coupled-clustertheory of Ref. [18]. Computational limitations include the restrictionof the pure quartet correlations to a relatively smallwindow around the Fermi level (discussed in detailin the main text) and also the need of a restrictedspace for particle-number projection. While the latterdoes not spoil the exact particle-number conservation(within the numerical accuracy) for moderate valuesof the pairing strength, alternative approaches needto be considered for the very strong pairing regime G (cid:29) G cr . One possibility would involve limiting theseparable CCD excitations to linear order within thevCCD sep PBCS approach. Computations with thissimpler wavefunction could then be performed effi-ciently without resorting to any approximations foran improvement over the already very good PBCS re-sults for this regime.The vCCD ( w )sep PBCS (17) wavefunction was shownto be quite flexible but it is still affected by structurallimitations leading to a visible maximum in the energyerror around G ∼ . G cr in all analyzed cases. At-tempts at mitigating these effects have included treat-ing as independent variational parameters the factori-als in Eq. (23) originating from the vCCD sep expan-sion (9), with only marginal benefits.Possible avenues to explore could involve the vari-ational treatment of beyond-ph CCD excitations ontop of the symmetry-restored mean-field state as in | Ψ (cid:105) = exp L (cid:88) i,j =1 z ij c † i c † ¯ i c ¯ j c j | PBCS (cid:105) , (32) which would act as a more natural choice (albeit morecomputationally challenging) for a regime lacking awell-defined Fermi sea. As a first step one couldenvision building a multi-separable beyond-ph CCDexcitation operator which could be optimized withthe very recently developed methods for constructinglinearly-independent PBCS states [19].Finally, we leave to future studies extensions of thetheory presented in this work that could incorporatethe effect of seniority breaking terms within a gener-alized variational wavefunction. G cr L =
60, M = PBCSvCCD sep vCCD sep
PBCSvCCD sep ( w = ) PBCS G / ϵ - Δ / Δ e x ac t ( % ) FIG. 5. Error in the canonical gap (31) relative to itsexact value (in percentages) versus the pairing strength G (in units of the level spacing (cid:15) ) for the fully varia-tional PBCS (6), vCCD sep (9), vCCD sep PBCS (13) andvCCD ( w =6)sep PBCS (17) wavefunctions in the case of an L = 60 level system at half filling. ACKNOWLEDGMENTS
This work was supported by a grant of the Roma-nian Ministry of Education and Research, CNCS -UEFISCDI, project number PN-III-P1-1.1-PD-2019-0346, within PNCDI III, and PN-19060101/2019-2022; and by the Spanish Ministerio de Ciencia eInnovaci´on, and the European regional developmentfund (FEDER) project N º PGC2018-094180-B-I00.We thank the authors of Ref. [18] for sharing withus their numerical results. [1] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Theoryof superconductivity, Phys. Rev. , 1175 (1957).[2] A. Bohr, B. R. Mottelson, and D. Pines, Possible anal-ogy between the excitation spectra of nuclei and thoseof the superconducting metallic state, Phys. Rev. ,936 (1958).[3] K. Dietrich, H. J. Mang, and J. H. Pradal, Conserva-tion of particle number in the nuclear pairing model,Phys. Rev. , B22 (1964).[4] J. A. Sheikh and P. Ring, Symmetry projected hartreefock bogoliubov equations, Nuclear Physics A , 71(2000).[5] G. E. Scuseria, C. A. Jimenez-Hoyos, T. M. Hender-son, K. Samanta, and J. K. Ellis, Projected quasipar-ticle theory for molecular electronic structure, TheJournal of Chemical Physics , 124108 (2011).[6] A. Khamoshi, T. M. Henderson, and G. E. Scuseria,Efficient evaluation of agp reduced density matrices,The Journal of Chemical Physics , 184103 (2019),https://doi.org/10.1063/1.5127850.[7] N. Sandulescu and G. F. Bertsch, Accuracy of bcs-based approximations for pairing in small fermi sys-tems, Phys. Rev. C , 064318 (2008).[8] J. Dukelsky, S. Pittel, and C. Esebbag, Structure ofthe number-projected bcs wave function, Phys. Rev.C , 034313 (2016).[9] F. Braun and J. von Delft, Fixed-n superconductivity:The crossover from the bulk to the few-electron limit,Phys. Rev. Lett. , 4712 (1998).[10] J. Dukelsky and G. Sierra, Density matrix renor-malization group study of ultrasmall superconductinggrains, Phys. Rev. Lett. , 172 (1999).[11] R. Richardson, A restricted class of exact eigenstatesof the pairing-force hamiltonian, Physics Letters ,277 (1963).[12] G. Sierra, J. Dukelsky, G. G. Dussel, J. von Delft, andF. Braun, Exact study of the effect of level statisticsin ultrasmall superconducting grains, Phys. Rev. B , R11890 (2000).[13] J. G. Hirsch, A. Mariano, J. Dukelsky, and P. Schuck,Fully self-consistent rpa description of the many levelpairing model, Annals of Physics , 187 (2002).[14] J. Dukelsky, G. Dussel, J. Hirsch, and P. Schuck,Comparison between exact and approximate treat-ments of the pairing interaction for finite fermi sys-tems, Nuclear Physics A , 63 (2003).[15] T. M. Henderson, G. E. Scuseria, J. Dukelsky, A. Sig-noracci, and T. Duguet, Quasiparticle coupled clus-ter theory for pairing interactions, Phys. Rev. C ,054305 (2014).[16] M. Degroote, T. M. Henderson, J. Zhao, J. Dukelsky,and G. E. Scuseria, Polynomial similarity transfor-mation theory: A smooth interpolation between cou-pled cluster doubles and projected bcs applied to thereduced bcs hamiltonian, Phys. Rev. B , 125124(2016).[17] A. Rubio-Garc´ıa, D. R. Alcoba, P. Capuzzi, andJ. Dukelsky, Benchmarking the variational reduceddensity matrix theory in the doubly occupied configu-ration interaction space with integrable pairing mod-els, Journal of Chemical Theory and Computation ,4183 (2018), pMID: 29906104. [18] Y. Qiu, T. M. Henderson, T. Duguet, and G. E. Scuse-ria, Particle-number projected bogoliubov-coupled-cluster theory: Application to the pairing hamilto-nian, Phys. Rev. C , 044301 (2019).[19] R. Dutta, G. P. Chen, T. M. Henderson, andG. E. Scuseria, Construction of linearly independentnon-orthogonal agp states (2021), arXiv:2101.08911[physics.chem-ph].[20] T. M. Henderson and G. E. Scuseria, Geminal-based configuration interaction, The Jour-nal of Chemical Physics , 051101 (2019),https://doi.org/10.1063/1.5116715.[21] T. M. Henderson and G. E. Scuseria, Correlatingthe antisymmetrized geminal power wave function,The Journal of Chemical Physics , 084111 (2020),https://doi.org/10.1063/5.0021144.[22] J. Dukelsky, G. Dussel, J. Hirsch, and P. Schuck,Comparison between exact and approximate treat-ments of the pairing interaction for finite fermi sys-tems, Nuclear Physics A , 63 (2003).[23] T. M. Henderson, I. W. Bulik, and G. E. Scuse-ria, Pair extended coupled cluster doubles, TheJournal of Chemical Physics , 214116 (2015),https://doi.org/10.1063/1.4921986.[24] J. Ripoche, D. Lacroix, D. Gambacurta, J.-P. Ebran,and T. Duguet, Combining symmetry breaking andrestoration with configuration interaction: A highlyaccurate many-body scheme applied to the pairinghamiltonian, Phys. Rev. C , 014326 (2017).[25] T. M. Henderson and G. E. Scuseria, Correlatingthe antisymmetrized geminal power wave function,The Journal of Chemical Physics , 084111 (2020),https://doi.org/10.1063/5.0021144.[26] A. Khamoshi, G. P. Chen, T. M. Henderson, andG. E. Scuseria, Exploring non-linear correlators onagp, (2021), arXiv:2012.02377 [cond-mat.str-el].[27] Y. Qiu, T. M. Henderson, J. Zhao, and G. E.Scuseria, Projected coupled cluster theory, TheJournal of Chemical Physics , 064111 (2017),https://doi.org/10.1063/1.4991020.[28] R. W. Richardson, Numerical study of the 8-32-particle eigenstates of the pairing hamiltonian, Phys.Rev. , 949 (1966).[29] H. Ui and G. Takeda, A Class of Simple Hamiltoni-ans with Degenerate Ground State. II A Model of Nu-clear Rotation: Spontaneous Breakdown of RotationSymmetry and Goldstone Theorem for Finite Dimen-sional System, Progress of Theoretical Physics ,176 (1983), https://academic.oup.com/ptp/article-pdf/70/1/176/5367313/70-1-176.pdf.[30] T. Koma and H. Tasaki, Symmetry breaking andfinite-size effects in quantum many-body systems, J.Stat. Phys , 745 (1994).[31] C. Yannouleas and U. Landman, Symmetry breakingand quantum correlations in finite systems: studiesof quantum dots and ultracold bose gases and relatednuclear and chemical methods, Reports on Progressin Physics , 2067 (2007).[32] T. Papenbrock and H. A. Weidenm¨uller, Effectivefield theory for finite systems with spontaneously bro-ken symmetry, Phys. Rev. C , 014334 (2014).[33] P.-O. L¨owdin, Quantum theory of many-particle sys- tems. iii. extension of the hartree-fock scheme toinclude degenerate systems and correlation effects,Phys. Rev. , 1509 (1955).[34] I. Mayer, The spin-projected extended hartree-fockmethod, Advances in Quantum Chemistry, , 189(1980).[35] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Self-consistent mean-field models for nuclear structure,Rev. Mod. Phys. , 121 (2003).[36] K. Schmid, On the use of general symmetry-projectedhartree-fock-bogoliubov configurations in variationalapproaches to the nuclear many-body problem,Progress in Particle and Nuclear Physics , 565(2004).[37] C. A. Jim´enez-Hoyos, T. M. Henderson, T. Tsuchi-mochi, and G. E. Scuseria, Projected hartree-fock the-ory, The Journal of Chemical Physics , 164109(2012), https://doi.org/10.1063/1.4705280.[38] P. Ring and P. Schuck,The Nuclear Many-Body Problem (Springer-Verlag,Berlin Heidelberg, 1980).[39] J.-P. Blaizot and G. Ripka,Quantum Theory of Finite Systems (The MITPress, Cambridge, MA, 1985).[40] A. J. Coleman, Structure of fermion density ma-trices. ii. antisymmetrized geminal powers, Jour-nal of Mathematical Physics , 1425 (1965),https://doi.org/10.1063/1.1704794.[41] J. Dukelsky and G. Sierra, Density matrix renor- malization group study of ultrasmall superconductinggrains, Phys. Rev. Lett. , 172 (1999).[42] M. Sambataro and N. Sandulescu, Treatment of thelike-particle pairing of quartets, J. Phys. G: Nucl.Part. Phys. , 055107 (2013).[43] T. Duguet, Symmetry broken and restored coupled-cluster theory: I. rotational symmetry and angularmomentum, J. Phys. G: Nucl. Part. Phys. , 025107(2015).[44] T. Duguet and A. Signoracci, Symmetry broken andrestored coupled-cluster theory: Ii. global gauge sym-metry and particle number, J. Phys. G: Nucl. Part.Phys. , 015103 (2017).[45] T. Tsuchimochi and S. Ten-no, Bridging single-and multireference domains for electron correla-tion: Spin-extended coupled electron pair ap-proximation, Journal of Chemical Theory andComputation , 1667 (2017), pMID: 28240896,https://doi.org/10.1021/acs.jctc.7b00073.[46] V. V. Baran and D. S. Delion, A quartet bcs-like the-ory, Physics Letters B , 135462 (2020).[47] A. Zee, Quantum field theory in a nutshell (Prince-ton University Press, 2010).[48] K. Peeters, Cadabra2: computer algebra for field the-ory revisited, Journal of Open Source Software3