Multi-reference many-body perturbation theory for nuclei

Abstract

In spite of missing dynamical correlations, the projected generator coordinate method (PGCM) was recently shown to be a suitable method to tackle the low-lying spectroscopy of complex nuclei. Still, describing absolute binding energies and reaching high accuracy eventually requires the inclusion of dynamical correlations on top of the PGCM. In this context, the present work discusses the first realistic results of a novel multi-reference perturbation theory (PGCM-PT) that can do so within a symmetry-conserving scheme for both ground and low-lying excited states. First, proof-of-principle calculations in a small ($e_{\mathrm {max}}=4$) model space demonstrate that exact binding energies of closed- (${}^{16}\mathrm {O}$) and open-shell (${}^{18}\mathrm {O}$, ${}^{20}\mathrm {Ne}$) nuclei are reproduced within 0.5–$1.5\%$ at second order, i.e. through PGCM-PT(2). Moreover, profiting from the pre-processing of the Hamiltonian via multi-reference in-medium similarity renormalization group transformations, PGCM-PT(2) can reach converged values within smaller model spaces than with an unevolved Hamiltonian. Doing so, dynamical correlations captured by PGCM-PT(2) are shown to bring essential corrections to low-lying excitation energies that become too dilated at leading order, i.e., at the strict PGCM level. The present work is laying the foundations for a better understanding of the optimal way to grasp static and dynamical correlations in a consistent fashion, with the aim of accurately describing ground and excited states of complex nuclei via ab initio many-body methods.

Appendices

Anti-symmetry reduction

From a technical viewpoint, and as extensively explained in Paper I, PGCM-PT(2) calculations rely on solving a large-scale linear problem of the form

$$\begin{aligned} \mathbf{M } \mathbf{a }&= -{\mathbf{h }_1} \, . \end{aligned}$$

(1)

The linear problem relates to excitations I of several non-orthogonal Hartree-Fock-Bogoliubov vacua, each of which is defined by a set of quasi-particle creation operators, i.e. a rank-n excitation is defined through a set of n quasi-particle labels $I\sim (k_{i_1}, \cdots , k_{i_n})$. Correspondingly, the problem is initially expressed in terms of unrestricted sets of quasi-particle indices. Still, anti-commutation relations of quasi-particle creation operators imply that ${\mathbf {M}},\ {\mathbf {a}}$ and ${\mathbf{h }_1}$ are anti-symmetric with respect to the permutations of quasi-particle indices. This can be exploited to reduce the effective dimensionality of the linear system.

