Excited states from eigenvector continuation: the anharmonic oscillator

Eigenvector continuation (EC) has recently attracted a lot attention in nuclear structure and reactions as a variational resummation tool for many-body expansions. While previous applications focused on ground-state energies, excited states can be accessed on equal footing. This work is dedicated to a detailed understanding of the emergence of excited states from the eigenvector continuation approach. For numerical applications the one-dimensional quartic anharmonic oscillator is investigated, which represents a strongly non-perturbative quantum system where the use of standard perturbation techniques break down. We discuss how different choices for the construction of the EC manifold affect the quality of the EC resummation and investigate in detail the results from EC for excited states compared to results from a full diagonalization as a function of the basis-space size.


Introduction
The solution of the nuclear many-body problem has seen significant progress over the last years, enabling a firstprinciples microscopic description based on nuclear Hamiltonians from chiral effective field theory combined with ab intio many-body approaches [1,2]. While the exact solution can be obtained from variational techniques such as configuration interaction (CI) or quantum Monte-Carlo (QMC) approaches, the exponential scaling of the underlying methods prevents its use in medium-mass systems. As an alternative, a diverse toolbox of many-body techniques have been designed that expand the exact solution around a suitably chosen A-body reference state at low polynomial cost at the price of sacrificing the variational character, e.g., many-body perturbation theory (MBPT) [3][4][5], the in-medium similarity renormalization group (IMSRG) [6][7][8], coupled-cluster (CC) theory [9][10][11], or self-consistent Green's function (SCGF) theory [12][13][14]. These developments allowed ab initio studies of atomic nuclei containing up to one hundred particles [14][15][16] as well as global calculations of the nuclear chart up to the iron region [17].
Even with powerful non-perturbative methods at hand, MBPT methods have undergone a major revival (see Ref. [5]) in nuclear structure theory due to the development of RG evolution techniques that make the nuclear Email addresses: margarida.companys@stud.tu-darmstadt.de (M. Companys Franzke), alexander.tichai@physik.tu-darmstadt.de (A. Tichai), kai.hebeler@physik.tu-darmstadt.de (K. Hebeler), schwenk@physik.tu-darmstadt.de (A. Schwenk) many-body problem computationally more tractable by softening nuclear forces to lower resolution [18,19]. Still, the convergence of perturbative expansions is not guaranteed and a reliable extraction of nuclear observables is additionally complicated in quasi-degenerate many-body systems, such as open-shell nuclei.
In the past, various resummation techniques were developed that allow to extract observables from a possibly divergent expansion. While such schemes have been applied successfully, their power usually relies on an a priori knowledge of the asymptotic behavior of the expansion, information that is rarely available in realistic applications. Consequently, practitioners desire a robust yet accurate tool that enables for extracting physical observables in a numerically reliable way. Recently, eigenvector continuation (EC) has been introduced for systems where the Hamiltonian H(c) admits a smooth dependence on some external parameter c ∈ R. The EC method is based on analytical continuation of the expansion outside its initial domain of convergence by performing several re-expansions of the Taylor series, thus, effectively shifting the reference point of the many-body expansion. This allows to explore the system at coupling values that are outside of the initial domain of convergence [20][21][22][23] (see also Ref. [24] for a similar resummation approach developed in quantum chemistry). The EC method has been applied in various manybody methods, but mainly restricted to the evaluation of ground-state properties. Moreover, EC has been used as a powerful and accurate emulator for quantifying theoretical uncertainties due to the nuclear Hamiltonian [25][26][27][28][29].
In this Letter, we investigate low-lying excited states in the EC framework. We consider the anharmonic oscillator as a strongly coupled benchmark system that is outside the range of applicability of standard perturbation theory calculations [30], and in all cases the low-lying spectrum can be benchmarked against the exact solution obtained from full diagonalization. We show that EC targeted at the ground state is able to access the excited states as well, and that this can be further improved by including in the EC information from perturbative corrections for the excited states.

