Improved many-body expansions from eigenvector continuation

Quantum many-body theory has witnessed tremendous progress in various fields, ranging from atomic and solid-state physics to quantum chemistry and nuclear structure. Due to the inherent computational burden linked to the ab initio treatment of microscopic fermionic systems, it is desirable to obtain accurate results through low-order perturbation theory. In atomic nuclei however, effects such as strong short-range repulsion between nucleons can spoil the convergence of the expansion and make the reliability of perturbation theory unclear. Mathematicians have devised an extensive machinery to overcome the problem of divergent expansions by making use of so-called resummation methods. In large-scale many-body applications such schemes are often of limited use since no a priori analytical knowledge of the expansion is available. We present here eigenvector continuation as an alternative resummation tool that is both efficient and reliable because it is based on robust and simple mathematical principles.

(Received 29 November 2019; accepted 12 March 2020; published 9 April 2020) Quantum many-body theory has witnessed tremendous progress in various fields, ranging from atomic and solid-state physics to quantum chemistry and nuclear structure. Due to the inherent computational burden linked to the ab initio treatment of microscopic fermionic systems, it is desirable to obtain accurate results through low-order perturbation theory. In atomic nuclei, however, effects such as strong short-range repulsion between nucleons can spoil the convergence of the expansion and make the reliability of perturbation theory unclear. Mathematicians have devised an extensive machinery to overcome the problem of divergent expansions by making use of so-called resummation methods. In large-scale many-body applications, such schemes are often of limited use since no a priori analytical knowledge of the expansion is available. We present here eigenvector continuation as an alternative resummation tool that is both efficient and reliable because it is based on robust and simple mathematical principles. Introduction. The quantum-mechanical treatment of many interacting particles from first principles poses a formidable formal and computational challenge. At the heart of the ab initio philosophy is the capacity to control the error induced by many-body approximations used to solve the stationary Schrödinger equation. Light systems, e.g., atomic nuclei with small mass number A, can be accurately described by largescale diagonalization approaches, such as configuration interaction (CI) [1][2][3] or by imaginary-time propagation methods via Quantum Monte Carlo techniques [4][5][6]. However, the exponential growth of the Hilbert-space dimension with increasing A constitutes an ordeal for many-body practitioners interested in systems containing tens or hundreds of particles. In this Rapid Communication, we present a novel approach to the many-body problem that uses the method of eigenvector continuation to accelerate the convergence of calculations using many-body perturbation theory. As we show, the method can produce accurate results even in cases where standard approaches fail to converge.
In fact, there are two characteristics of the two-nucleon interaction making the many-body problem hard to solve, i.e., a priori nonperturbative [8,[37][38][39]. The first one relates to strong high-momentum correlations induced by the short-range repulsive interactions noted above. Fortunately, this problem can be tamed via the use of renormalization-group transformations. The second issue stems from the large scattering lengths corresponding with resonant scattering in the S-wave channels, which can produce strong many-body correlations. Although Pauli-blocking and effective-range effects ameliorate this problem to some extent, one of the motivations of approaches, such as Bogoliubov many-body perturbation theory (BMBPT), is that the use of symmetry-breaking reference states allows one to include more long-range correlations into the reference state, which can help to tame this problem further in nuclei.
Although we have reduced some of the problems of strong correlations, we may find that the perturbative expansion is still not convergent and requires using resummation tools. However, since the analytical properties of the many-body expansion are usually unknown, conventional resummation methods often cannot be applied with confidence. Therefore, an alternative framework that does not require such knowledge is highly desirable. In this Rapid Communication, we present eigenvector continuation (EC) as a technique that allows to achieve this goal. Two numerical applications dedicated to the closed-and open-shell nuclei 3 H and 18 O, respectively, are provided to illustrate the power of the method.
Eigenvector continuation. Recently, the EC technique was introduced [40,41] to treat physical systems the Hamiltonian of which is a particular instance, e.g., H ≡ H (1) of an operator H (c) depending smoothly on a parameter c. The power of EC relies on (1) the fact that there exists a regime 0 c c e < 1 for which the many-body problem is easier to solve than for the target value c = 1, (2) the stability of the eigenvectors of H (c) against variations of c, i.e., the targeted many-body state remains in a low-dimensional manifold of the full A-body Hilbert space when changing c from [0, c e ] to 1.
Based on these principles and targeting a particular eigenstate of H, e.g., the nuclear ground state, EC works in two successive steps, In practice, it is typically advantageous to first diagonalize the norm matrix and eliminate eigenvectors associated with its very small eigenvalues that arise when choosing the c i from a narrow range.
Perturbation theory. The EC method is, in principle, agnostic with respect to the particular computational method used to generate the N EC auxiliary states {| (c i ) }. In practice, the performance of EC does, of course, depend on the ability to describe the ground state of H (c i ) with sufficient accuracy. In this Rapid Communication, the low-dimensional EC manifold is built from perturbative corrections on top of a well-chosen reference state. Although the perturbative expansion is, in general, not suited to reach directly H (1), EC can be understood as effectively performing a sequence of analytic continuations and is, thus, able to go beyond the radius of convergence of the perturbative expansion whereas using the same inputs. It has been observed empirically that EC even works in cases where the radius of convergence of perturbation theory is zero and we use a variational subspace spanned by derivatives of the eigenvector at zero coupling [42]. The only requirement is that the dependence of H (c) on c is sufficiently smooth.
The procedure starts from the partitioning of the Hamiltonian according to H ≡ H 0 + H 1 such that the reference state is the ground state of the unperturbed Hamiltonian, The eigenstates of H 0 obtained through elementary, e.g., particle-hole, excitations of | (0) provide an orthonormal basis of the many-body Hilbert space. Scaling the residual interaction H 1 by a parameter c to introduce the parameterdependent Hamiltonian H (c) ≡ H 0 + c H 1 , perturbation theory (PT) generically parametrizes the exact ground state of the latter via an infinite power series, characterized by an (unknown) radius of convergence c ∈ [0, R c ]. In Eq. (2), | (p) , p 1 denotes the perturbative state correction of order p, which is independent of c and typically given as a specific linear combination of the eigenstates of H 0 [43].
Choosing 0 c i c e R c , i = 1, . . . , N EC and working at PT order P, the approximate low-dimensional manifold can be related to the reference state and the first P state corrections via the transformation, ⎛ . . .
such that the EC eigenvalue problem is equivalently formulated within the manifold {| (p) ; p = 1, . . . , P}. Eventually, the actual values {c i ; i = 1, . . . , N EC } do not matter, and the dimensionality of the manifold is effectively set by the perturbative order P. One, thus, needs to compute the (P + 1)(P + 1) matrices, where the diagonal matrix E gathers P + 1 eigenenergies, the lowest of which relates to the ground state. 1 Eventually, the procedure reduces to calculating PT state corrections at a chosen order P, computing the (P + 1) × (P + 1) Hamiltonian and norm matrices as well as solving the associated secular equation. Last but not least, the eigenvectors can further be used to compute other observables of interest. While the present Rapid Communication provides proofof-concept results up to high order P, future work will target fully realistic calculations at low orders, e.g., P = 3, that are currently within reach of state-of-the-art numerical codes implemented in large model spaces.
Applications. We now consider the MBPT-based EC scheme applied to two situations of increasing complexity with the goal of illustrating the difficulties posed by the two sources of nonperturbative behavior mentioned at the 1 Although the ground state is targeted, Eq. (4) actually delivers P + 1 vectors and energies. It will be interesting to investigate to what extent the excited states within the EC subspace can be associated with states in the spectrum of H carrying the same quantum numbers as the ground state. In fact, the excited states of a Hermitian operator in the EC subspace are variational with respect to the corresponding large-space spectrum [44]. outset. In both cases, the calculations are performed up to high perturbative orders based on a recursive scheme, and the interaction employed is the EM500 interaction [45] with three-nucleon forces omitted for simplicity. Furthermore, a similarity renormalization group (SRG) transformation characterized by a continuous flow parameter λ is applied to soften the nuclear Hamiltonian [46].
The first application is dedicated to the very light nucleus 3 H. The exact result and the recursive scheme to perform MBPT are based on a CI code built within the harmonicoscillator Jacobi-coordinate no-core shell-model formalism [47]. The large three-body Hilbert (sub)space employed includes configurations up to N max = 12 excitations above the harmonic-oscillator Slater-determinant reference state. In Fig. 1, the ground-state energy is displayed as a function of the MBPT order P for the unevolved EM500 interactions and for the SRG-evolved potential with λ = 2.0 fm −1 . Results from the CI diagonalization are compared to MBPT as well as to EC and (diagonal and superdiagonal) Padé resummations built on top of it [48].
When employing a "hard" interaction displaying large lowto-high momentum coupling (left panels), the perturbative series exhibits a divergent behavior. Using EC on top of it yields a rapid and monotonic convergence towards the exact CI result. The EC and Padé results are very similar. The effect of "softening" the short-range behavior of the interaction via the SRG transformation is seen in the right panels: The perturbation series oscillates much more mildly and seems to 041302-3 converge towards the CI limit, although it, in fact, diverges at even higher orders. Even in this case, EC improves upon MBPT by exhibiting a monotonic and rapid convergence providing, at each order, a variational upper bound to the CI result. Again, EC and Padé results are very similar. At P = 4, their deviation from the CI result is already below 5%. The same behavior is found for 4 He and for various SRG parameters.
Eigenvector continuation based on the first few MBPT orders manages to reveal the exact result even for the unevolved Hamiltonian where the series displays a divergence of ultraviolet character. It must be noted that, although the above calculations have been performed using a reference Slater determinant built from the spherical harmonic-oscillator onebody basis, the divergence of the MBPT series can actually be overcome by starting from an optimized (e.g., Hartree-Fock) mean-field reference state [13]. In open-shell nuclei, the infrared divergence associated with the emergence of nuclear superfluidity and/or deformation requires an even more drastic optimization of the reference state, i.e., of H 0 . This can be achieved by using symmetry-breaking references states, e.g., using a Bogoliubov extension of MBPT as we will discuss in the following.
Neither EC nor any other resummation technique is useful when the reference state is degenerate with respect to elementary excitations since this prevents the construction of the first few perturbative orders. A powerful way to resolve the issue is to allow H 0 and | (0) to break symmetries of H. For example, nuclear superfluidity can be handled at the price of breaking U (1) symmetry associated with the conservation of particle number. In this case, the reference state becomes a Bogoliubov vacuum that already captures (most of) the so-called static correlations originally responsible for the infrared divergence. Building perturbation theory on top of such a reference state led recently to the introduction of BMBPT [49]. In BMBPT, the average particle number must be monitored and adjusted to the targeted value A at each perturbative order. Thus, BMBPT qualifies as a perturbation theory under constraint [50]. Recent large-scale second-and third-order BMBPT calculations of midmass open-shell nuclei demonstrated the merits of the method and provided ground-state energies in agreement with the most sophisticated nonperturbative many-body schemes on a few percents level at a small fraction of the computational cost [16].
The second application is, thus, dedicated to the open-shell nucleus 18 O and employs the soft interaction characterized by the SRG parameter λ = 2.0 fm −1 . Compared to 3 H, the CI result and the recursive scheme to perform BMBPT require a more drastic limitation of the model space. First, a one-body harmonic-oscillator basis utilized in the present application is severely truncated 2 at e max = 2n + l = 4. A symmetry-broken Hartree-Fock-Bogoliubov (HFB) reference state is obtained in that model space by solving the variational mean-field HFB equations [51]. Subsequently, the many-body basis of H 0 used to expand perturbative state corrections is 2 A converged ab initio calculation with respect to the size of the one-body basis would typically require e max = 2n + l = 13 [16]. limited to configurations obtained via two-, four-, and selected [52] six-quasiparticle excitations of the reference state. 3 The latter truncation leads to performing approximate BMBPT calculations at order P 3 where basis states associated with eight quasiparticles and beyond do contribute to state corrections. The CI result acting as a benchmark is obtained in the same configuration basis. Figure 2 displays the ground-state energy of 18 O as a function of the perturbative order P. Results are shown for BMBPT as well as for BMBPT-based (diagonal) Padé resummation and EC. One first observes that variationally optimizing the U (1)-breaking Bogoliubov reference state does not prevent the perturbative series to diverge. Indeed, best coping with the most dramatic, i.e., infrared, divergence is achieved at the price of inducing an ultraviolet one [50] that was not present in 3 He for λ = 2.0 fm −1 . Still, the tremendous benefit of the approach is to make low MBPT orders meaningful in open-shell nuclei such that they actually provide a decent account of the exact result. Performing a parametric Padé resummation slightly reduces the oscillations for intermediate orders, but the resummed expansion still diverges. Only with EC, a rapid and monotonically converging series towards the CI result is obtained.
As a matter of fact, EC based on a low-order BMBPT calculation provides an accurate upper bound of the CI result with a subpercent error for P 1. The method, thus, constitutes a very promising tool for fully realistic 4 calculations where low-order BMBPT can be applied at moderate computational cost. One must, of course, beware that, by continuing approximate eigenvectors | P (c) obtained at small values of c, EC can, at best, account for the CI result at c = 1 within the eigensubspace of H 0 covered by the approximate eigenvectors. Realistic BMBPT(2) calculations probe, at most, zero-, two, four-, and (selected) six-(and, in principle, eight-) quasiparticle excitations. Thus, the corresponding EC result can capture nonperturbative correlations associated with up to six-quasiparticle (three-particle and three-hole) excitations, which is superior to what available state-of-the-art nonperturbative methods can deliver today for midmass (open-shell) nuclei.
Conclusions and outlook. In this Rapid Communication, a novel nonparametric resummation method was applied to perturbative expansions appropriate to closed-and open-shell nuclei. Calculations were based on a realistic nuclear twobody interaction derived from chiral effective field theory. Although the perturbative expansion is divergent in most cases, EC provides a robust framework to obtain a rapidly converging resummed sequence.
In particular, EC elegantly overcomes divergences arising from possible strong high-momentum couplings in nuclear interactions and/or induced by the interference between them and strong infrared correlations as long as the perturbation theory exploits a symmetry-breaking reference state. Benchmark calculations reveal that EC based on low-order perturbation theory is vastly superior to parametric resummation schemes, such as Padé approximants that are widespread among many-body practitioners. Consequently, EC provides an efficient and accurate tool to improve many-body expansions independently of the origin of the Hamiltonian, thus, being very promising for various fields of many-body physics and chemistry.
While the results presented in this Rapid Communication are of proof-of-principle character, EC will soon be applied to fully realistic ab initio large-scale calculations of midmass closed-and open-shell nuclei. This will be performed on top of low-order perturbative calculations performed in a large single-particle basis set and starting from realistic two-and three-nucleon interactions. Since the EC method yields a variational upper bound for binding energies, it also provides a first step towards the quantification of many-body uncertainties from nonvariational many-body techniques. Notably, EC also gives access to other groundstate observables [41] and low-lying excitations at the same time. Furthermore, the method can be applied to infinite nuclear matter calculations.
Although the present Rapid Communication implements EC on top of perturbation theory, it can be similarly applied on top of more sophisticated nonperturbative expansion methods, such as coupled cluster theory [53]. It will be interesting to see if a further optimization of the solution can be obtained in this way. Be it initially through perturbation theory or perturbatively corrected coupled cluster, the fully nonperturbative treatment of three-particle and three-hole excitations in midmass (and eventually heavy-mass) nuclei via EC will be superior to presently available implementations of many-body expansion methods.