Eigenvector continuation for the pairing Hamiltonian

The development of emulators for the evaluation of many-body observables has gained increasing attention over the last years. In particular the framework of eigenvector continuation (EC) has been identified as a powerful tool when the Hamiltonian admits for a parametric dependence. By training the emulator on a set of training data the many-body solution for arbitrary parameter values can be robustly predicted in many cases. Furthermore, it can be used to resum perturbative expansions that otherwise diverge. In this work, we apply EC to the pairing Hamiltonian and show that EC-resummed perturbation theory is in qualitative agreement with the exact solution and that EC-based emulators robustly predict the ground-state energy once the training data are chosen appropriately. In particular the phase transition from the normal to the superfluid regime is quantitatively predicted using a very low number of training points.


I. INTRODUCTION
In ab initio nuclear structure calculations new challenges emerged over the last years.The combination of basis-expansion methods with chiral effective field theory interactions has enabled routine computations of medium-mass nuclei [1,2].This has led to the quest to develop accurate interactions for medium-mass nuclei, e.g., by including medium-mass observables in the construction of new chiral interactions [3,4], and to explore uncertainty estimates from the effective field theory truncation and uncertainties in the low-energy couplings [5].However, both of these aspects require repeated solutions of the many-body problems when varying the underlying interactions.To cope with these significant computational demands, emulators have emerged as powerful tools to efficiently lower the computational burden of performing explicit calculations by mimicking the true many-body solution.
For this purpose reduced basis methods (RBM) and more general model order reduction (MOR) techniques have been widely applied in the natural sciences to reduce the computational cost by projecting the problem onto an appropriately chosen subspace [6][7][8][9][10][11].In nuclear theory, RBMs have been applied to density functional theory [12,13], to the solution of the Schrödinger equation for single-particle Hamiltonians [14] and to scattering problems.For an overview over different MOR methods see Refs.[15,16].Eventually, MOR methods also enable uncertainty quantifications.A particular variant of RBMs is the so-called eigenvector continuation (EC) approach.The EC framework has emerged as a particularly versatile method to emulate solutions of the many-body problem governed by parametric dependencies [17].Over the last years, the EC method has been applied to various few-and many-body problems, e.g., for uncertainty quantification in few-nucleon systems [18] and mediummass nuclei [19], for scattering and reactions [20][21][22][23], for finite-volume extrapolation [24] and nuclear matter computations [25].Moreover, EC has been used as a resummation tool for perturbative expansions enabling a robust extraction of many-body observables when the perturbation series diverges [26][27][28].
Practically, the emulator construction is based on the solution of a generalized eigenvalue problem on a small subspace spanned by a set of non-orthogonal many-body states defining the training vectors.From that training manifold, high-dimensional parameter spaces can be exhausted at a tractable computational cost.While in few-body applications the training vectors are virtually exact, the development for many-body calculations is accompanied by the introduction of approximation schemes that facilitate the emulator construction.While various applications from exact training vectors exist, the interplay of emulator quality and truncated training vectors has not been studied in detail.
This work is dedicated to a detailed study of the exactly solvable pairing Hamiltonian [29,30] that serves as a model for nuclear superfluidity and is a common test ground for novel many-body frameworks [31].We will first establish EC as a robust resummation tool for many-body perturbation theory (MBPT) even though the bare perturbative expansion breaks down.In the second part of this work, we investigate the accuracy of an EC-based emulator and its sensitivity to the selected training points1 used in the construction.Recently, the same system has been studied using RBMs built from approximate density matrix renormalization group (DMRG) wave functions [32].