Eigenvector continuation
The EC method is used to target physical systems that smoothly depend on an external parameter H(c). In many applications there exists a regime 0 ≤ c ≤ c crit < 1, where the many-body problem is easier to solve than for the target value c = 1. While physical observables, e.g., energy eigenvalues, can change significantly when c is varied, the eigenstates themselves often are less sensitive and remain in a low-dimensional manifold of the many-body Hilbert space when changing c to its physical value c = 1.
Consequently, the EC framework is performed in two successive steps. First a reference manifold containing N EC auxiliary states is constructed, (1) does not assume specific properties of the defining manybody wave functions. Following the previous work in Refs. [21,22], M EC is here constructed from PT wave functions for both ground and excited states on top of a non-degenerate reference state.
Perturbation theory starts from a partitioning of the initial Hamiltonian H = H 0 + H 1 into an unperturbed part H 0 and a perturbation H 1 . The 0-th order reference state is an eigenstate of the unperturbed Hamiltonian where the subscript n labels the ground and excited states. Introducing a parameter-dependent Hamiltonian H(c) ≡ H 0 + cH 1 the exact many-body wave function is written as an infinite Taylor expansion in terms of an auxiliary parameter c, where |Φ (p) n denotes the p-th order state correction.
Using a perturbative ansatz for the many-body wave functions the EC manifold can be re-expressed via the transformation . . .
thus, yielding an equivalent representation of the EC manifold where P is the maximum perturbation order considered in the wave function [21,22]. Finally, the explicit set of parameters {c i } does not enter in this setup and the dimension of M EC is set by the maximum perturbative order P . Consequently, the EC approach amounts to the numerical solution of a generalized eigenvalue problem [20][21][22], where H and N are the Hamiltonian and norm matrices expanded in the basis of PT state corrections, The quantity E denotes the P + 1 generalized eigenvalues and corresponds to the EC energies. By employing intermediate normalization, i.e., Φ n |Ψ = 1 with all basis states, the norm matrix fulfils N 0p = δ 0p and the Hamiltonian matrix entries within the first row (and column) correspond to the MBPT energy corrections, i.e., H 0p = H p0 = E (p+1) . However, when p, q = 0 the matrix element H pq ∼ c p+q contains many-body correlations up to order p + q, thus, going beyond a simple PT evaluation. Due to the final diagonalization the EC approach is intrinsically non-perturbative and resums correlations not present in a simple PT approach. Moreover, the diagonalization ensures the EC framework to be manifestly variational, such that going to higher orders guarantees an improvement in accuracy as opposed to most mediummass many-body frameworks applicable.
In practice, high-order energy and state corrections can be accessed using a recursive formulation of PT [4,31]. Computationally, all quantities are obtained from matrixvector multiplications of the interacting Hamiltonian represented in the unperturbed basis, i.e., in terms of eigenstates of H 0 .
For targeting excited states, we investigate two different strategies for EC. First, perturbative state corrections are evaluated for the many-body ground state and excited states are accessed as eigenvalues from the solution of the generalized eigenvalue problem by forming Hamilton and norm matrices in Eq. (7) with n = 0. Second, the acronym EC n is used to indicate that the PT state corrections are evaluated for the n-th excited state in the unperturbed spectrum |Φ (p) n and the EC matrices are constructed from these states. For EC 0 , this coincides with the first strategy.