Given a rank-n excitation $I\sim (k_{i_1}, \cdots , k_{i_n})$ on a given Bogoliubov state, the set ${\mathcal {I}} \equiv \{\tau ( I)\}_{\tau \in {\mathcal {S}}_n}$ of $|{\mathcal {I}}| \equiv n!$ permutations of the quasi-particle indices of $I$ needs to be considered. For a pair $(I,J)$ of excitations and two permutations $(\tau ,\tau ')$ applicable to $I$ and $J$, the antisymmetry properties are given by

$$\begin{aligned} M _{pIqJ}&= \epsilon (\tau ) \epsilon (\tau ') M_{p\tau (I) q\tau '(J)} \, , \end{aligned}$$

(2a)

$$\begin{aligned} a^{J}(q)&=\epsilon (\tau ) a^{\tau (J)}(q) \, , \end{aligned}$$

(2b)

$$\begin{aligned} h_1^{I}(p)&= \epsilon (\tau ') h_1^{\tau '(I)}(p) \, , \end{aligned}$$

(2c)

where $\epsilon (\tau )$ denotes the signature of the permutation $\tau$. First, these antisymmetry properties trivially imply that excitations with repeated quasi-particle indices can be excluded from the basis. Second, the set of excitations ${\mathcal {I}}$ corresponding to one another via a change of the quasi-particle ordering can all be tracked through the one representative ${{\bar{I}}}$ of ${\mathcal {I}}$ characterized by a strictly increasing ordering of the quasi-particle indices $k_1< \cdots < k_n$. Writing Eq. (1) for such an external ordered excitation ${{\bar{I}}}$

$$\begin{aligned} \sum _{q} \sum _{{\mathcal {J}}} \sum _{J\in {\mathcal {J}}} M _{p{{\bar{I}}} qJ} a^{J}(q)&= - h_1^{{{\bar{I}}}}(p) \, , \end{aligned}$$

(3)

the internal sum is split such that, with the help of Eq. (2), $|{\mathcal {J}}|$ equivalent terms are generated that eventually yield the reduced form

$$\begin{aligned} \sum _{q} \sum _{{\mathcal {J}}} |{\mathcal {J}}| M _{p{{\bar{I}}}\ q {{\bar{J}}}} a^{{{\bar{J}}}}(q)&= - h_1^{{{\bar{I}}}}(p) \, . \end{aligned}$$

(4)

In order to maintain the Hermiticity of the reduced matrix one further left-multiplies the equation by $\sqrt{|}{\mathcal {I}}|$ such that the final form

$$\begin{aligned} \sum _{q} \sum _{{\mathcal {J}}} \sqrt{|{\mathcal {J}}|}\sqrt{|{\mathcal {I}}|} M _{p{{\bar{I}}} q{{\bar{J}}}} \sqrt{|{\mathcal {J}}|} a^{{{\bar{J}}}}(q)&= -\sqrt{|{\mathcal {I}}|} h_1^{{{\bar{I}}}}(p) \, , \end{aligned}$$

(5)

naturally leads to a trivial redefinition of the reduced matrix and vectors through the inclusion of the combinatorial factors. In the following, the above reduction process is assumed such that the effective working basis only includes excitations characterized by quasi-particle indices in a strictly increasing order. For example, exploiting the anti-symmetry for the dominant double (i.e. 4 quasi-particle) excitations reduces the number of associated matrix elements by a factor of $24^2$.

Solution of the linear problem

Finding the numerical solution of Eq. (1) is delicate due to the non-orthogonality of the many-body states used to represent it. Thus, a careful handling of zeros in the norm eigenvalues is typically necessary to avoid instabilities while solving the equation. In the following, techniques of increasing sophistication are progressively introduced in order to eventually motivate the use of the iterative MINRES-QLP algorithm.

1.1 Exact SVD-based solution

The pedestrian way to solve the linear system can be summarized in three steps: (i) diagonalize the norm matrix to transform the equation into an orthonormal basis, (ii) diagonalize the Hamiltonian matrix in that basis and (iii) finally invert the problem. This strategy is essentially the same as the one used in PGCM to solve the HWG equation (see Appendix A of Paper II).

The norm matrix (see Paper I) is first decomposed by projecting on the range of ${\mathbf {N}}$ via a singular-value decomposition (SVD)

$$\begin{aligned} {\mathbf {N}} = {\mathbf {X}}{\mathbf {I}} {\mathbf {X}}^\dag \,, \end{aligned}$$

(6)

where ${\mathbf {X}}$ is unitary. Matrix I gathers the singular values whose smallest representatives can eventually be discarded. Correspondingly, ${\mathbf {M}}$, ${\mathbf {a}}$ and ${\mathbf {h}}_1$ are transformed into the resulting orthogonal basis

$$\begin{aligned} {\tilde{\mathbf {M}}}&\equiv {\mathbf {X}} {\mathbf {M}} {\mathbf {X}}^\dag \, , \end{aligned}$$

(7a)

$$\begin{aligned} {\tilde{\mathbf {a}}}&\equiv {\mathbf {X}}^\dag {\mathbf {a}} \, , \end{aligned}$$

(7b)

$$\begin{aligned} {\tilde{\mathbf {h}}}_1&\equiv {\mathbf {X}}^\dag {\mathbf {h}}_1 \, , \end{aligned}$$

(7c)

such that the linear problem equivalently reads

$$\begin{aligned} {\tilde{\mathbf {M}}} {\tilde{\mathbf {a}}} = - {\tilde{\mathbf {h}}}_1 \,. \end{aligned}$$

(8)

The solution of this system is then found by diagonalizing ${\tilde{\mathbf {M}}}$

$$\begin{aligned} {\varvec{\Delta }} = {\mathbf {Y}}^\dag {{\tilde{\mathbf {M}}}} {\mathbf {Y}} \,, \end{aligned}$$

(9)

such that, similarly to canonical MBPT, the system is inverted in the basis where ${{\tilde{\mathbf {M}}}}$ is diagonal to obtain the second-order energy in the form

$$\begin{aligned} E^{(2)} = - \mathbf { h_1}^\dag {\mathbf {X}}{\mathbf {Y}} {\varvec{\Delta }}^{-1}{\mathbf {Y}}^\dag {\mathbf {X}}^\dag \mathbf { h_1} \,. \end{aligned}$$

(10)

In principle, the projection on the range of ${\mathbf {N}}$ is not necessary to solve the system. However, in numerical applications, the coupling between spurious eigenvalues of ${\mathbf {N}}$ and large eigenvalues of ${\mathbf {M}}$ can arise and the explicit removal of the redundancies is often necessary. Eventually, full diagonalization is anyway not feasible for the large matrices encountered in realistic PGCM-PT(2) calculations (contrary to the PGCM step where the typical dimensions are sufficiently small) such that other methods need to be designed to solve the problem.

As an example, the distribution of the eigenvalues of ${\mathbf {N}}$ and ${\mathbf {M}}$ obtained from a PHFB-PT(2) calculation of ${}^{20}\mathrm {Ne}$ in a small model-space is displayed in Fig. 17. Two observations can be made

• The eigenvalue distributions of both matrices are very close up to a scaling factor. In particular, as expected, their (numerical) kernels have the same dimension.

• The kernel’s dimension is small compared to the matrices’ dimension, and all eigenvalues outside the kernel have the same magnitude. This prevents us from using truncated SVD approaches in larger model spaces.

Although the GCM mixing enlarges the kernel of the PGCM-PT(2) matrices compared to PHFB-PT(2) due to the partial linear dependencies of the added HFB states, a large number of independent configurations is still present in that case too.

Fig. 17

Distribution of eigenvalues of ${\mathbf {N}}$ and ${\mathbf {M}}$ matrices for ${}^{20}\mathrm {Ne}$. The calculation is performed with a two-body $\chi$EFT Hamiltonian, $\lambda _{\mathrm {srg}}=1.88\,\mathrm {fm}^{-1}$, $\hbar \omega =20\,\mathrm {MeV}$ and $e_{\mathrm {max}}=2$

1.2 Pivoted factorizations

Factorization algorithms can be applied in order to remove spurious eigenvalues without paying the price of fully diagonalizing the norm and Hamiltonian matrices. Typical examples are pivoted QR [30] and QLP [31] factorizations, which are briefly discussed in the following.

1.2.1 Pivoted QR

An arbitrary $n\times m$ complex matrix $A$ can be decomposed according to

$$\begin{aligned} A D = Q R \,, \end{aligned}$$

(11)

where $D$ is obtained via a permutation of the columns of $A$, $Q$ is a unitary matrix, and $R$ is an upper triangular matrix. The permutation $D$ is used to sort the diagonal entries of $R$ in decreasing order of magnitude. In this way, the kernel of $A$ corresponds to the last columns of $R$.

1.2.2 Pivoted QLP

Matrix A can be decomposed further by performing two successive pivoted QR decompositions, yielding the form

$$\begin{aligned} D' A D = QLP \end{aligned}$$

(12)

where $D, D'$ are permutation matrices, $Q, \ P$ are unitary matrices and $L$ is a lower triangular matrix. In particular, $L$ has the block-diagonal form

$$\begin{aligned} L = \begin{pmatrix} {{\tilde{L}}} &{}\quad 0\\ 0 &{}\quad 0 \end{pmatrix}, \end{aligned}$$

(13)

such that $Q$ and $P$ naturally block factorize $A$ into a full-rank part and its null-space.

Fig. 18

PHFB-PT correlation energy of ${}^{20}\mathrm {Ne}$ obtained for SVD, QR and QLP decompositions as function of the size of the excluded kernel in the decomposition. The calculation is performed with a two-body $\chi$EFT Hamiltonian, $\lambda _{\mathrm {srg}}=1.88\,\mathrm {fm}^{-1}$, $\hbar \omega =20\,\mathrm {MeV}$ and $e_{\mathrm {max}}=2$

1.2.3 Illustration of non-iterative solvers

In our case, pivoted QR/QLP factorizations can either be used directly on ${\mathbf {M}}$ to solve Eq. (1) or on the norm matrix in order to remove redundancies in the basis. In both cases, the symmetry of the matrices guarantees that the range and the kernel of both matrices are in direct sum. QLP factorization can thus be seen as a way to re-express the original problem in the range of ${\mathbf {M}}$ or ${\mathbf {N}}$. In practical applications, some tolerance must be used (as with SVD) to discard numerically small eigenvalues and disentangle the numerical kernel from the numerical range. The QLP factorization, although twice more expensive than the single QR decomposition, is found to be more stable and to better discard spurious eigenvalues.

Figure 18 shows a comparison of SVD, QR and QLP decompositions in a PHFB-PT(2) calculation. Errors estimated via $\delta E\equiv \Vert {\mathbf {M}} {\mathbf {a}} + {\mathbf {h}}_1\Vert \Vert {\mathbf {a}}\Vert$ (see Sect. B.3.2) are also represented on the figure. The three methods are in very good agreement with vanishing errors when nearly all the space is kept in the calculation. However, discrepancies arise when the truncation is performed according to the magnitude of the diagonal elements of the decomposition. While the SVD is by far the most reliable method, the QLP decomposition significantly improves the correlation energy with respect to the simpler QR decomposition and reduces the corresponding error for a fixed kernel dimension. Eventually, the reduced cost of QLP/QR decompositions compared to SVD, especially in their sparse version, make them well suited to large-scale calculations.

1.3 Iterative solvers

The QLP decomposition introduced above is still not applicable to very large matrices due to the runtime and memory requirements. An alternative solution is to use an iterative method, preferably exploiting the symmetry of the input matrix. Among various available solvers, the MINRES algorithm [32] finds the minimum-residual solution to $||{\mathbf {M}} {\mathbf {a}} + {\mathbf {h}}_1||$ via QR factorizations in the Krylov space of ${\mathbf {M}}$. In the case of ill-conditioned problems, QR factorizations are replaced by QLP factorizations, and the corresponding algorithm is referred to as MINRES-QLP [33].

The benefit of iterative solvers compared to exact decompositions is that, in the former, QR and QLP factorizations are performed in a Krylov subspace of the matrix. At iteration $k$, the problem is of dimension $k\times k$, where $k$ is usually much smaller than the original matrix dimensions. This results in both runtime and memory savings, at the cost of solving the system only approximately.

1.3.1 Preconditioning of the linear system

The number of iterations required by the solvers strongly depends on the eigenvalue distribution of the linear system under consideration. Typically, systems where the eigenvalues are clustered will have a faster convergence than systems with a spread-out spectrum. The spread of the eigenvalues can be altered with preconditioning techniques that amounts to finding equivalent systems with different (generally much smaller) condition numbers.

In this subsection, the compatible symmetric system

$$\begin{aligned} {\mathbf {M}} {\mathbf {x}} = -{\mathbf {h}}_1 \,, \end{aligned}$$

(14)

is considered. Let ${\mathbf {A}} = {\mathbf {C}} {\mathbf {C}}^T$ be a positive definite matrix. The solution of the initial system can be deduced from the solution of the preconditioned system

$$\begin{aligned} {\mathbf {C}}^{-1} {\mathbf {M}} {\mathbf {C}}^{-T} {\mathbf {y}} = - {\mathbf {C}}^{-1} {\mathbf {h}}_1 \,. \end{aligned}$$

(15)

Whereas various techniques are available to build an appropriate matrix ${\mathbf {A}}$, designing efficient preconditioners is still an active field of research [34]. There is no perfect preconditioner, and finding the trade-off between effectiveness and computational cost heavily relies on heuristics. Furthermore, for systems only known up to a given precision, preconditioners can artificially magnify eigenvalues that are numerically close to zero. Thus, a slower convergence with the possibility to stop the iterations before the appearance of spurious divergences might be preferable. Eventually, solving several equivalent systems simultaneously can make it easier to identify problematic features, and discrepancies between different solutions can be used as uncertainty estimates in the resolution.

Matrix scaling. Matrix scaling is a type of preconditioning where the preconditioner is a diagonal matrix ${\mathbf {A}} \equiv {\mathbf {D}}$ such that the equivalent system reads

$$\begin{aligned} {\mathbf {D}}^{-\frac{1}{2}} {\mathbf {M}} {\mathbf {D}}^{-\frac{1}{2}} {\mathbf {y}} = - {\mathbf {D}}^{-\frac{1}{2}} {\mathbf {h}}_1 \,. \end{aligned}$$

(16)

If the scaled matrix

$$\begin{aligned} {{\tilde{\mathbf {M}}}} \equiv {\mathbf {D}}{\mathbf {M}}{\mathbf {D}} \end{aligned}$$

(17)

is better conditioned than ${\mathbf {M}}$, then the solution to the initial system can be found in fewer iterations. For a diagonally dominant matrix ${\mathbf {M}}$, scaling the matrix with its own diagonal elements will reduce its condition number. The binormalization method detailed in Ref. [35] amounts to scaling all rows and columns to unit norm, which can yield significantly better results at low cost. A stochastic matrix-free variant [36] allows one to efficiently apply this method for abstract linear operators that are not necessarily defined explicitly by a matrix–vector product. Below, the stochastic binormalization preconditioner is denoted as SBIN.

Incomplete Cholesky decomposition. For a sparse positive definite matrix ${\mathbf {N}}$, an approximate Cholesky factorization preserving the sparsity pattern of the original matrix can be computed as

$$\begin{aligned} {\mathbf {N}} \sim {\mathbf {L}} {\mathbf {L}}^T, \end{aligned}$$

(18)

with ${\mathbf {L}}$ a (sparse) lower triangular matrix. A variant of Cholesky factorization applicable to positive indefinite matrices can be applied directly on the norm matrix ${\mathbf {N}}$. Since ${\mathbf {N}}$ and ${\mathbf {M}}$ have similar eigenvalue spread, eigenvalues of the system preconditioned with ${\mathbf {L}}{\mathbf {L}}^T$ will be much more clustered than those of the original system, hence separating the useful directions of the problem from the rest of the Hilbert space. Below, the incomplete Cholesky preconditioner is denoted as IC0.

Norm preconditioning. In some cases, spurious eigenvalues in the linear system can couple to physical modes and prevent any convergence of the iterative solvers. In this case, clustering the eigenvalues via preconditioning techniques is counterproductive as spurious modes are given an equivalent amplitude to physical ones. When this happens, it is preferable to amplify the separation of scale between numerically small and large eigenvalues. Instead of manually removing redundancies in the norm matrix, there exists a simple way to reach the image of ${\mathbf {N}}$ without resorting to a decomposition: Instead of solving

$$\begin{aligned} {\mathbf {M}} {\mathbf {a}} = - {\mathbf {h}}_1 \,, \end{aligned}$$

(19)

one directly solves for ${\mathbf {N}}^{-1} {\mathbf {a}}$ inside the range of ${\mathbf {N}}$ via

$$\begin{aligned} {\mathbf {N}} {\mathbf {M}} {\mathbf {N}} \left( {\mathbf {N}}^{-1} {\mathbf {a}}\right) = - {\mathbf {N}} {\mathbf {h}}_1 \,. \end{aligned}$$

(20)

Even if ${\mathbf {N}}$ is singular, the fact that ${\mathbf {h}}_1 \text { and } {\mathbf {a}}$ live in the range of ${\mathbf {N}}$ by construction ensures that ${\mathbf {N}}^{-1} {\mathbf {a}}$ is well-defined. The procedure ensures that small numerical eigenvalues of ${\mathbf {N}}$, originating from collinear many-body basis vectors, are tamed down in ${\mathbf {N}} {\mathbf {M}} {\mathbf {N}}$. Furthermore, numerical errors in ${\mathbf {h}}_1$ are suppressed as well. Of course, in exact arithmetic, the two systems are equivalent. This method corresponds in fact to preconditioning the system with ${\mathbf {N}}^{-2}$. As mentioned, this slows down the convergence of the iterative procedure and must be kept for cases where the direct approach or the complex shift method (see below) do not provide accurate solutions. Below, the norm precondition is denoted as ${\mathbf {N}} {\mathbf {M}} {\mathbf {N}}$.

1.3.2 Error evaluation

Iterative methods may require a large number of iterations or even diverge due to numerical errors. In this subsection, a conservative bound to estimate the error on the computed second-order energy is developed.

Given an approximate solution of the system

$$\begin{aligned} {\mathbf {M}} {\mathbf {a}} = -{\mathbf {h}}_1 +{\mathbf {b}} \,, \end{aligned}$$

(21)

the second-order energy evaluated with Hylleraas’ functional reads

$$\begin{aligned} E^{(2)}&= {\mathbf {a}}^\dag {\mathbf {M}} {\mathbf {a}} + {\mathbf {h}}_1^\dag {\mathbf {a}} + {\mathbf {a}}^\dag {\mathbf {h}}_1 \nonumber \\&= {\mathbf {a}}^\dag {\mathbf {b}} + {\mathbf {h}}_1^\dag {\mathbf {a}} \, . \end{aligned}$$

(22)

The difference between this expression and the directly evaluated second-order energy reads

$$\begin{aligned} \delta E^{(2)} = {\mathbf {a}}^\dag \left( {\mathbf {M}}{\mathbf {a}} + {\mathbf {h}}_1\right) \,. \end{aligned}$$

(23)

Thus, a conservative error estimate on the second-order energy is given by

$$\begin{aligned} |\delta E|^{(2)} \equiv \Vert {\mathbf {a}}\Vert \Vert {\mathbf {M}} {\mathbf {a}} + {\mathbf {h}}_1\Vert \,. \end{aligned}$$

(24)

The quantity $|\delta E|^{(2)}$ vanishes for an exact solution and grows whenever $\Vert {\mathbf {a}} \Vert$ becomes too large, which generally occurs if ${\mathbf {M}}$ is badly conditioned. When the norm-preconditioning is used, the error estimate is obtained as

$$\begin{aligned} |\delta E|^{(2)} \equiv \Vert {\mathbf {N}}^{-1} {\mathbf {a}}\Vert \Vert {\mathbf {N}} {\mathbf {M}} {\mathbf {N}} \left( {\mathbf {N}}^{-1}{\mathbf {a}}\right) + {\mathbf {N}} {\mathbf {h}}_1\Vert \,. \end{aligned}$$

(25)

1.3.3 Stopping condition for the iterative solver

MINRES-QLP already implements by default its own stopping criterion based on the relative norm of the residuals

$$\begin{aligned} r\equiv \frac{\Vert {\mathbf {M}} {\mathbf {a}} + {\mathbf {h}}_1\Vert }{\Vert {\mathbf {a}}\Vert }. \end{aligned}$$

(26)

In the present case, elements of ${\mathbf {M}}$ and ${\mathbf {h}}_1$ are obtained after several computational steps such that round-off and discretization errors will alter the quality of the input matrices. Furthermore, a threshold on the magnitude of the matrix elements of ${\mathbf {M}}$ is employed to enforce the sparsity of the matrix. As such, iterations should be stopped when the residual errors are of the same order as the precision of the input matrix elements.

1.3.4 Illustration of iterative solvers

In order to illustrate the use of iterative solvers and preconditioning techniques, results obtained in ${}^{20}\mathrm {Ne}$ with $e_{\mathrm {max}}=2$ are shown in Fig. 19. One observes that the IC0 preconditioning significantly reduces the number of iterations needed to reach the converged value. Contrarily, the norm preconditioning tends to spread the eigenvalues of the system and therefore slows down the convergence. For a well-behaved system, applying the IC0 preconditioning to the original matrix is therefore the method of choice. In contrast, whenever spurious eigenvalues prevent the convergence of the iterative process, the IC0 preconditioning amplifies the problem. Such a case is shown Fig. 20 for the ground state of ${}^{18}\mathrm {O}$. Here, applying a combination of SBIN and norm preconditioning is necessary to converge the system to the SVD solution.

Fig. 19

Correlation energy (top) and corresponding error (bottom) at each MINRES-QLP iteration for the ground state of ${}^{20}\mathrm {Ne}$. Results obtained with combination of IC0 and norm preconditionings are compared to the exact solution obtained via SVD. Calculations are performed with a two-body $\chi$EFT Hamiltonian [13, 14], $\lambda _{\mathrm {srg}}=1.88\,\mathrm {fm}^{-1}$, $\hbar \omega =20\,\mathrm {MeV}$ and $e_{\mathrm {max}}=2$

Fig. 20

Correlation energy (top) and corresponding error (bottom) at each MINRES-QLP iteration for the ground state of ${}^{18}\mathrm {O}$. Results obtained with a combination of SBIN and norm preconditionings are compared to the exact solution obtained via SVD. Calculations are performed with a two-body $\chi$EFT Hamiltonian [13, 14], $\lambda _{\mathrm {srg}}=1.88\,\mathrm {fm}^{-1}$, $\hbar \omega =20\,\mathrm {MeV}$ and $e_{\mathrm {max}}=2$

Complex-shift method

1.1 Motivations

As it appears in Eq. (10), the second-order energy relies on the invertibility of ${\varvec{\Delta }}$ to generate non-zero energy denominators. However, the eigenvalues of ${\varvec{\Delta }}$ can vanish, which makes the calculation of the second-order energy unstable or even ill-defined. Multi-reference methods are indeed susceptible to this so-called intruder-state problem [37, 38].

One way to regularize these zeros is to introduce a diagonal imaginary shift in the eigenbasis of ${\mathbf {M}}$. The eigenvalues are thus replaced by

$$\begin{aligned} \bar{{\varvec{\Delta }}} \equiv {\varvec{\Delta }} + \imath \gamma I \,, \end{aligned}$$

(27)

or, equivalently, working in the original basis

$$\begin{aligned} \bar{{\mathbf {M}}} \equiv {\mathbf {M}} + \imath \gamma {\mathbf {N}} \,, \end{aligned}$$

(28)

which simply corresponds to adding a complex term to the unperturbed Hamiltonian $H_0$. The imaginary shift moves zero eigenvalues of ${\varvec{\Delta }}$ into the complex plane and provides a robust way to remove intruder-state divergences.

In this context, the second-order energy is eventually evaluated by simply taking the real part of the Hylleraas functional,

$$\begin{aligned} E^{(2)} = \mathfrak {R}\left[ {\mathbf {a}}^\dag \bar{{\mathbf {M}}} {\mathbf {a}} +{\mathbf {a}}^\dag {\mathbf {h}}_1 + {\mathbf {h}}_1^\dag {\mathbf {a}}\right] . \end{aligned}$$

(29)

1.2 Implementation in real arithmetic

Although an extension of MINRES-QLP has been developed to handle complex symmetric matrices [39, 40], it is possible to rewrite the complex PGCM-PT(2) equations as an enlarged real-valued system, for which the original MINRES-QLP algorithm can be applied directly. The system of equations

$$\begin{aligned} ({\mathbf {M}} + \imath \gamma {\mathbf {N}}) ({\mathbf {a}} + \imath {\mathbf {b}}) = - {\mathbf {h}}_1 \end{aligned}$$

(30)

is recast into a blockwise 2x2 real symmetric system

$$\begin{aligned} \begin{pmatrix} {\mathbf {M}} &{}\quad -\gamma {\mathbf {N}}\\ -\gamma {\mathbf {N}} &{}\quad -{\mathbf {M}} \end{pmatrix} \begin{pmatrix} {\mathbf {a}}\\ {\mathbf {b}} \end{pmatrix} = \begin{pmatrix} -{\mathbf {h}}_1\\ {\mathbf {0}} \end{pmatrix} \,. \end{aligned}$$

(31)

Since the matrices are real by default after projection, implementing the imaginary shift via an augmented real system is profitable to make use of the MINRES-QLP real symmetric solver instead of variants designed for complex matrices.

Note that preconditioning techniques such as SBIN and ${\mathbf {N}}$-IC0 are still applicable within the augmented system.

1.3 Illustration

Proceeding with the ${}^{20}\mathrm {Ne}$ test case ($e_{\text {max}}=2$), the effect of the complex shifts on the MINRES-QLP iterations with the IC0 preconditioning is illustrated in Fig. 21. In general, the complex shift tends to lower the correlation energy — in the limit of an infinite shift, the correlation energy vanishes. Thus, a bias is introduced that must be monitored. Eventually, the larger the complex shift, the faster the iterative procedure will converge (towards a biased result).

Fig. 21

Correlation energy (top panel) and corresponding estimated error for different values of the complex shift $\gamma$ as a function of the number of MINRES-QLP iterations in ${}^{20}\mathrm {Ne}$. Calculations are performed with a two-body $\chi$EFT Hamiltonian [13, 14], $\lambda _{\mathrm {srg}}=1.88\,\mathrm {fm}^{-1}$, $\hbar \omega =20\,\mathrm {MeV}$ and $e_{\mathrm {max}}=2$

In the case of ${}^{18}\mathrm {O}$, we can combine the norm and SBIN preconditioning (cf.  Fig. 20), with the complex shift, as pictured in Fig. 22. In contrast to ${}^{20}\mathrm {Ne}$, the complex shift with the norm preconditioning decreases the convergence speed in this case. Applying the complex shift without norm preconditioning is also possible, but does only lead to a stabilization of the result before the occurrence of a divergence.

In practical applications, the optimal shift depends on the interaction, the model space and the system under consideration. As the model space is enlarged, encountering small eigenvalues becomes more probable and the complex shift becomes necessary to smear out the contaminations. A shift $\gamma \in [10,20]$ MeV is well suited to remove spurious behaviors, with an estimated error of around $4\%$ on the correlation energy, as can be seen in Fig. 23. The difference between the results obtained with $\gamma =15$ MeV and $\gamma =4$ MeV²⁸ is used to estimate the bias due to the shift. Note that PHFB-PT(2) is more sensitive to intruder-state problems than PGCM-PT(2), hence the need to employ a larger shift to smooth out singularities on the energy curve. In practice, it is essential to use the same shift for all quantum states of a given nucleus to obtain a consistent bias in absolute binding energies that will largely cancel out in the excitation spectrum. The development of an extrapolation method towards $\gamma \rightarrow 0$ to correct for the bias due to the complex shift is left for a future study.

Fig. 22

Correlation energy (top panel) and corresponding estimated error for different values of the complex shift $\gamma$ as a function of the number of MINRES-QLP iterations in ${}^{18}\mathrm {O}$. Calculations are performed with a two-body $\chi$EFT Hamiltonian [13, 14], $\lambda _{\mathrm {srg}}=1.88\,\mathrm {fm}^{-1}$, $\hbar \omega =20\,\mathrm {MeV}$ and $e_{\mathrm {max}}=2$

Fig. 23

Correlation energy in ${}^{20}\mathrm {Ne}$ for a complex shift $\gamma =15\,\mathrm {MeV}$. Error bars associated to the effect of the shift correspond to the correlation energy with $\gamma =4\,\mathrm {MeV}$. The calculation is performed with a two-body $\chi$EFT Hamiltonian [13, 14], $\lambda _{\mathrm {srg}}=1.88\,\mathrm {fm}^{-1}$, $\hbar \omega =20\,\mathrm {MeV}$ and $e_{\mathrm {max}}=2$

Table 1 Runtime complexity and storage requirements for various resolution methods of the many-body problem. $n_{\text {proj}}$ denotes the number of gauge angles used for projections and $n_\text {gcm}$ the number of states used in the mixing

Discussion on numerics

1.1 Scaling

Any method to solve A-body Schrödinger equation’s comes with its numerical complexity and memory requirement. For a given basis size $N$ of the one-body Hilbert space, the naive polynomial scaling of runtime and storage of the methods discussed in the present work is displayed in Table 1. These asymptotic values are to be revised when particular symmetries are exploited in the many-body bases (e.g. spherical or axial symmetry reducing the size of the many-body tensors at play). Moreover, prefactors (ignored here) may play a significant role. Nevertheless, the table gives a fair idea of the asymptotic cost of the different many-body techniques.

As an example, the computational cost of each individual matrix element at play in PGCM-PT(2), which requires about 1000 vectorized elementary operations, makes the construction of the matrix the most time-consuming step in the calculation. In other words, the computation of the $N^8$ matrix elements dictates the overall complexity, given that the (approximate) sparsity of the matrix makes the cost of solving the linear system subleading. Similarly, even if BMBPT(3) has the same storage cost as HFB in principle, symmetry properties of the density matrices are used to drastically reduce the number of matrix elements that are needed at the HFB level. In general, the nominal complexity and storage requirement have to be balanced with the optimizations (vectorization, parallelization, compression techniques) that can be applied for a specific method, and they can play a decisive role in practical applications. Also, shape mixing through PGCM scales quadratically with the number of reference states, i.e. a PGCM (PGCM-PT(2)) calculation with 10 states is 100 times more costly than a PHFB (PHFB-PT(2)) calculation.

A selection of runtimes as a function of the one-body basis dimension is displayed in Fig. 24. For BMBPT(2,3) and PHFB, symmetry properties lower the effective complexity to $O(N^4)$. The main differences reside in the prefactor, which is, intiutively, larger for BMBPT(3). Note that the normal ordering of the Hamiltonian and the transformation to the quasi-particle basis are included in the runtime estimate.

Fig. 24

Timing of many-body methods as a function of the basis size of the one body Hilbert space. Projections are performed on total angular momentum $J$ with 24 gauge angles

1.2 Complexity reduction in PGCM-PT(2)

The multi-reference PGCM-PT(2) calculation is, in its naive formulation, significantly more costly than its single-reference counterparts. This is mainly due to the redundancies in the visited Hilbert space: many projected quasi-particle configurations play little to no role in the correlation energy or are redundant. This is even more true for large-scale applications where the multi-reference unperturbed state mixes many product HF(B) states. This naturally leads to the idea of reducing the dimensionality of the problem by selecting only relevant configurations (see Refs. [41, 42] for recent applications of this idea in nuclear physics and chemistry). In particular, the application of importance-truncation techniques in the context of non-perturbative methods [43] shows promising results that should be applicable to the present problem.

Several procedures to reduce the number of configurations in a controlled way are now briefly introduced, although not all of them have been implemented yet.

1.2.1 Norm-based importance truncation

Exact arithmetic The norm of a projected configuration is

$$\begin{aligned} n^I(p) \equiv \langle \varOmega ^I(p)|\varOmega ^I(p)\rangle \in [0,1] \end{aligned}$$

(32)

such that a configuration I for which $n^I(p)=0$ satisfies

$$\begin{aligned} |\varOmega ^I(p)\rangle =0 \,. \end{aligned}$$

(33)

Trivially, a null vector does not contribute to the linear system and can be safely removed from the calculation.

Approximate zeros Given a threshold $\epsilon _n > 0$, the norm-based importance-truncated problem is introduced by removing configurations $I$ with $n^I(p) < \epsilon _n$. The exact problem is obtained in the limit $\epsilon _n=0$. For now, this is the only method that has been implemented and applied to discard configurations in ${}^{18}\mathrm {O}$ at $e_{\text {max}}=6$. Although the number of configurations was divided by two (from $10^6$ to $5 {\cdot } 10^5$) by only keeping configurations whose norm reaches 2% of the maximal value, the induced error was shown to be less than 1%. A systematic study of the results obtained via this procedure still remains to be performed.

1.2.2 Hamiltonian-based importance truncation

Exact arithmetic The Hamiltonian matrix element of a projected configuration reads

$$\begin{aligned} h_1^I(p) \equiv \langle \varOmega ^I(p)| H_1 |\varTheta ^{(0)}\rangle \,. \end{aligned}$$

(34)

A configuration I for which $h_1^I(p)=0$ does not contribute to the linear system nor to the second-order energy, hence it can be safely removed from the calculation.

Approximate zeros The Hamiltonian-based importance-truncated problem is introduced by removing configurations $I$ satisfying $|h_1^I(p)| < \epsilon _h$, with $\epsilon _h > 0$. The exact problem is obtained in the limit $\epsilon _h=0$.

1.2.3 Energy-based importance truncation

The contribution of a configuration $I$ associated with the vacuum $|\varPhi (q)\rangle$ to the second-order correlation energy is

$$\begin{aligned} e^{(2)I}(q) = h_1^{I*}(q) a^{I}(q) \,. \end{aligned}$$

(35)

A configuration $(q,I)$ for which $e^{(2)I}(q)=0$ does not contribute to the correlation energy²⁹. Removing configurations based on the size of their contribution to the correlation energy corresponds to the method advocated in Refs. [41, 43]. The method is expected to lead to a substantial gain for a negligible error on the energy, although the impact on other observables must be checked as well. Of course, computing $e^{(2)I}(q)$ requires to solve the problem in the first place and is thus impractical. The idea is thus to evaluate the importance of a given configuration (q, I) by calculating an approximation to $e^{(2)I}(q)$ at a significantly reduced cost, which can typically be achieved by using BMBPT(2) based on the HFB vacuum $| \varPhi (q) \rangle$.