II. EIGENVECTOR CONTINUATION
Eigenvector continuation provides a powerful framework for the construction of emulators whenever the Hamiltonian admits for a parametric dependence H(g 1 , ..., g N ) with all g ∈ R [17].By training the EC emulator on a small set of training points in the parameter space, it may robustly predict many-body observables for different parameter values without the need of performing explicit many-body simulations.In nuclear structure applications such a tool is particularly powerful, because the number of low-energy couplings in chiral effective field theory interactions may be systematically varied with an accurate emulator [19,25].
The construction of the EC emulator is based on the definition of a set of where |Ψ⟩ denotes the many-body state of interest, e.g., the ground state.Energies in the EC framework are evaluated at the target coupling g • ∈ R by solving the generalized eigenvalue problem where the Hamiltonian and norm kernel are obtained for the training vectors as Since the number of (non-orthogonal) training vectors {|Ψ(g i )⟩} is much smaller than the size of the configuration basis in diagonalization approaches, the EC diagonalization is of very limited cost.The solution of the generalized eigenvalue problem admits for N EC solutions, so that the EC approach provides also access to excited states belonging to the same symmetry class, e.g., spin or parity.This was recently studied for the case of the anharmonic oscillator [28].Finally, we note that the EC approach can be used for any other operator O instead of H, thus enabling calculations of other observables [18,19].

III. PAIRING HAMILTONIAN
In this work we study the pairing Hamiltonian where p denotes the time-reversed state of p, Ω denotes the number of two-fold degenerate levels with singleparticle energies ϵ p , and g the (real-valued) strength of the two-body interaction.The single-particle spectrum is taken to be equidistantly spaced ϵ p = p ∆ϵ, with a level spacing of ∆ϵ = 1 (in natural units).The pairing Hamiltonian has an exact Richardson solution, which is accessible for arbitrary coupling values and system size [29,30].The Richardson solution is obtained by solving a set of non-linear coupled equations for the pair energies E α , where α = 1, . . ., N occ with the number of occupied pair states N occ , Here, the first sum runs over all possible single-particle states, while the second sum is restricted to occupied pair states.In the following, we will study the system at half filling, i.e., half of the shells are doubly occupied yielding a pair number of N pair = Ω/2.The final ground-state energy is given by α E α .In our numerical benchmark we employ the formalism described in Ref. [33] to solve Eq. ( 5).
It is well known that the pairing Hamiltonian undergoes a transition from a normal to a superfluid phase once a critical (positive) coupling value g crit is exceeded.For g > g crit the system is no longer well approximated by low-rank particle-hole excitations but dominated by the formation of Cooper pairs resembling the superfluid character of the system.This transition can be qualitatively captured by using Bardeen-Cooper-Schrieffer (BCS) mean-field theory, where U (1) particle-number symmetry is explicitly broken at the mean-field level to account for the pairing correlations in the system [34,35].Therefore, the corresponding BCS vacuum |Φ BCS ⟩ is not an eigenstate of the particle-number operator A, but only restricted to have correct average particle number, i.e., N occ = ⟨Φ BCS |A|Φ BCS ⟩.The emergence of pairing correlations is further accompanied by a breakdown of standard symmetry-conserving correlation expansions and requires the use of quasiparticle reformulations to account for the collective effects induced by the static correlations [36][37][38].
In our benchmarks we employ Rayleigh-Schrödinger perturbation theory using the partitioning where H 0 is taken to be the (diagonal) one-body part of the pairing Hamiltonian.The ground-state wave function is expanded in the series where |Φ (p) ⟩ denotes the p-th order state correction [39].
Corrections to the ground-state energy are calculated via ⟩.The reference state |Φ (0) ⟩ is given by the Hartree-Fock (HF) state |Φ HF ⟩ [40,41].Moreover, we employ the n-particle-n-hole pair configuration interaction (pCI-npnh) framework, where the A-body configuration space is spanned by all pair excitations up to n-particles-n-holes from of the HF state |Φ HF ⟩.The low-lying spectrum is obtained from diagonalizing the many-body Hamiltonian using Lanczos techniques.While this does not provide a scalable solution for large systems, it is numerically tractable for Ω ≲ 20.In the case the maximum number of n-particle-n-hole excitations is given by n = 2N pair , the exact result (full CI, FCI) is recovered.