Anharmonic oscillator
As an example we consider the one-dimensional harmonic oscillator with a quartic anharmonic term where m is the mass of the particle, ω the oscillator frequency, and p and x the momentum and position operators, respectively. The strength of the perturbation is controlled by c. In the following, we work in natural units with m and ω set to unity. For the perturbative treatment the anharmonic oscillator (aHO) Hamiltonian is partitioned according to where the unperturbed system corresponds to the harmonic oscillator (HO) that can be solved exactly, with an equidistant spectrum n = n + 1/2 with quantum number n = 0, 1, 2, . . .. In the following n max defines the maximum excitations n ≤ n max in the model space.
The seminal work by Bender and Wu revealed that the application of perturbation theory for the aHO Hamiltonian leads only to an asymptotic expansion for arbitrary coupling strength due to the presence of a branch-cut singularity at the origin [30]. The divergence of the energy corrections for different reference states is shown in Fig. 1. It is clear that the energy corrections are exponentially growing, thus, making an extraction of observables from the bare perturbation series impossible.
Consequently, the aHO Hamiltonian provides a nontrivial testcase to study the low-lying spectrum within the EC framework. We note that the spectrum of the unperturbed Hamiltonian is non-degenerate such that excited states can be targeted on equal footing without resorting to (quasi-)degenerate PT extensions that further require a diagonalization within a degenerate subspace. The reference states for the recursive PT evaluation are chosen as |Φ Since the aHO Hamiltonian is parity-conserving, eigenstates with different parities Π do not mix, and we can therefore consider the Hamiltonian and EC for positive and negative parity separately. For completeness, matrix elements of the perturbation operator in the unperturbed basis are given by

Sensitivity of the norm matrix
The quality of the numerical solution of a generalized eigenvalue problem is strongly sensitive to the condition κ of the norm matrix, where the {s i } correspond to the non-negative singular values of the norm matrix. In our cases the norm matrix arises from the overlap of perturbative state corrections that potentially admit (near) linear dependences. The larger the condition number, the more singular is the norm matrix, such that the presence of (near) linear dependencies make the inversion of the norm matrix numerically challenging. As a consequence, special care needs to be taken to resolve this at the level of a desired accuracy. Figure 2 shows the norm matrix in the basis of PT state corrections for model-space dimension n max = 50. While at low order PT basis vectors point in independent directions, at higher order in the PT expansion linear dependences occur, i.e., Φ norm matrix has large entries near the diagonal. The presence of such states induces very small singular values of the norm matrix as can be seen from the lower panel of Fig 2. We observe a rapid falloff of singular values such that s i < 10 −10 for p > 25. Therefore, high-order state corrections contain a lot of redundant information manifesting in off-diagonal norm matrix entries close to unity (in absolute value) or, equivalently, very small singular values. In the case of n max = 50 the observed condition number is given by κ 10 16 hinting at a numerically illconditioned problem. In our calculations this problem was resolved by performing an explicit re-orthogonalization of the PT statecorrection basis using Householder transformations yielding a new basis of PT corrections. Within this basis the norm matrix reduces to the unit matrix and the generalized eigenvalue problem in Eq. (6) reduces to a standard eigenvalue problem for the transformed Hamiltonian. Trivially, the condition number is given by κ = 1 in that case. In principle alternative orthogonalization techniques such as (modified) Gram-Schmidt or QR-decompositions can be applied. In our applications we did not find a strong sensitivity on the employed method for the range of desired accuracy. Figure 3 displays the lowest five eigenstates obtained from EC 0 for model-space dimension n max = 50 in each parity subspace. The low-lying exact eigenenergies are indicated by black lines. We observe in both parity subspaces similar convergence properties. In all cases the exact energy is recovered in a numerically stable way for ground and excited states. However, the number of EC basis states required to reach convergence changes drastically for different states in the spectrum. While good accuracy for the ground state is obtained within p ≈ 5, excited states require larger dimensions, e.g., p ≈ 20 for k = 1. Note that the subspace probed by M PT at order p = n max is identical to the full basis space such that the EC energies have to coincide with the full diagonalization results. For k ≥ 2 we observe a significantly slower convergence and some untypical convergence patterns. In particular, already for k ≥ 3 the full model space needs to be exhausted to reproduce the exact energies. Nevertheless, convergence is in all cases monotonic as imposed by the underlying variational principle. Next the EC k framework is benchmarked for k = 0. Figure 4 shows the low-lying spectrum in both parity blocks obtained from EC 1 , i.e., using the first excited state in the unperturbed spectrum as reference for the subsequent PT treatment. We observe that for p 40 the lowest EC state corresponds to the first excited state of the exact spectrum, i.e., the selected subspace does not contain information on the ground-state wave function. Moreover, in the EC k -basis convergence for the k-th excited state is most rapid. This becomes apparent when comparing Figs. 3 and 4, where the EC 1 convergence of the first excited states is faster than in the case of EC 0 . This is consistent with the fact that the construction of the basis was optimized for a particular configuration. We have checked that these observations hold for higher excited states as well. Eventually, for p → n max the EC manifold exhausts the full configuration space. Consequently, for p ≈ 45 the EC states become lower in energy since the variational principle enforces the reproduction of the exact spectrum. At p = n max the full diagonalization energies and results from EC k coincide for arbitrary k since EC provides a reparametrization of the same Hilbert space in a transformed basis. We expect the particular value of p < n max where the EC spectrum starts approaching the groundstate energy to be connected to the particular form of the anharmonic perturbation employed in this work. While