IV. EIGENVECTOR CONTINUATION AS A RESUMMATION TOOL
Eigenvector continuation can be used to resum the perturbation series by performing the EC in the space of MBPT state corrections [26,27].The Hamiltonian and norm matrices that enter Eq. ( 2) are then given by where N p0 = δ p0 by intermediate normalization.In this work, MBPT state corrections up to first order are included; therefore we denote this approach as EC-PT(1).Explicit expressions for the matrix elements in Eqs.(8a) and (8b) are evaluated by virtue of Wick's theorem [42].
The simplest expression involves the second-order energy correction where ϵ a i = f i − f a is a shorthand notation for the energy denominators (with i, j referring to occupied hole states and a, b to particle states) and f p = ϵ p − gn p denotes the normal-ordered one-body part.The overlap is given by The most complicated expression corresponds to the matrix element between two first-order state corrections and consists of three contributions where the superscript indicates from which 0-, 1-, or 2body part of the Hamiltonian the contribution originates.For the next order EC-PT(2), second-order state corrections that contain 4p4h-excited components are needed.
Figure 1 shows a comparison of different many-body methods for the solution of the pairing Hamiltonian.In the normal phase (g < g crit ) the HF and BCS solution coincide, since the coupling strength is too weak to support a superfluid ground state [35,43].For g > g crit the BCS state explores the enlarged variational space, yielding a lower energy than the symmetry-conserving HF solution.Both HF/BCS mean-field approaches lack all dynamical particle-hole correlations leading to the deviations from the exact solution with increasing coupling.
As shown by the second-and third-order MBPT results, the perturbative expansion is only reliable for moderate couplings and fails for larger values.However, the EC-PT(1) results are in good agreement with the exact solution for a wide range of couplings.The EC framework thus effectively resums particle-hole correlations emerging from leading order perturbative corrections to the ground state.In addition, EC-PT(1) and pCI-2p2h give very similar results for all coupling values.This is to be expected since both frameworks probe up to two-particle-two-hole correlations in the configuration space, either by a large-scale expansion in Slater determinants in pCI-2p2h or by a small-scale expansion of multi-reference states arising from the perturbative state amplitudes in EC-PT (1).Here we emphasize that in this work the evaluation of the EC kernel is performed based on a diagrammatic approach and is hence only of polynomially scaling complexity.This makes it possible to compute EC-PT(1) corrections for arbitrary system sizes.In contrast, in Ref. [26] the perturbative corrections at different orders were computed recursively using an exponentially large configuration basis, which restricts this approach to limited basis spaces Ω ≤ 20.
As shown in Fig. 2, the EC-PT(1) results perform at a similar level as state-of-the-art nonperturbative coupledcluster doubles (CCD) [44,45] and in-medium similarity renormalization group (IMSRG) [46,47] calculations when restricted to the normal phase, which is relevant for ab initio calculations of medium-mass nuclei.Still the CCD and IMSRG results are consistently closer to the exact solution which we attribute to a more elaborate resummation of higher-order effects that are absent from the simpler EC-PT(1) estimate.Once further corrections are included, e.g., using EC-PT(n ≥ 2), the variational character of the EC ensures improvement of our results towards the Richardson solution.
The failure of conventional basis-expansion approaches can be understood in terms of the norm matrix shown in Fig. 3 obtained by sampling n = 200 couplings from the interval g ∈ [−2, 2] and evaluating the overlap of the resulting ground-state eigenvectors.Matrix elements |N pq | ≈ 1 indicate similar structures in the ground states, |Ψ(g p )⟩ and |Ψ(g q )⟩.In the case where g p and g q belong to the same regime (either normal or superfluid) the eigenvectors show strong linear dependencies, thus, indicating similar correlations for the two couplings.However, for g p < g crit < g q the overlap is significantly reduced yielding almost orthogonal eigenvectors from the two regimes.As a consequence, one cannot expect a Slater-determinant-based correlation expansion to give accurate results in the superfluid regime [48].This also holds for non-perturbative expansions since the breakdown is related to the use of an improper reference state and not a lack of non-perturbative ladder-type resummations of particle-hole correlations.Figure 3 also nicely illustrates the weak dependence of the eigenstates on the parameters in each regime, which is fundamental to the EC.

V. EIGENVECTOR CONTINUATION AS AN EMULATION TOOL
Next, we use EC to construct a many-body emulator from a selected set of training data.Since the groundstate dynamics in the normal and superfluid regimes are fundamentally different, we will investigate the performance of the EC approach as a function of the training data that enter the definition of the EC manifold in Eq. (1).Therefore, three different scenarios characterized by different training sets will be investigated: Naturally, this gives rise to three emulator types: EC(D < train ), EC(D > train ), and EC(D mixed train ), that have been trained using the corresponding regimes.
Figure 4 shows the EC emulator predictions for the ground-state energy for the three different training sets as well as the deviation from the exact result.By construction, the error at the training points vanishes because the exact many-body state was used as training vector.We observe that the EC(D < train ) emulator (dashed line) performs well in the normal regime while incorrectly predicting the superfluid regime.Once the critical coupling is approached the error rapidly increases and no quantitative prediction is possible.Similarly, EC(D > train ) (dotted  line) accurately predicts the superfluid regime while being incapable of capturing the correct trend in the normal regime.These limitations are overcome by employing a mixed emulator EC(D mixed train ): now the normal and superfluid regimes are correctly described and only small deviations appear in between training points.
We conclude that the inclusion of both dynamical particle-hole correlations and static pairing correlations are crucial in the training stage since otherwise the EC manifold probes an insufficient subspace for capturing the Hamiltonian's dynamics at arbitrary coupling values.These findings are consistent with the form of the norm matrix (see Fig 3): while in the same regime two vectors have a large overlap, there is nearly no overlap between the normal and the superfluid regime.Therefore, it is not surprising that only selections of training points from both regimes can approximate the solution across both ranges of couplings.Finally, we compare the EC as an emulator to the EC-PT(1) prediction.Here large errors have to be expected, since only 2p2h-excitations span the EC-PT(1) manifold, lacking important higherbody excitations in the normal phase.The error can be systematically improved by including higher-order state corrections in EC-PT(n).However, the breakdown beyond the critical coupling is associated to the insufficient HF reference state in the presence of superfluid pairing correlations.
Finally, we comment on the impact of enlarging the training set.While for two training points it was not possible to approximate the solution qualitatively over the considered coupling range, three training points were sufficient to approximate the solution.Once a mixed training set is chosen, a further increase in the number of training vectors does not significantly increase the quality of the emulator.This is consistent with the observation from the norm kernel, which reveals strong linear dependencies among vectors from the same regime.Therefore, adding (almost) co-linear vectors will not further improve the emulated results.The location of training points seems to be more important than the total number of training points.Once the various physical regimes are covered in the training stage, the emulator quantitatively predicts the exact solution very well.Our general findings agree with the results from Ref. [32] that employ an iterative update of the RBM manifold.

VI. TRUNCATION OF THE TRAINING VECTORS
The availability of the exact many-body solution for a given parameter set is almost never fulfilled in realistic applications and hence the design of emulators naturally depends on the many-body approximations.Therefore, the performance of EC emulators has to be validated when approximate training vectors are employed.To this end, we employ the n-particle-n-hole pair configuration interaction (pCI-npnh) to obtain training data at the 2p2h and 4p4h truncation.In general, the quality of truncated EC emulators is limited by the quality of the many-body approximations used for the training vectors.Therefore, improvements of the truncated EC results against the exact many-body solution need to be considered with care.
In Fig. 5 we show the EC ground-state energies for 2p2h-truncated and 4p4h-truncated training data and for the exact (10p10h or FCI) training data for the case of the mixed training set.The training-point truncations in between are not shown, as they differ little from the FCI case.In the top panel of Fig. 5, we observe that the difference between the EC and pCI-npnh energies increases as the quality of the many-body approximation improves.Comparing truncated EC and pCI-npnh results against the exact results shows at small couplings a probably accidental improvement at lower particle-hole truncations, but overall the better the many-body approximation the better is the truncated EC and pCI-npnh.This is especially pronounced in the superfluid regime.To explore this further, we consider in Fig. 6 the superfluid training set.Also for this training set, the difference between the EC and pCI-npnh energies is smallest at the 2p2h level.Moreover, for this case of the superfluid training set, the EC at the 4p4h and FCI level does not reproduce the CI (or the exact) solution well when extrapolated to the normal regime.
As can be seen both in the lower panels of Fig. 5 and Fig. 6, EC better emulates lower CI truncations.This is not surprising, since the size for the EC subspace is three for all three training point truncations.However, while for pCI-2p2h the three-dimensional subspace is used to approximate a 26-dimensional subspace, for full CI it tries to emulate the full 252-dimensional space, which is much larger.This makes emulating full CI a more complex problem.Therefore, while the more accurate full CI approximation gives the better approximation close to the training points, further away from the training points, it is more difficult for EC to emulate the larger model space.