Comparison of configuration bases
As discussed previously, the EC framework provides a way of defining a low-dimensional manifold of configurations which is used for a subsequent diagonalization. Naturally, the quality of the EC results will depend on the quality of the selected subspace for a given observable. Therefore, we compare the convergence of EC with the where |Φ k = e k denotes the k-th unit vector. In manybody applications the states |Φ k will correspond to simple Slater determinants whereas the PT corrections are of multi-configurational character, i.e., superposition of several Slater determinants. Figure 6 displays the convergence with the model space dimension of the low-lying spectrum in both parity subspaces. The convergence as a function of basis size is much faster and more natural than in the case of an EC-selected subspace. In particular, the two lowest states for each parity can be accurately extracted using p 10. Even for the fourth excited state less than 20 configurations are needed to reproduce the full diagonalization result. Comparing this to the quality of low-order EC calculations indicates that in the present case the canonical basis performs much better than any of the choices in the EC approach. This is not surprising, since the quality of the EC manifold is directly connected to the quality of the PT state corrections. Since we are working in a non-perturbative environment, i.e., exponentially divergent perturbation series, PT-based wave functions will not accurately approximate the true eigenstates. In fact, PT state corrections seem not to pick optimal subspaces of the Hilbert space in our case and, therefore, give slower convergence in EC applications. This is to be contrasted with previous applications using softer Hamiltonians in many-body problems, where the EC was able to efficiently select the relevant low-dimensional subspace out of a large set of Slater determinants [21].

Conclusion
In this Letter, EC was used to extract the energies of ground and excited states for the case of the anharmonic oscillator with quartic perturbation. Even in presence of strongly divergent PT expansions the EC-resummed energies robustly converges towards the full diagonalization results where the convergence is fastest for the ground-state expansion. The convergence of the EC approach for excited states can be improved by using PT state corrections for the excited states themselves, thus, informing the EC basis on the properties of the target wave function. The EC 0 approach is limited to excited states in the same symmetry class as the fully interacting ground state. While in aHO applications parity conservation separates the spectrum into two parts, most large-scale ab initio applications are performed in a symmetry-restricted way, such that only states with the same total angular-momentum are accessible from such calculations, e.g., low-lying J Π = 0 + states in the case of even-even nuclei. We note that the design of EC k in realistic applications is technically more complicated due to the presence of degeneracies in the unperturbed spectrum making the formulation of PT for excited states more complex.
The EC approach can be seen as a selected CI approach employing a contracted basis from PT-based wave function corrections. As such it suffers from the same limitations as truncated CI such as the lack of size-extensivity at finite truncation order. In the future, we will investigate size-extensivity properties, which will require the evaluation of low-order EC approximations for large system sizes. These developments will be supported by deriving a diagrammatic expansion of the EC formalism at low order which circumvents the formation of the exponentially large configuration space.