VII. SUMMARY AND CONCLUSIONS
In this Letter we applied the EC framework to the pairing Hamiltonian.In the first part we used EC to resum the perturbative expansion based on a HF state.In the second part, we designed an EC emulator and tested its sensitivity on the location of training data and manybody approximations for the training vectors.The ECresummed perturbative results were in qualitative agreement with the exact Richardson solution even though the underlying perturbative corrections indicated a divergent perturbative expansion.The employed EC-PT(1)truncation gave similar results as the CI results truncated at 2p2h-level.We expect that our EC-PT(1) results can be efficiently improved by designing a quasiparticle extension built on a BCS reference state [41,48].
For the emulator design, the proper choice of training data was crucial.If the training set was entirely located in the normal (superfluid) regime the predictions for the superfluid (normal) regime gave rise to large errors.By incorporating training points from both regimes (even just one from the respective other regime) a quantitative prediction of the phase transition was achieved and only moderate errors were encountered.While these errors can be systematically reduced by including training points close to the critical value, the specific uncertainties will generally depend on details like level spacing and system size.In addition, we explore the CI truncation of the training vectors.This showed that, for a good choice of training points, the higher the truncation of the training vectors the better the EC approximation.Surprisingly, further away from the training data, there is a tendency that lower truncated training points give better results.This could be traced to the EC working better for lower truncated CI in these cases.
Future studies will focus on the power of EC emulators for performing large-scale sensitivity studies with respect to variations of the low-energy couplings in chiral Hamiltonians.While pioneering work has already been performed [19,25], more work is needed to develop emulators for other ab initio methods and to explore the quality of the emulator in the presence of many-body approximations.Ultimately, EC provides a powerful framework for global statistical analyses that otherwise remain intractable when exhausting high-dimensional parameter spaces through explicit many-body calculations from first principles.

Figure 1 .Figure 2 .
Figure 1.Comparison of different basis-expansion methods (for details see text) for the ground-state energy of the pairing Hamiltonian as a function of the coupling g for Ω = 10 at halffilling Nocc = 5.The critical coupling gcrit = 0.34 is indicated by the vertical line.
i) a normal training set D < train , where all training points are in the normal regime, ii) a superfluid training set D > train , where all training points are in the superfluid regime, iii) a mixed set D mixed train , where training points from both regimes are present.

Figure 4 .
Figure 4. EC results for the ground-state (top panel) of the pairing Hamiltonian for three different training sets (dashed, dash-dotted, and dotted lines with three training points indicated in the legend and by the crosses) in comparison to the exact and EC-PT(1) results.The bottom panel shows the deviation from the exact result.Results are given for Ω = 10 at half-filling Nocc = 5, and gcrit = 0.34 is indicated by the vertical line.

Figure 5 .
Figure 5. Difference between the EC and pCI-npnh groundstate energies (dashed lines in the lower panel) for different truncations (top panel) of the pairing Hamiltonian as well as their deviation from the exact result (bottom panel).The EC is based on the case of the mixed training points gi ∈ {−0.9, 0.3, 0.9}.Results are given for Ω = 10 at half-filling Nocc = 5.