11institutetext:  IRFU, CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette,
France 22institutetext:  KU Leuven, Department of Physics and Astronomy,
Instituut voor Kern- en Stralingsfysica, 3001 Leuven, Belgium 33institutetext:
CEA, DAM, DIF, 91297 Arpajon, France 44institutetext:  Université Paris-
Saclay, CEA, Laboratoire Matière en Conditions Extrêmes, 91680, Bruyères-le-
Châtel, France 55institutetext:  Departamento de Física Teórica, Universidad
Autónoma de Madrid, 28049 Madrid, Spain 66institutetext:  Facility for Rare
Isotope Beams, Michigan State University, East Lansing, MI 48824-1321, USA
77institutetext:  Department of Physics and Astronomy, Michigan State
University, East Lansing, MI 48824-1321, USA 88institutetext:  Centro de
Investigación Avanzada en Física Fundamental-CIAFF-UAM, 28049 Madrid, Spain
99institutetext:  Institut für Kernphysik, Technische Universität Darmstadt,
64289 Darmstadt, Germany 1010institutetext:  Helmholtz Forschungsakademie
Hessen für FAIR, GSI Helmholtzzentrum, 64289 Darmstadt, Germany
1111institutetext:  School of Physics and Astronomy, Sun Yat-sen University,
Zhuhai 519082, P.R. China

em:mf em:td em:jpe em:bb em:hh em:tr em:rr em:jmy em:vs

# Multi-reference many-body perturbation theory for nuclei

III. Ab initio calculations at second order in PGCM-PT

M. Frosiniad:saclay,em:mf  T. Duguetad:saclay,ad:kul,em:td  J.-P.
Ebranad:dam,ad:fakedam,em:jpe  B. Ballyad:dft,em:bb  H.
Hergertad:msu1,ad:msu2,em:hh  T.R. Rodríguezad:dft,ad:cdiaf,em:tr  R.
Rothad:tud,ad:hfhf,em:rr  J. M. Yaoad:sysu,em:jmy  V. Somàad:saclay,em:vs
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(Received: / Revised version: date)

###### 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.

##  1 Introduction

The recent breaking of ab initio calculations away from p-shell nuclei into
the realm of mid-mass nuclei has been made possible by the formulation and
implementation of so-called many-body expansion methods. Because of their
polynomial scaling with system size, expansion methods provide the best
candidates yet to extend the reach of ab initio calculations to even heavier
nuclei. However, and as explained in the introduction to the first paper of
the present series frosiniI , hereafter referred to as Paper I, a current
challenge concerns the optimal way to consistently capture both static and
dynamical correlations within such methods. While doubly closed-shell nuclei
are dominated by (weak) dynamical correlations that are efficiently grasped
through a coherent sum of (mostly low-rank) particle-hole excitations of a
symmetry-conserving unperturbed product state, open-shell nuclei display
strong static correlations that cannot be conveniently accounted for in this
way. This results in the necessity to design expansion methods based on more
general unperturbed states that can already capture static correlations.

A good candidate to provide appropriate unperturbed states is the projected
generator coordinate method (PGCM). The main conclusion of the second paper of
the present series frosiniII  (Paper II in the following) is that the PGCM is
suitable to address the low-lying spectroscopy of complex nuclei within
reasonable theoretical uncertainties in spite of missing dynamical
correlations. For instance, the energy spectrum and electric multipole
transition strengths of the low-lying parity-doublet bands in 20Ne were
reproduced by taking into account both quadrupole and octupole collective
fluctuations.

Still, describing absolute binding energies, accounting consistently for a
wide range of spectroscopic observables, tackling a large class of nuclei
displaying different characteristics and achieving high accuracy eventually
requires the inclusion of dynamical correlations on top of the PGCM. This
coherent incorporation is made possible by expanding the wave operator
$\Omega$ connecting the PGCM state to the exact eigenstate via the novel
multi-reference perturbation theory (PGCM-PT) formulated in Paper I. Doing so,
PGCM-PT embeds, for the first time, the PGCM within a systematic symmetry-
conserving expansion method.

The objective of the present work, the third paper of the series, is to
discuss first proof-of-principle results of second-order PGCM-PT, i.e. PGCM-
PT(2), calculations in three selected nuclei, namely the doubly closed-shell
${}^{16}\mathrm{O}$, the singly open-shell ${}^{18}\mathrm{O}$ and the doubly
open-shell ${}^{20}\mathrm{Ne}$ that was studied at length at the PGCM level
in Paper II.

In addition to displaying the first set of PGCM-PT(2) results, the goal of the
present work is to do so while exploiting an additional degree of freedom at
our disposal in quantum many-body calculations, i.e. the possible pre-
processing of the Hamiltonian, e.g., via unitary transformations generated by
nucleus-dependent in-medium similarity renormalization group (IMSRG)
evolution. While nucleus-independent vacuum similarity renormalization group
(VSRG) transformations of the Hamiltonian have already become a standard tool
to pre-sum ultra-violet (UV) dynamical correlations by decoupling low- and
high-momentum modes, nucleus-dependent transformations can be exploited more
systematically to pre-sum infra-red (IR) dynamical correlations.

Figure 1: (color online) Schematic workflow of expansion many-body methods
(vertical axis) versus potential pre-processings of the Hamiltonian
(horizontal axis). Unitary vacuum (in-medium) similarity renormalization group
transformations denote a nucleus-independent (nucleus-dependent) pre-
processing of the Hamiltonian.

The single-reference IMSRG (SR-IMSRG) method Tsukiyama:2010rj ;
Hergert:2015awm  applicable to closed-shell systems can in principle fully
decouple the unperturbed product state from the rest of the Hilbert space,
i.e. from the ${\cal Q}$ space, and thereby make it the actual ground-state of
the pre-processed Hamiltonian at the end of the flow. In this case, the wave
operator eventually becomes nothing but the identity operator and the
expansion method on top of the unperturbed state is trivial. The more general
multi-reference IMSRG (MR-IMSRG) method Hergert:2014iaa ; Hergert:2016iju ;
Hergert:2016etg  applicable to all nuclei cannot, even in principle, fully
decouple the PGCM unperturbed state from the associated ${\cal Q}$ space such
that non-zero dynamical correlations remain to be included via a non-trivial,
e.g. PGCM-PT, wave operator111While the IMSRG constitutes per se a method to
solve Schrödinger’s equation when the SR-IMSRG implementation can be applied,
it is not the case for the MR-IMSRG approach that can only be seen as a pre-
processing of the Hamiltonian on top of which an appropriate many-body method
must be applied.. While the impact of these remaining dynamical correlations
on absolute energies may be small, we will show that their inclusion is
important for a proper description of low-lying excitation spectra.

Eventually, a clear picture will emerge that is schematically illustrated in
Fig. 1. Three complementary levers must be consistently exploited to
efficiently capture correlations within many-body expansion methods in order
to describe (complex) nuclei

  1. 1.

the pre-processing of the Hamiltonian,

  2. 2.

the possibly non-trivial nature of the unperturbed state,

  3. 3.

the rationale of the expansion.

While each lever is best suited to capture a certain category of correlations,
the latter are not orthogonal to one another such that the ideal way to share
the load is unclear and will require extensive trial-and-error in the future.
The present work wishes to contribute to this long-term endeavor.

The present paper is organized as follows. Section 2.1 details the results
obtained in ${}^{16}\mathrm{O}$, ${}^{18}\mathrm{O}$ and ${}^{20}\mathrm{Ne}$
on the basis of (the two-body part of) a chiral effective field theory
($\chi$EFT) Hamiltonian evolved via VSRG. Section 3 then elaborates on the
impact of further pre-processing the Hamiltonian via MR-IMSRG transformations
on the results. Conclusions and future perspectives are eventually discussed
in Sec. 4. A set of technical appendices provides additional details about the
numerical solution of the large-scale linear system of equations at play in
PGCM-PT(2) calculations.

##  2 Calculations with VSRG pre-processing

The reader is referred to Paper I for all necessary details about the PGCM-PT
formalism as well as to Refs. Bogner:2009bt ; PhysRevLett.107.072501 ;
PhysRevC.90.024325  and Refs. Hergert:2015awm ; Hergert:2016iju ;
Hergert:2016etg  for vacuum and in-medium IMSRG methods, respectively. The
PGCM-PT(2) solver is built on top of an axially-deformed Hartree Fock
Bogoliubov (HFB) code frosini21e  and a consistent PGCM solver frosini21f
allowing for the projections on good particle number, angular momentum and
parity. All notations used below are consistent with those introduced in
Papers I and II that should be consulted for reference.

###  2.1 Numerical setting

Proof-of-principle calculations are performed using the spherical harmonic
oscillator (HO) basis of the one-body Hilbert space ${\cal H}_{1}$
characterized by an oscillator frequency $\hbar\omega=20\text{ MeV}$ and $5$
oscillator shells ($e_{\text{max}}=4$).

The next-to-next-to-next-to-leading order (N3LO) $\chi$EFT Hamiltonian
introduced in Refs. Huther_2020 ; Entem:2017gor  and evolved via VSRG to the
low-momentum resolution scale $\lambda_{\text{vsrg}}=1.88\text{ fm}^{-1}$ is
employed. Thus, UV dynamical correlations are already processed via the VSRG
decoupling of low- and high-momentum modes.

In these proof-of-principle calculations, only the two-body part of the
evolved Hamiltonian is actually retained. Thus, the goal is not to reproduce
experimental data but rather to benchmark PGCM-PT(2) results against those
obtained from full configuration interaction (FCI) calculations in the same
$e_{\text{max}}=4$ space. The FCI calculations rely on a sequence of
$N_{\mathrm{max}}$-truncated many-body Hilbert spaces up to
$N_{\mathrm{max}}=8$ embedded into the FCI space defined by the
$e_{\text{max}}=4$ truncation on the single-particle basis. The results are
extrapolated to the full $e_{\text{max}}=4$ model space limit such that FCI
results come with an uncertainty associated with this extrapolation222The
uncertainties on excitation energies do not originate from this extrapolation
but are taken from the difference between the results obtained for the largest
$N_{\mathrm{max}}=8$ and the smallest space. Excitation energies are more
accurate than absolute ones because they converge faster with
$N_{\mathrm{max}}$..

Additional many-body methods are also considered for comparisons. First, the
sub-cases of PGCM and PGCM-PT obtained by only using one ”seed” HFB state,
i.e. omitting the GCM part of the calculation, are considered and referred to
as PHFB and PHFB-PT methods. The case where the projection part is further
omitted is utilized as well. This single-reference limit of PGCM-PT has been
formally elaborated on in Paper I and denotes a symmetry-breaking scheme in
case the seed state (i.e. becoming the unperturbed state) is itself symmetry
breaking. This limit will be also compared to the standard single-reference
symmetry-breaking Bogoliubov many-body perturbation theory (BMBPT)
Duguet:2015yle ; Tichai18BMBPT ; Arthuis:2018yoo ; Demol:2020mzd ;
Tichai2020review .

Our study focuses on three nuclei of increasing complexity. In each case, a
different collective coordinate $q$ is employed at the constrained HFB (cHFB)
level and for the subsequent GCM mixing333In each case, the employed interval
of $q$ values ensures the convergence of the PGCM calculation with respect to
that degree of freedom.. The characteristics of the associated mean-field,
PGCM and PGCM-PT calculations are

  1. 1.

Doubly closed-shell ${}^{16}\mathrm{O}$

     * •

spherically-symmetric Hartree-Fock (HF) states,

     * •

constraint on the root-mean-square (rms) matter radius ($q\equiv
r_{\text{rms}}$),

     * •

no symmetry projection needed.

  2. 2.

Singly open-shell ${}^{18}\mathrm{O}$

     * •

spherically-symmetric HFB states,

     * •

constraint on the pairing gap ($q\equiv\delta$) Duguet:2020hdm ,

     * •

projection of neutron number $N$.

  3. 3.

Doubly open-shell ${}^{20}\mathrm{Ne}$

     * •

axially-deformed HFB states,

     * •

constraint on the axial quadrupole moment ($q\equiv q_{20}$),

     * •

projections on neutron $N$ and proton $Z$ numbers as well as on angular
momentum $J$.

####  2.1.1 ${}^{16}\mathrm{O}$

Figure 2: (Color online) Ground-state energy of ${}^{16}\mathrm{O}$ as a
function of $r_{\text{rms}}$ of the (underlying) HF vacua.

In doubly closed-shell systems, the mean-field solution is nothing but a
spherical HF state. Since the rms radius operator employed to perform
constrained calculations commutes with the total angular momentum $J^{2}$, all
mean-field states involved in the ${}^{16}\mathrm{O}$ calculation carry good
symmetry quantum numbers and no symmetry projection is necessary in the
subsequent PGCM and PGCM-PT calculations, hence they are simply referred to as
GCM and GCM-PT, respectively.

The ground-state total energy curves (TECs) of ${}^{16}\mathrm{O}$ are
displayed in Fig. 2 as a function of the rms radius $r_{\text{rms}}$ of the
(underlying) HF vacua. One first observes that cHF and GCM results are
underbound by about 20 MeV ($12\%$) with respect to FCI, missing significant
IR dynamical correlations. In the present case444The GCM and GCM-PT(2)
calculations are performed on the basis of the nine cHF states visible on the
TEC., the GCM adds almost no energy (specifically, $165$ keV) to the HF
minimum, which signals that static IR collective correlations are marginal in
such a doubly closed-shell nucleus.

Given the negligible character of static correlations, ${}^{16}\mathrm{O}$
acts as a good benchmark for the (P)GCM-PT formalism. First, its single-
reference reduction HF-PT(2) is, as formally demonstrated in Paper I,
identical to canonical MBPT(2), i.e. Møller-Plesset MBPT based on the
unconstrained HF solution at the minimum of the TEC ($r_{\text{rms}}=2.03$
fm). While both single-reference partitionings of the Hamiltonian provide
slightly different results away from the minimum of the HF TEC, they are
qualitatively and quantitatively similar. The minima of the two TECs are close
to the FCI result. However, perturbation theories are not variational such
that it is difficult to argue that these values are to be preferred to
canonical ones. As a matter of fact MBPT(3) (not shown) does not flatten the
curve in the vicinity of the lowest MBPT(2) value555Empirically, the MBPT
expansion shows less sign of convergence away from the canonical point such
that the corresponding MBPT(2) values should not be preferred..

Figure 3: (Color online) Upper panel: collective ground-state GCM wave-
function probability distribution
($|\breve{f}^{0^{+}_{1}}(r_{\text{rms}})|^{2}$) in ${}^{16}\mathrm{O}$ as a
function of the rms radius of the underlying HF vacua. Lower panel:
contributions to GCM ($e^{(0+1)}_{0}(r_{\text{rms}})$) and GCM-PT(2)
($e^{(2)}_{S}(r_{\text{rms}})+e^{(2)}_{D}(r_{\text{rms}})$) ground-state
energies as a function of $r_{\text{rms}}$. GCM-PT(2) constributions are split
into single (one-particle/one-hole) and double (two-particle/two-hole)
excitations.

Focusing on the canonical point, one observes that GCM-PT(2) is consistent
with MBPT(2)/HF-PT(2), adding only $146$ keV static correlation energy. This
consistency constitutes a validation of GCM-PT(2), knowing that it is formally
very different from MBPT(2) and relies on a completely different numerical
procedure as can be appreciated from the various appendices to the present
paper.

Furthermore, this consistency sheds some light on single-reference MBPT(2)/HF-
PT(2) results. The upper panel of Fig. 3 shows that, while the GCM ground-
state collective wave-function spreads over a large interval of
$r_{\text{rms}}$ values due to nuclear-size fluctuations, the Hamiltonian
dictates that the contributions to the left of the HF minimum (i.e. for
$r_{\text{rms}}\leq 2.03$ fm) dominate it. From the energetic viewpoint, the
lower panel of Fig. 3, which shows the decomposition666The PGCM collective
wave function and the contribution $e^{(0+1)}_{0}(r_{\text{rms}})$ of each
value of the collective coordinate to the PGCM energy are introduced in App. A
of Paper II. of the GCM energy as a function of $r_{\text{rms}}$ demonstrates
that the largest contributions originate from configurations centered around
the HF minimum. Next777The decomposition of the PGCM-PT(2) correlation energy
is provided in Sec. 3.3.2 of Paper I., the lower panel also illustrates that
the physically-informed weights in the GCM unperturbed state propagate to GCM-
PT(2) such that configurations around the HF minimum contribute the most to
the second-order correction whereas those associated with the lowest
MBPT(2)/HF-PT(2) values around $r_{\text{rms}}\in[2.1,2.2]$ fm are largely
subleading. Eventually, the total GCM-PT(2) energy is nearly identical to
canonical MBPT(2)/HF-PT(2) results. This definitely gives more credit to low-
order MBPT(2)/HF-PT(2) energies obtained at the canonical point than to those
obtained at smaller and larger values of $r_{\text{rms}}$. Interestingly, one
also observes that the GCM-PT(2) energy correction is dominated by double
(two-particle/two-hole) excitations given that the energy contribution of
single (one-particle/one-hole) excitations is negligible at all values of
$r_{\text{rms}}$. While this feature is expected at the canonical point given
that single excitations do not contribute to MBPT(2)/HF-PT(2)888This a
consequence of Brillouin’s theorem that implies a decoupling between the HF
reference state and its singles excitations at the canonical point., it is not
evident away from it.

Eventually, the GCM-PT(2) binding energy differs by $0.8\%$ from the FCI
result. A common theme throughout the paper regards the best way to achieve
even greater accuracy. At this point, one can either hope to enrich the PGCM
unperturbed state by selecting a potentially pertinent additional collective
coordinate $q$ and/or go to PGCM-PT(3)999The benefit of going to PGCM-PT(3)
(see Ref. Ladoczki20  for a similar situation) given the associated numerical
scaling makes probably more efficient to seek for a further improvement of the
PGCM unperturbed state.. A third (complementary) option to achieve such a goal
will be introduced in Sec. 3.

####  2.1.2 ${}^{18}\mathrm{O}$

Figure 4: Ground-state energy of ${}^{18}\mathrm{O}$ as a function of the
pairing constraint $\delta$ calculated from various many-body methods.

The singly open-shell ${}^{18}\mathrm{O}$ constitutes the first nucleus of the
present study in which static correlations are expected to be significant. In
this particular case, static correlations relate to superfluidity and thus
translate first at the HFB level into the spontaneous breaking of the $U(1)$
global-gauge symmetry associated with particle number conservation.
Correspondingly, the pairing gap operator is used here as a constraint to vary
the amount of pairing correlations in the HFB seeds Duguet:2020hdm . As a next
step, further static correlations are captured via the restoration of neutron
number and the inclusion of pairing fluctuations through the PGCM.

Figure 5: (Color online) Upper panel: collective ground-state PGCM wave-
function probability distribution ($|\breve{f}^{0^{+}_{1}}(\delta)|^{2}$) in
${}^{18}\mathrm{O}$ as a function of the pairing constraint $\delta$
characterizing the underlying HFB vacua. Lower panel: contributions to the
PGCM ($e^{(0+1)}_{0}(\delta)$) ground-state energy and to the PGCM-PT(2)
($e^{(2)}_{S}(\delta)+e^{(2)}_{D}(\delta)$) correlation energy. The latter is
split into single (two quasi-particle) and double (four quasi-particle)
contributions.

The ground-state TECs of ${}^{18}\mathrm{O}$ are displayed in Fig. 4 as a
function of the pairing constraint $\delta$ of the (underlying) HFB
vacua101010See Ref. Duguet:2020hdm  for the definition of $\delta$.. By
definition $\delta=1$ corresponds to the canonical, i.e. unconstrained, HFB
solution. While the PHFB TEC follows the HFB one, it is less bound, e.g. by
$1.2$ MeV at the canonical point. The fact that the particle-number projection
after variation (PNPAV) decreases the binding reflects the fact that the
distribution of particle numbers in the HFB state around the average is
distorted towards heavier systems. In the next step, the GCM mixing associated
with the inclusion of pairing fluctuations yields negligible correlation
energy compared to the PNPAV that provides the essential IR correlations.

Similarly to ${}^{16}\mathrm{O}$, PGCM underbinds the FCI result by about $25$
MeV ($\sim 13\%$), thus missing significant IR dynamical correlations. While
formally not identical to canonical BMBPT(2)111111See the related discussion
in Paper I., the single-reference reduction of PGCM-PT(2), here denoted by
HFB-PT(2), captures dynamical correlations on top of HFB. Results from both
single-reference methods are very similar and agree with FCI within
uncertainties.

However, this close agreement is accidental and somewhat spurious. Indeed,
PHFB-PT(2), which actually corrects for the $U(1)$ breaking of HFB-PT(2),
pushes the energy up by about $1.5$ MeV away from the FCI result at the
canonical point. This number is close to the difference between HFB and PHFB
mentioned above. While the present calculation constitutes the first example
demonstrating the impact of exactly restoring symmetries within (perturbative)
expansion methods, it is seen below that the $1.5$ MeV shift shall not be
fully attributed to the symmetry restoration. Adding the GCM mixing into the
unperturbed state, the PGCM-PT(2) result remains consistent with PHFB-PT(2) at
the canonical point within uncertainties.

Going away from the canonical point, BMBPT(2) and HFB-PT(2) differ. This
behavior reflects the different nature of the partitionings used by the two
expansions, which is magnified as one departs from the canonical point. At the
same time, PHFB-PT(2) becomes less (more) bound than PGCM-PT(2) as $\delta$
becomes smaller (greater) than $1$. Once again, these behaviors do not instill
trust in perturbative results away from the canonical point. Thankfully, PGCM-
PT(2) is better controlled given that the configurations associated with
different values of the collective coordinate $\delta$ enter the unperturbed
PGCM state with weights dictated by the physical Hamiltonian. As shown in Fig.
5, the collective PGCM wave-function spreads significantly on both sides of
the canonical point with a maximum located to the left of it $(\delta=0.7)$.
While the decomposition of the PGCM energy reflects this distribution, the
second-order correction is flatter with $\delta$ but slightly favors values
smaller than 1. Eventually, the PGCM-PT(2) binding energy is very close to the
PHFB-PT(2) energy at the canonical point and lies $1.5$ MeV (0.8%) above the
FCI result.

Still, PGCM-PT(2) and PHFB-PT(2) results carry error bars associated with the
approximate solution of the linear system at play in the formalism. Due to to
its large dimension, this linear system is solved iteratively as discussed in
App. B, introducing an uncertainty that can be evaluated through Eq. (24). The
solution can also be affected by linear redundancies and intruder problems
that are dealt with via the simultaneous use of a norm preconditioning and a
complex shift $\gamma$ as detailed in Apps. B.3.1 and C, respectively. While
increasing the precision, the use of an overly large complex shift may degrade
the accuracy by generating a bias in the extracted value.

In ${}^{18}\mathrm{O}$, which qualifies as a difficult case, the iterative
procedure can be converged in a stable fashion with a complex shift
$\gamma=10$ MeV, eventually leading to a $\pm 0.3$ MeV precision on the PGCM-
PT(2) energy121212The precision on the PHFB-PT(2) is better ($\pm 0.1$ MeV)
thanks to the lower dimension and the near diagonal character of the linear
system. that is visualized by a band in Fig. 4. While the central value
reported in Fig. 4 is obtained for $\gamma=10$ MeV, the bias (not reported on
the figure) due to this complex shift131313The bias is estimated by varying
the shift over the interval $\gamma\in[5,15]$ MeV, see App. C.3 for an
illustration. pushes the PGCM-PT(2) energy up by about $1$ MeV. Eventually,
the bias accounts for two thirds of the $1.5$ MeV (0.8%) disagreement with the
FCI result and for two thirds of the shift upward compared to HFB-PT(2) that
was fully attributed to the symmetry restoration at first.

Figure 6: Ground-state energy of ${}^{20}\mathrm{Ne}$ as a function of the
axial quadrupole deformation $\beta_{2}$ of the (underlying) HFB states.

####  2.1.3 ${}^{20}\mathrm{Ne}$

The doubly open-shell ${}^{20}\mathrm{Ne}$ displays strong static correlations
that manifest through the breaking of $SU(2)$ rotational symmetry associated
with angular momentum conservation at the HFB level. Accordingly, the axial
quadrupole moment operator is used as a constraint to vary the deformation of
the HFB seeds. As a next step, further static correlations are captured via
the restoration of angular momentum and the inclusion of shape fluctuations
through the PGCM. As demonstrated in Paper II, the description of
${}^{20}\mathrm{Ne}$ strongly benefits from breaking and restoring parity as
well as the inclusion of octupole shape fluctuations. Our present calculations
are however restricted to axial quadrupole deformation, leaving some room for
further improvement in the future. While $U(1)$ global gauge symmetry is also
allowed to break spontaneously, it does not do so with the presently employed
Hamiltonian, hence all HFB states actually reduce to (deformed) HF Slater
determinants.

Figure 7: (Color online) Upper panel: collective ground-state PGCM wave-
function probability distribution ($|\breve{f}^{0^{+}_{1}}(\beta_{2})|^{2}$)
in ${}^{20}\mathrm{Ne}$ as a function of the axial quadrupole deformation
$(\beta_{2})$ of the underlying HFB vacua. Lower panel: contributions to PGCM
($e^{(0+1)}_{0}(\beta_{2})$) and PGCM-PT(2)
($e^{(2)}_{S}(\beta_{2})+e^{(2)}_{D}(\beta_{2})$) ground-state energies as a
function of the axial quadrupole deformation $(\beta_{2})$ of the underlying
HFB vacua. The PGCM-PT(2) constributions are split into singles (two quasi-
particle) and doubles (four quasi-particle) contributions. Figure 8: (Color
online) Excitation energy in ${}^{20}\mathrm{Ne}$ as a function of the axial
quadrupole deformation $(\beta_{2})$ of the underlying HFB vacua. Top panel:
first $4^{+}$ state. Bottom panel: first $2^{+}$ state. Calculations are
performed with $\hbar\omega=20\text{ MeV}$, $e_{\text{max}}=4$ and employing
the two-body part of the N3LO $\chi$EFT Hamiltonian evolved to
$\lambda_{\text{vsrg}}=1.88$ fm-1.

The ground-state TECs of ${}^{20}\mathrm{Ne}$ are displayed in Fig. 6 as a
function of the axial quadrupole deformation141414See Paper II for a precise
definition. $\beta_{2}$ of the (underlying) HFB vacua. One first observes that
the projection on $J$ provides a significant energy gain of $5.5$ MeV and
moves the minimum of the PHFB TEC to larger deformation ($\beta_{2}=0.35$)
than the canonical HFB minimum ($\beta_{2}=0.3$). The GCM mixing only adds
$80$ keV correlation energy given that the TEC is rather stiff along the axial
quadrupole direction151515As shown in Paper II, the energy is softer against
axial octupole deformations.. Once again, static correlations are dominated by
the symmetry restoration. Having included essential static correlations, the
PGCM energy is still $21.7$ MeV ($10\%$) away from the FCI result, and misses
significant dynamical correlations.

Stepping back to canonical HFB and adding dynamical correlations via BMBPT(2)
lowers the energy by $24.6$ MeV, yielding a result that is $2.6$ MeV ($1.2\%$)
underbound compared to FCI161616Canonical BMBPT(2) is the closest point to FCI
along the TEC in the present example. Note that canonical BMBPT(3) only
provides an extra $0.3$ MeV correlation energy compared to canonical
BMBPT(2)..

On the other hand, starting from the PHFB TEC and adding dynamical
correlations via PHFB-PT(2) lowers the energy by $25.1,24.9$ and $25.7$ MeV at
the HFB, PHFB, and PHFB-PT(2) minima, respectively. These energies overshoot
the FCI result by about $2.5/3.2/3.4$ MeV ($1.2/1.5/1.6\%$). While the
difference between BMBPT(2) and PHFB-PT(2) TECs is similar to the difference
between HFB and PHFB TECs, one observes that a consistent angular-momentum
restoration favors larger deformations when adding dynamical correlations.

The mixing of quadrupole shapes in PGCM-PT(2) only adds $310$ keV to the PHFB-
PT(2) minimum. The PGCM-PT(2) result keeps a close memory of the PHFB-PT(2)
minimum ($\beta_{2}=0.4$) rather than the PHFB-PT(2) value at the canonical
HFB minimum ($\beta_{2}=0.3$). All in all, the PGCM-PT(2) energy171717Present
PGCM-PT(2) and PHFB-PT(2) results were obtained with a complex shift
$\gamma=15$ MeV. The precision error associated with solving the linear system
is shown through an error band in Fig. 23. overshoots the FCI result by
$1.7\%$. This discrepancy is expected to decrease after the inclusion of the
octupole degree of freedom into the PGCM.

Figure 9: Absolute energies of the first $0^{+}$, $2^{+}$ and $4^{+}$ states
in ${}^{20}\mathrm{Ne}$ computed via PGCM, PGCM-PT(2) and FCI.

In order to further analyse the theoretical content of the above results, Fig.
7 shows that the collective PGCM ground-state wave-function and the associated
energy contributions are distributed rather symmetrically around the
$J^{\pi}=0^{+}$ PHFB minimum ($\beta_{2}=0.35$) of the TEC visible in Fig. 6
and spread over a large interval of $\beta_{2}$ values. Interestingly,
dynamical correlations captured via PGCM-PT(2) favor configurations181818Once
again, single excitations bring negligible contributions to the correlation
energy. to the left of the HFB minimum ($\beta_{2}=[0.25,0.30]$). As a result,
dynamical correlations could counterbalance the overestimated radii obtained
at the PGCM level (see Paper II) due to the opposite predilection of the
latter for deformations larger than the HFB minimum. This interesting and non-
trivial finding will have to be confirmed by an explicit calculation of rms
radii at the PGCM-PT(2) level in the future.

In addition to providing accurate absolute energies in complex systems, e.g.
in doubly open-shell nuclei displaying strong collective static correlations,
a key advantage of the multi-reference PGCM-PT formalism over BMBPT is that it
provides natural access to the low-lying spectroscopy within a symmetry-
conserving scheme by correcting each PGCM eigenstate for dynamical
correlations.

The first $2^{+}$ and $4^{+}$ excitation energies in ${}^{20}\mathrm{Ne}$ are
shown in Fig. 8 as a function of the axial quadrupole deformation. First, one
observes that the PGCM $2^{+}_{1}$ and $4^{+}_{1}$ excitation energies differ
from the FCI results by $300$ keV ($27\%$) and $560$ keV ($13\%$),
respectively. This is consistent with the results displayed in Paper II. One
also sees that PHFB results at the canonical deformation ($\beta_{2}=0.3$) are
very close to PGCM ones, but the differences grow for smaller or large
deformations. Adding dynamical correlations, PHFB-PT(2) flattens the
excitation energies as a function of $\beta_{2}$ compared to PHFB,
systematically going into the direction of PGCM-PT(2) for each deformation.
Given that exact results would be independent of the deformation of the
underlying vacuum, this feature is an empirical sign that PHFB-PT(2) results
are better converged than PHFB ones. It also implies that the PGCM-PT(2)
spectrum converges with fewer states than the PGCM one. Still, at the
canonical deformation ($\beta_{2}=0.3$) dynamical correlations are small,
which remains true even when shape mixing is added, given that PGCM-PT(2)
excitation energies are essentially identical to PGCM ones.

Overall, the PGCM-PT(2) $2^{+}_{1}$ and $4^{+}_{1}$ excitation energies differ
by $24\%$ and $15\%$ from FCI results respectively, which seems to indicate
that missing correlations are beyond two-particle/two-hole excitations of
axially deformed HF states. While going to PGCM-PT(3) will help reduce this
difference, it might be numerically less costly and more relevant in this case
to enrich the PGCM unperturbed state via, e.g., the inclusion of octupole,
triaxial and/or pairing degrees of freedom, or to start from HFB states
obtained via a variation after particle-number-projection (VAPNP) calculation,
in order to compress the spectrum. In the future, another possibility would be
to design a non-perturbative extension of the multi-reference PGCM-PT
formalism to more efficiently capture higher-rank particle-hole excitations.

Our ${}^{20}\mathrm{Ne}$ results are summarized in Fig. 9 where the combined
benefits of PGCM-PT are clearly apparent. Although a slight overbinding of
about $3$ MeV ($\sim 1.5\%$) is observed, PGCM-PT(2) brings down absolute
energies to the right range of values without degrading their relative
position. This latter feature is far from trivial given that the PGCM-PT
formalism is state specific, i.e. calculations are performed separately on top
of each PGCM eigenstate, and considering that each PGCM energy is corrected by
about $25$ MeV while their relative distance is on the MeV scale. In
particular, the (non-trivial) numerical techniques used to solve the PGCM-
PT(2) equations must be well controlled to maintain the consistency of the
spectra. For example, it is essential to use the same complex shift $\gamma$
for all states belonging to a given nucleus in order for the bias on absolute
energies to be consistent and to largely cancel out in the excitation
spectrum.

##  3 Adding the MR-IMSRG pre-processing

In the present part, PGCM-PT(2) calculations are performed in a larger model
space with $e_{\text{max}}=6$ (and $\hbar\omega=16$ MeV)191919The limitation
to $e_{\text{max}}=4$ was due to the wish to benchmark PGCM-PT(2) calculations
against FCI results.. We use a Hamiltonian consisting of an SRG-evolved chiral
N3LO nucleon-nucleon interaction with $\lambda_{\text{vsrg}}=1.8\text{
fm}^{-1}$, supplemented with an N2LO three-nucleon interaction with cutoff
$\Lambda=2.0\,\mathrm{fm}^{-1}$ whose low-energy constants are adjusted to
$A=3,4$ observables, as described in Refs. PhysRevC.83.031301 ; Nogga:2004il .
The Hamiltonian is further pre-processed via the MR-IMSRG unitary
transformation based on the $J^{\pi}=0^{+}$ canonical PHFB state. The
evolutions are based on the MR-IMSRG(2) truncation scheme, employing the so-
called Brillouin generator — see Refs. Hergert:2015awm ; Hergert:2016etg  for
details. The MR-IMSRG transformation is parametrized by the flow parameter
$s\in[0,20]\,\mathrm{MeV}^{-1}$, where $s=0$ means that no transformation is
applied and the upper limit is chosen such that the transformed Hamiltonian no
longer exhibits significant evolution. In the MR-IMSRG(2) truncation scheme,
the three-nucleon interaction included at the beginning of the flow (along
with higher-body operators for $s>0$) is approximated via the normal ordering
with respect to the J=0 PHFB state.

In closed-shell nuclei (not shown here), the PHFB reference state reduces to a
spherically invariant Slater determinant such that MR-IMSRG is nothing but the
simpler SR-IMSRG method. In this case, pushing the transformation to
$s=\infty$ ($s$ sufficiently large in practice) leads to a complete
resummation of dynamical correlations into the pre-processed Hamiltonian such
that the unperturbed HF Slater determinant becomes its exact ground state,
i.e. no further correlations need to be added. While dynamical correlations
are largely resummed in open-shell nuclei via the MR-IMSRG pre-processing, the
decoupling of the reference state cannot be complete Hergert:2016etg  such
that an additional step is always needed to grasp the remaining correlations
as illustrated below202020While the exact decoupling is formally realized in
the simpler SR-IMSRG method, it can usually be achieved to good accuracy in
actual MR-IMSRG calculations when pushing the transformation far enough.
Strictly speaking, however, the decoupling cannot be complete in MR-IMSRG, at
least for generators that have an explicit and manageable second-quantized
representation. One may be able to define a generator that formally achieves
the decoupling but such a generator would have no practical, i.e. low-rank,
representation..

###  3.1 ${}^{18}\mathrm{O}$

Figure 10: (Color online) Absolute binding energy of ${}^{18}\mathrm{O}$ as a
function of the flow parameter $s$ associated with the MR-IMSRG pre-processing
of the Hamiltonian.

The absolute binding energy of ${}^{18}\mathrm{O}$ is displayed in Fig. 10 as
a function of $s$. Due to the PNPAV in ${}^{18}\mathrm{O}$, the number of
single and double excitations of the HFB vacuum required to perform a PHFB-
PT(2) calculation is already very large ($10^{6}$ states) for
$e_{\text{max}}=6$. The numerical implementation will be optimized in the
future, but in the mean-time, the calculation is made faster by discarding
configurations based on their norm as specified in App. D. For the same
reason, only PHFB-PT(2) calculations on top of the spherical
${}^{18}\mathrm{O}$ canonical HFB vacuum have been performed, leaving a PGCM-
PT(2) calculation for the future.

Figure 11: (Color online) Correlation energy, i.e. difference to the canonical
HFB result, in ${}^{18}\mathrm{O}$ as a function of the flow parameter $s$
associated with the MR-IMSRG pre-processing of the Hamiltonian.

In spite of the change of model space and Hamiltonian, the situation
encountered at $s=0$ is qualitatively similar to the one discussed in Sec.
2.1. Indeed, while HFB and PHFB are largely underbound, BMBPT(2) and PHFB-
PT(2) bring in the dominant fraction of dynamical correlations212121Contrarily
to the results obtained in Sec. 2.1.2 with a two-body interaction only and
$e_{\text{max}}=4$, PHFB-PT(2) is very close to BMBPT(2) at $s=0$. At the same
time, the contribution of BMBPT(3) is enlarged., with BMBPT(3) adding an extra
$2$ MeV. Switching on the MR-IMSRG pre-processing, HFB and PHFB energies drop
dramatically for small values of $s$ and flatten out very quickly beyond
$s=1\,\mathrm{MeV}^{-1}$. At the same time, BMBPT(2), PHFB-PT(2) and BMBPT(3)
drop towards a similar value, about $1.5$ MeV below the original BMBPT(3)
result, which happens to be also similar to the PHFB value. Eventually, PHFB-
PT(2) is about $4$ MeV ($2.9\%$) away from experiment. No convergence analysis
as a function of the model space has been performed and reaching a converged
absolute binding energy clearly requires (an extrapolation to) a larger model
space.

To better appreciate the impact of the MR-IMSRG evolution, the correlation
energy, i.e. the difference to the HFB result, is shown in Fig. 11. Having
already absorbed the bulk of dynamical correlations, pre-processed
Hamiltonians become more and more perturbative with increasing $s$ such that
BMBPT(2,3) and PHFB-PT(2) corrections become less important with the flow,
i.e. one goes from $38.8$ MeV and $36.2$ MeV for BMBPT(3) and PHFB-PT(2) at
$s=0$ to $288$ keV and $369$ keV at $s=10\,\mathrm{MeV}^{-1}$, respectively,
with an inversion of the two results. At the same time, the particle number
projection that is repulsive at $s=0$ ($-394$ keV) brings in additional
binding for $s\geq 1\,\mathrm{MeV}^{-1}$ ($+327$ keV at
$s=10\,\mathrm{MeV}^{-1}$). These results demonstrate that correlations are
reshuffled through the MR-IMSRG flow, such that the importance of dynamical
correlations is strongly reduced whereas static correlations are somewhat
enhanced.

Dynamical correlations added on top of PHFB via PHFB-PT(2)222222The numerical
solution of the PHFB-PT(2) linear system is very stable in the present example
such that a small complex shift ($\gamma=1$ MeV) can be used safely. The
precision error on PHFB-PT(2) energies is essentially invisible in Fig. 11
whereas the bias generated for $\gamma=1$ MeV is negligible compared to the
$42$ keV difference between PHFB and PHFB-PT(2) energies at
$s=10\,\mathrm{MeV}^{-1}$. become as small as $42$ keV at
$s=10\,\mathrm{MeV}^{-1}$. Thus, the PHFB state used as a reference for the
MR-IMSRG pre-processing is, for all practical purposes, decoupled from the
${\cal Q}$ space at the end of the transformation in the present calculation.
Although the decoupling cannot be exact in principle, ${}^{18}\mathrm{O}$
behaves similarly to a closed-shell nucleus such that the dynamical
correlations left to be captured after PNPAV are very small.

Figure 12: (Color online) $J^{\pi}=0^{+},2^{+},4^{+}$ PHFB TECs in
${}^{20}\mathrm{Ne}$ as a function of the axial quadrupole deformation
$\beta_{2}$ for $s=0\,\mathrm{MeV}^{-1}$ (upper panel),
$s=10\,\mathrm{MeV}^{-1}$ (middle panel) and $s=20\,\mathrm{MeV}^{-1}$ (lower
panel). Figure 13: (Color online) Absolute PGCM and PGCM-PT(2) binding
energies of ${}^{20}\mathrm{Ne}$ as a function of the MR-IMSRG flow parameter
$s$.

###  3.2 ${}^{20}\mathrm{Ne}$

The doubly open-shell ${}^{20}\mathrm{Ne}$ constitutes a richer and more
instructive example. Figure 12 shows the $J^{\pi}=0^{+},2^{+}$ and $4^{+}$
PHFB TECs as a function of the axial quadrupole deformation $\beta_{2}$ for
three values $(s=0,10,20)\,\mathrm{MeV}^{-1}$ of the MR-IMSRG flow parameter
Hergert:2016etg . The TECs are strongly lowered with $s$, e.g. the PHFB
minimum gains $45.4$ MeV going from $s=0$ to $s=20\,\mathrm{MeV}^{-1}$, with
most of the effect occuring for $0\leq s\leq 10\,\mathrm{MeV}^{-1}$. At the
same time, the deformation of the PHFB minimum is lowered from
$\beta_{2}=0.55$ to $\beta_{2}=0.52$ while the TECs become stiffer.

In Fig. 13, PGCM and PGCM-PT(2) binding energies are displayed as a function
of the flow parameter. Starting from $J^{\pi}=0^{+}$ PHFB TECs, PGCM and PGCM-
PT(2) calculations mix five HFB configurations with axial quadrupole
deformations $\beta_{2}=(0.3,0.4,0.5,0.6,0.7)$. Unlike in ${}^{18}\mathrm{O}$,
the convergence of PGCM energies is not fully reached yet for
$s=20\,\mathrm{MeV}^{-1}$. Still, the bulk of dynamical correlations has
already been resummed into the pre-processed Hamiltonian at $s=10$, which
suggests a convergent behavior. Eventually, the PGCM energy is lowered by
$45.2$ MeV between $s=0$ and $s=20\,\mathrm{MeV}^{-1}$. At the same time,
PGCM-PT(2) systematically lowers the PGCM value, the added dynamical
correlations reducing from $42.5$ MeV at $s=0$ to only $2.0$ MeV at
$s=20\,\mathrm{MeV}^{-1}$. Similarly, the difference between PHFB and PGCM-
PT(2) is drastically reduced as $s$ grows but does not vanish, i.e. it still
amounts to $2.03$ MeV with the most pre-processed Hamiltonian232323PHFB and
PGCM energies differ by less than $200$ keV all throughout the interval
$s\in[0,20]\,\mathrm{MeV}^{-1}$.. This indicates that, while very effective,
the decoupling of the PHFB state from the ${\cal Q}$ space is not complete and
thus less effective than in the singly open-shell ${}^{18}\mathrm{O}$. This
feature points to the stronger multi-reference character of
${}^{20}\mathrm{Ne}$ associated with the breaking and restoration of $SU(2)$
symmetry.

Figure 14: (Color online) $2_{1}^{+}$ (bottom row) and $4_{1}^{+}$ (top row)
excitation energies as a function of $\beta_{2}$ for $s=0$ (left column),
$s=10\,\mathrm{MeV}^{-1}$ (middle column) and $s=20\,\mathrm{MeV}^{-1}$ (right
column).

The PGCM-PT(2) energy changes by less than $5$ MeV over the interval
$s\in[0,20]\,\mathrm{MeV}^{-1}$, thus strongly reducing the flow parameter
dependence compared to PGCM results. The residual dependence of the ground-
state energy on the flow parameter results both from the breaking of unitarity
associated with the truncation of the flow equations at the MR-IMSRG(2) level
and from the approximations to the solution of the $A$-body Schrödinger’s
equation at the PGCM-PT(2) level. Under the hypothesis that the PGCM-PT is
convergent and given that the second-order correction reduces to $2$ MeV at
$s=20\,\mathrm{MeV}^{-1}$, one can speculate that the PGCM-PT(2) energy is
eventually better converged than the $5$ MeV spread over the interval
$s\in[0,20]\,\mathrm{MeV}^{-1}$, i.e. by better than $3\%$. The fact that the
experimental value is consistent with the PGCM-PT(2) prediction within
estimated uncertainty must not be overinterpreted given that improving over
the presently used $e_{\text{max}}$ truncation is expected to lower the
ground-state energy by several MeVs.

Figure 15: (Color online) Low lying spectrum of ${}^{20}\mathrm{Ne}$ as a
function of the MR-IMSRG flow parameter.

Turning to the low-lying spectroscopy, Fig. 14 displays the first $2^{+}$ and
$4^{+}$ excitation energies as a function of $\beta_{2}$ for the three values
of the flow parameter. Focusing first on $s=0$, the conclusions drawn in Sec.
2.1.3 remain valid, i.e. PHFB-PT(2) flattens the excitation energies as a
function of $\beta_{2}$ compared to PHFB whereas dynamical correlations
brought in through PGCM-PT(2) do not modify the low-lying part of the PGCM
ground-state rotational band. However, the picture changes drastically when
pre-processing the Hamiltonian via MR-IMSRG. Indeed, the PGCM spectrum becomes
more dilated with increasing $s$. This feature is already visible from the
$J^{\pi}=0^{+},2^{+},4^{+}$ PHFB TECs displayed in Fig. 12 that become more
distant with increasing $s$. Based on this trend, one observes that PHFB-
PT(2), while always flattening the dependence on $\beta_{2}$, systematically
corrects for this dilatation of the rotational spectrum. This non-trivial
feature is confirmed at the PGCM-PT(2) level.

Figure 16: (Color online) Total ground-state energy of 20Ne computed within
various many-body methods for three different values of the MR-IMSRG flow
parameter. Numbers next to downward arrows denote the corresponding gain in
correlation energy (in MeV).

This key result can be better appreciated in Fig. 15 where PGCM and PGCM-PT(2)
spectra are compared to experiment and to the spectrum obtained from a richer
PGCM calculation including additional axial states along with triaxially
deformed ones. Although the PGCM calculation based on five axial states is
rudimentary, the observed dilatation of spectra is not compensated for by such
an enrichment of the PGCM unperturbed state. Correspondingly, the systematic
compensation of that dilatation via PGCM-PT(2) corresponds to a genuine action
of the perturbation that captures dynamical correlations lying outside the
reach of the presently used PGCM ansatz. In the end, the PGCM-PT(2)
$2^{+}_{1}$ excitation energy is independent of $s$ within uncertainties.
While reduced compared to PGCM, the $s$ dependence of the PGCM-PT(2)
$4^{+}_{1}$ excitation energy is still significant and would probably benefit
from being performed on top of a richer PGCM state and/or by going to PGCM-
PT(3).

The global picture that emerges for pre-processed Hamiltonians is illustrated
for ${}^{20}\mathrm{Ne}$ in Fig. 16. The MR-IMSRG evolution largely reshuffles
the hierarchy of correlations at play. As $s$ grows, one observes that

  1. 1.

static correlations captured through the breaking of symmetries at the HFB
level as well as by their restoration and the inclusion of collective
fluctuations at the PGCM level slightly increase,

  2. 2.

dynamical correlations brought either on top of HFB via BMBPT(2) or on top of
PGCM via PGCM-PT(2) are drastically reduced.

Overall, dynamical correlations go from being highly dominant to being largely
subleading. Still, their inclusion on top of PGCM via PGCM-PT(2) remains
mandatory, in particular when dealing with low-lying excitation energies.

Eventually, the great benefit of the pre-processing relates to the fact that
many-body calculations performed with evolved Hamiltonians become numerically
gentler as $s$ increases, i.e. the numerical solution of the PGCM-PT(2) linear
system is more precise, corrections beyond PGCM-PT(2) are minimized and the
convergence with the model-space size ($e_{\mathrm{max}}$) is probably faster,
although this latter point remains to be studied242424See Ref.
Gebrerufael:2017fk  for an accelerated convergence in so-called in-medium no-
core shell model calculations.. Given that PGCM-PT(2) is numerically more
costly than the MR-IMSRG(2) step (see App. D.1), the optimal combination of
both methods is of great interest. Of course, this optimal point must be such
that the error due to the breaking of unitarity through the MR-IMSRG(2) pre-
processing is not larger than the error associated with PGCM-PT(2) results.

##  4 Conclusions

This work, the third paper of the series on PGCM-PT, presented the first
realistic results for the novel multi-reference perturbation theory built on
top of an unperturbed state generated through the projected generator
coordinate method. While the unperturbed state captures crucial static
correlations via the breaking and restoration of symmetries along with
collective fluctuations, the perturbative expansion brings in complementary
dynamical correlations in a consistent fashion within a symmetry conserving
scheme. Furthermore, being a state-specific multi-reference many-body
perturbation theory, PGCM-PT accesses ground and low-lying excited states on
an equal footing.

First, the novel many-body formalism was shown to be both versatile and
accurate by benchmarking proof-of-principle results for the doubly closed-
shell ${}^{16}\mathrm{O}$, singly open-shell ${}^{18}\mathrm{O}$ and doubly
open-shell ${}^{20}\mathrm{Ne}$ nuclei in a small ($e_{\mathrm{max}}=4$)
harmonic oscillator model space against full configuration interaction
results. Binding energies obtained at second order, i.e. through PGCM-PT(2),
were shown to be typically $0.5-1.5\%$ away from FCI results.

The second focus of the present paper was to demonstrate the benefit of
combining low-order PGCM-PT with a pre-processing of the Hamiltonian via
multi-reference in-medium similarity renormalization group transformations.
The rather low cost of MR-IMSRG(2) calculations makes it possible to
efficiently capture the bulk of dynamical correlations in large model spaces
(cf. Refs. Gebrerufael:2017fk ; Yao:2018qjv ; PhysRevLett.124.232501 ). Based
on such a pre-processed Hamiltonian, PGCM-PT(2) can bring in crucial static
correlations and any remaining dynamical correlations while working in a
smaller model space. The present work showed that, after the MR-IMSRG(2) pre-
processing, dynamical correlations included on top of the PGCM via PGCM-PT(2)
are indeed essential for a satisfactory description of low-lying spectra.

Eventually, it emerges from the present work that a versatile and accurate
description of complex mid- and heavy-mass nuclei will probably rely on the
combination of three levers whose complementarity needs to be further studied
and optimized:

  1. 1.

the pre-processing of the Hamiltonian via, e.g., MR-IMSRG to efficiently
capture the bulk of dynamical correlations,

  2. 2.

the use of a, e.g., PGCM unperturbed state capturing collective static
correlations via a low-dimensional diagonalization problem that is thus
scalable to heavy nuclei252525This can be particularly useful for
implementations on GPUs or other accelerators with limited memory.,

  3. 3.

the low-order truncation of a systematic expansion on top of the multi-
reference unperturbed state via, e.g., PGCM-PT to bring in remaining dynamical
correlations.

Each of the three steps comes with its own flexibility that can be exploited
in order to optimize their combination262626It is worth mentioning that the
combination of the three steps is always consistent, i.e. there is no double
counting given that each step automatically adapts to the other two.. First,
the pre-processing is a function of a flow parameter $s$ that must be
optimized to resum the bulk of dynamical correlations without inducing a large
breaking of unitarity272727In this context, the truncation order of the MR-
IMSRG(n) procedure plays a critical role. MR-IMSRG(3), for instance, would
allow us to reduce any violation of the unitarity, but it comes with a
significantly higher numerical cost.. Second, the PGCM depends on a choice of
suitable collective coordinates that must be rich enough to capture all non-
perturbative static correlations at play, only leaving weak perturbative
corrections to the subsequent PGCM-PT step, while maintaining a low-enough
dimensionality to retain its advantage over large-scale diagonalization
methods. For example, while adding the triaxial degree of freedom did not
impact the dilated PGCM spectrum of 20Ne at $s=10,20\,\mathrm{MeV}^{-1}$, the
use of HFB states obtained while adding a cranking constraint breaking time-
reversal invariance PhysRevLett.113.162501 ; Borrajo:PLB2015  typically
compresses the PGCM spectrum as demonstrated in MR-EDF calculations and in
recent ab initio studies PhysRevLett.124.232501 ; Hergert:2020FP . Obtaining
such a compression at the PGCM level is expected to correlate with a further
suppression of dynamical correlations on top of the PGCM step. Still, if
needed, the PGCM-PT can in principle be implemented at various perturbative
orders $n$. In practice, however, going beyond PGCM-PT(2) shall probably be
avoided due to the prohibitive numerical scaling.

While the present work has laid the foundations of such an optimal scheme,
future studies will allow us to better understand the way many-body
correlations can be most efficiently captured in complex heavy nuclei within
an ab initio setting. For example, describing nuclei displaying strong shape
coexistence via ab initio many-body calculations constitute an interesting
milestone to achieve in the years to come.

###### Acknowledgements.

The authors thank M. Saunders and T. Choi for sharing their iterative MINRES-
QLP solver and providing useful insights, as well as A. Tichai for interesting
discussions. This project has received funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Skłodowska-
Curie grant agreement No. 839847. T. R. R. acknowledges the support of the the
Spanish MICINN under PGC2018-094583-B-I00. R. R. was supported by the DFG
Sonderforschungsbereich SFB 1245 (Project ID 279384907) and the BMBF
Verbundprojekt 05P2021 (ErUM-FSP T07, Contract No. 05P21RDFNB). H. H. was
supported by the U.S. Department of Energy, Office of Science, Office of
Nuclear Physics under Grant No. de-sc0017887. J.M.Y. was supported by the
Fundamental Research Funds for the Central Universities, Sun Yat-sen
University. The MR-IMSRG code uses optimizations done in collaboration with
the U.S. Department of Energy, Office of Science, Office of Advanced
Scientific Computing Research, Scientific Discovery through Advanced Computing
(SciDAC) Program FASTMath and RAPIDS2 Institutes. Calculations were performed
by using HPC resources from GENCI-TGCC (Contract No. A009057392), CETA-Ciemat
(FI-2021-2-0013), the Lichtenberg high performance computer of the Technische
Universität Darmstadt and the Institute for Cyber-Enabled Research (ICER) at
Michigan State University.

##  Appendix A 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

$\displaystyle{\mathbb{M}}{\mathbb{a}}$ $\displaystyle=-{\mathbb{h}_{1}}\,.$
(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
${\mathbb{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^{\prime})$ applicable to $I$ and $J$, the antisymmetry properties
are given by

$\displaystyle M_{pIqJ}$
$\displaystyle=\epsilon(\tau)\epsilon(\tau^{\prime})M_{p\tau(I)q\tau^{\prime}(J)}\,,$
(2a) $\displaystyle a^{J}(q)$ $\displaystyle=\epsilon(\tau)a^{\tau(J)}(q)\,,$
(2b) $\displaystyle h_{1}^{I}(p)$
$\displaystyle=\epsilon(\tau^{\prime})h_{1}^{\tau^{\prime}(I)}(p)\,,$ (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}$

$\displaystyle\sum_{q}\sum_{\mathcal{J}}\sum_{J\in\mathcal{J}}M_{p\bar{I}qJ}a^{J}(q)$
$\displaystyle=-h_{1}^{\bar{I}}(p)\,,$ (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

$\displaystyle\sum_{q}\sum_{\mathcal{J}}|\mathcal{J}|M_{p\bar{I}\
q\bar{J}}a^{\bar{J}}(q)$ $\displaystyle=-h_{1}^{\bar{I}}(p)\,.$ (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

$\displaystyle\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)$
$\displaystyle=-\sqrt{|\mathcal{I}|}h_{1}^{\bar{I}}(p)\,,$ (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}$.

##  Appendix B 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.

###  B.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 App. 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)

$\mathbf{N}=\mathbf{X}\mathbf{I}\mathbf{X}^{\dagger}\,,$ (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

$\displaystyle\mathbf{\tilde{M}}$
$\displaystyle\equiv\mathbf{X}\mathbf{M}\mathbf{X}^{\dagger}\,,$ (7a)
$\displaystyle\mathbf{\tilde{a}}$
$\displaystyle\equiv\mathbf{X}^{\dagger}\mathbf{a}\,,$ (7b)
$\displaystyle\mathbf{\tilde{h}}_{1}$
$\displaystyle\equiv\mathbf{X}^{\dagger}\mathbf{h}_{1}\,,$ (7c)

such that the linear problem equivalently reads

$\mathbf{\tilde{M}}\mathbf{\tilde{a}}=-\mathbf{\tilde{h}}_{1}\,.$ (8)

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

$\mathbf{\Delta}=\mathbf{Y}^{\dagger}\mathbf{\tilde{M}}\mathbf{Y}\,,$ (9)

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

$E^{(2)}=-\mathbf{h_{1}}^{\dagger}\mathbf{X}\mathbf{Y}\mathbf{\Delta}^{-1}\mathbf{Y}^{\dagger}\mathbf{X}^{\dagger}\mathbf{h_{1}}\,.$
(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.

Figure 17: (Color online) 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$.

###  B.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 Businger1971  and QLP
QLP  factorizations, which are briefly discussed in the following.

####  B.2.1 Pivoted QR

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

$AD=QR\,,$ (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$.

####  B.2.2 Pivoted QLP

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

$D^{\prime}AD=QLP$ (12)

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

$L=\begin{pmatrix}\tilde{L}&0\\\ 0&0\end{pmatrix},$ (13)

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

####  B.2.3 Illustration of non-iterative solvers

Figure 18: (Color online) 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$.

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\|\mathbf{M}\mathbf{a}+\mathbf{h}_{1}\|\|\mathbf{a}\|$ (see Sec. 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.

###  B.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 MINRES  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 MINRES-QLP .

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.

####  B.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

$\mathbf{M}\mathbf{x}=-\mathbf{h}_{1}\,,$ (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

$\mathbf{C}^{-1}\mathbf{M}\mathbf{C}^{-T}\mathbf{y}=-\mathbf{C}^{-1}\mathbf{h}_{1}\,.$
(15)

Whereas various techniques are available to build an appropriate matrix
$\mathbf{A}$, designing efficient preconditioners is still an active field of
research WATHEN . 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

$\mathbf{D}^{-\frac{1}{2}}\mathbf{M}\mathbf{D}^{-\frac{1}{2}}\mathbf{y}=-\mathbf{D}^{-\frac{1}{2}}\mathbf{h}_{1}\,.$
(16)

If the scaled matrix

$\mathbf{\tilde{M}}\equiv\mathbf{D}\mathbf{M}\mathbf{D}$ (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. BINORM
amounts to scaling all rows and columns to unit norm, which can yield
significantly better results at low cost. A stochastic matrix-free variant
Bradley2010  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

$\mathbf{N}\sim\mathbf{L}\mathbf{L}^{T},$ (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

$\mathbf{M}\mathbf{a}=-\mathbf{h}_{1}\,,$ (19)

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

$\mathbf{N}\mathbf{M}\mathbf{N}\left(\mathbf{N}^{-1}\mathbf{a}\right)=-\mathbf{N}\mathbf{h}_{1}\,.$
(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}$.

####  B.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

$\mathbf{M}\mathbf{a}=-\mathbf{h}_{1}+\mathbf{b}\,,$ (21)

the second-order energy evaluated with Hylleraas’ functional reads

$\displaystyle E^{(2)}$
$\displaystyle=\mathbf{a}^{\dagger}\mathbf{M}\mathbf{a}+\mathbf{h}_{1}^{\dagger}\mathbf{a}+\mathbf{a}^{\dagger}\mathbf{h}_{1}$
$\displaystyle=\mathbf{a}^{\dagger}\mathbf{b}+\mathbf{h}_{1}^{\dagger}\mathbf{a}\,.$
(22)

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

$\delta
E^{(2)}=\mathbf{a}^{\dagger}\left(\mathbf{M}\mathbf{a}+\mathbf{h}_{1}\right)\,.$
(23)

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

$|\delta
E|^{(2)}\equiv\|\mathbf{a}\|\|\mathbf{M}\mathbf{a}+\mathbf{h}_{1}\|\,.$ (24)

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

$|\delta
E|^{(2)}\equiv\|\mathbf{N}^{-1}\mathbf{a}\|\|\mathbf{N}\mathbf{M}\mathbf{N}\left(\mathbf{N}^{-1}\mathbf{a}\right)+\mathbf{N}\mathbf{h}_{1}\|\,.$
(25)

####  B.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

$r\equiv\frac{\|\mathbf{M}\mathbf{a}+\mathbf{h}_{1}\|}{\|\mathbf{a}\|}.$ (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.

####  B.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.

Figure 19: (Color online) 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 Huther_2020 ;
Entem:2017gor , $\lambda_{\mathrm{srg}}=1.88\,\mathrm{fm}^{-1}$,
$\hbar\omega=20\,\mathrm{MeV}$ and $e_{\mathrm{max}}=2$. Figure 20: (Color
online) 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 Huther_2020 ; Entem:2017gor ,
$\lambda_{\mathrm{srg}}=1.88\,\mathrm{fm}^{-1}$,
$\hbar\omega=20\,\mathrm{MeV}$ and $e_{\mathrm{max}}=2$.

##  Appendix C Complex-shift method

###  C.1 Motivations

As it appears in Eq. (10), the second-order energy relies on the invertibility
of $\mathbf{\Delta}$ to generate non-zero energy denominators. However, the
eigenvalues of $\mathbf{\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 burton20a ;
FORSBERG1997 .

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

$\bar{\mathbf{\Delta}}\equiv\mathbf{\Delta}+\imath\gamma I\,,$ (27)

or, equivalently, working in the original basis

$\bar{\mathbf{M}}\equiv\mathbf{M}+\imath\gamma\mathbf{N}\,,$ (28)

which simply corresponds to adding a complex term to the unperturbed
Hamiltonian $H_{0}$. The imaginary shift moves zero eigenvalues of
$\mathbf{\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,

$E^{(2)}=\Re\left[\mathbf{a}^{\dagger}\bar{\mathbf{M}}\mathbf{a}+\mathbf{a}^{\dagger}\mathbf{h}_{1}+\mathbf{h}_{1}^{\dagger}\mathbf{a}\right].$
(29)

###  C.2 Implementation in real arithmetic

Although an extension of MINRES-QLP has been developed to handle complex
symmetric matrices CS-MINRES-QLP ; MINRES-IMPL , 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

$(\mathbf{M}+\imath\gamma\mathbf{N})(\mathbf{a}+\imath\mathbf{b})=-\mathbf{h}_{1}$
(30)

is recast into a blockwise 2x2 real symmetric system

$\begin{pmatrix}\mathbf{M}&-\gamma\mathbf{N}\\\
-\gamma\mathbf{N}&-\mathbf{M}\end{pmatrix}\begin{pmatrix}\mathbf{a}\\\
\mathbf{b}\end{pmatrix}=\begin{pmatrix}-\mathbf{h}_{1}\\\
\mathbf{0}\end{pmatrix}\,.$ (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.

###  C.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).

Figure 21: (Color online) 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 Huther_2020 ;
Entem:2017gor , $\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$ MeV282828In the present case, $\gamma=4$ MeV constitutes the lowest
value that is empirically found to deliver a controlled numerical result.
Below this value, the energy curve displays an erratic behavior due to the
occurrence of an intruder state. 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.

Figure 22: (Color online) 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 Huther_2020 ;
Entem:2017gor , $\lambda_{\mathrm{srg}}=1.88\,\mathrm{fm}^{-1}$,
$\hbar\omega=20\,\mathrm{MeV}$ and $e_{\mathrm{max}}=2$. Figure 23: (Color
online) 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 Huther_2020 ;
Entem:2017gor , $\lambda_{\mathrm{srg}}=1.88\,\mathrm{fm}^{-1}$,
$\hbar\omega=20\,\mathrm{MeV}$ and $e_{\mathrm{max}}=2$.

##  Appendix D Discussion on numerics

###  D.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 Tab. 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.

Method | HFB | PGCM | BMBPT(2) | BMBPT(3) | PGCM-PT(2) | FCI
---|---|---|---|---|---|---
Runtime | $O(N^{4})$ | $O(n_{\text{proj}}n_{\text{gcm}}^{2}N^{4})$ | $O(N^{5})$ | $O(N^{6})$ | $O(n_{\text{proj}}n_{\text{gcm}}^{2}N^{8})$ | $O(N^{\text{A}})$
Storage | $O(N^{4})$ | $O(N^{4})$ | $O(N^{4})$ | $O(N^{4})$ | $O(n_{\text{gcm}}^{2}N^{8})$ | $O(N^{\text{A}})$
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.

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.

Figure 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.

###  D.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. Tichai_2019 ; Garniron_2018
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 porro2021importance  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.

####  D.2.1 Norm-based importance truncation

##### Exact arithmetic

The norm of a projected configuration is

$n^{I}(p)\equiv\langle\Omega^{I}(p)|\Omega^{I}(p)\rangle\in[0,1]$ (32)

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

$|\Omega^{I}(p)\rangle=0\,.$ (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.

####  D.2.2 Hamiltonian-based importance truncation

##### Exact arithmetic

The Hamiltonian matrix element of a projected configuration reads

$h_{1}^{I}(p)\equiv\langle\Omega^{I}(p)|H_{1}|\Theta^{(0)}\rangle\,.$ (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$.

####  D.2.3 Energy-based importance truncation

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

$e^{(2)I}(q)=h_{1}^{I*}(q)a^{I}(q)\,.$ (35)

A configuration $(q,I)$ for which $e^{(2)I}(q)=0$ does not contribute to the
correlation energy292929Such a configuration might still contribute indirectly
by influencing the value of the other coefficients $\\{a^{J}(q^{\prime})\\}$..
Removing configurations based on the size of their contribution to the
correlation energy corresponds to the method advocated in Refs. Tichai_2019 ;
porro2021importance . 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
$|\Phi(q)\rangle$.

## References

  * (1) M. Frosini, T. Duguet, J.-P. Ebran, V. Somà, Multi-reference many-body perturbation theory for nuclei I – Novel PGCM-PT formalism (2021).  arXiv:2110.15737.
  * (2) M. Frosini, T. Duguet, J.-P. Ebran, B. Bally, T. Mongelli, T. R. Rodríguez, R. Roth, V. Somà, Multi-reference many-body perturbation theory for nuclei II – Ab initio study of neon isotopes via PGCM and IM-NCSM calculations (2021).  arXiv:2111.00797.
  * (3) K. Tsukiyama, S. K. Bogner, A. Schwenk, In-Medium Similarity Renormalization Group for Nuclei, Phys. Rev. Lett. 106 (2011) 222502.  arXiv:1006.3639, doi:10.1103/PhysRevLett.106.222502.
  * (4) H. Hergert, S. K. Bogner, T. D. Morris, A. Schwenk, K. Tsukiyama, The In-Medium Similarity Renormalization Group: A Novel Ab Initio Method for Nuclei, Phys. Rept. 621 (2016) 165–222.  arXiv:1512.06956, doi:10.1016/j.physrep.2015.12.007.
  * (5) H. Hergert, S. K. Bogner, T. D. Morris, S. Binder, A. Calci, J. Langhammer, R. Roth, Ab initio multireference in-medium similarity renormalization group calculations of even calcium and nickel isotopes, Phys. Rev. C 90 (4) (2014) 041302\.  arXiv:1408.6555, doi:10.1103/PhysRevC.90.041302.
  * (6) H. Hergert, S. K. Bogner, J. G. Lietz, T. D. Morris, S. Novario, N. M. Parzuchowski, F. Yuan, In-Medium Similarity Renormalization Group Approach to the Nuclear Many-Body Problem, Lect. Notes Phys. 936 (2017) 477–570.  arXiv:1612.08315, doi:10.1007/978-3-319-53336-0_10.
  * (7) H. Hergert, In-Medium Similarity Renormalization Group for Closed and Open-Shell Nuclei, Phys. Scripta 92 (2) (2017) 023002.  arXiv:1607.06882, doi:10.1088/1402-4896/92/2/023002.
  * (8) S. K. Bogner, R. J. Furnstahl, A. Schwenk, From low-momentum interactions to nuclear structure, Prog. Part. Nucl. Phys. 65 (2010) 94–147.  doi:10.1016/j.ppnp.2010.03.001.
  * (9) R. Roth, J. Langhammer, A. Calci, S. Binder, P. Navrátil, Similarity-transformed chiral $nn+3n$ interactions for the ab initio description of ${}^{12}\mathbf{C}$ and ${}^{16}\mathbf{O}$, Phys. Rev. Lett. 107 (2011) 072501.  doi:10.1103/PhysRevLett.107.072501.
URL https://link.aps.org/doi/10.1103/PhysRevLett.107.072501

  * (10) R. Roth, A. Calci, J. Langhammer, S. Binder, Evolved chiral $nn+3n$ hamiltonians for ab initio nuclear structure calculations, Phys. Rev. C 90 (2014) 024325.  doi:10.1103/PhysRevC.90.024325.
URL https://link.aps.org/doi/10.1103/PhysRevC.90.024325

  * (11) M. Frosini, J.-P. Ebran, N. Dubray, A. Porro, T. Duguet, V. Somà, unpublished (2021).
  * (12) M. Frosini, J.-P. Ebran, A. Porro, T. Duguet, V. Somà, unpublished (2021).
  * (13) T. Hüther, K. Vobig, K. Hebeler, R. Machleidt, R. Roth, Family of chiral two- plus three-nucleon interactions for accurate nuclear structure studies, Physics Letters B 808 (2020) 135651.  doi:10.1016/j.physletb.2020.135651.
URL http://dx.doi.org/10.1016/j.physletb.2020.135651

  * (14) D. R. Entem, R. Machleidt, Y. Nosyk, High-quality two-nucleon potentials up to fifth order of the chiral expansion, Phys. Rev. C 96 (2) (2017) 024004.  arXiv:1703.05454, doi:10.1103/PhysRevC.96.024004.
  * (15) T. Duguet, A. Signoracci, Symmetry broken and restored coupled-cluster theory. II. Global gauge symmetry and particle number, J. Phys. G 44 (1) (2017) 015103, [Erratum: J.Phys.G 44, 049601 (2017)].  arXiv:1512.02878, doi:10.1088/0954-3899/44/1/015103.
  * (16) A. Tichai, P. Arthuis, T. Duguet, H. Hergert, V. Somà, R. Roth, Bogoliubov Many-Body Perturbation Theory for Open-Shell Nuclei, Phys. Lett. B 786 (2018) 195.  doi:10.1016/j.physletb.2018.09.044.
  * (17) P. Arthuis, T. Duguet, A. Tichai, R.-D. Lasseri, J.-P. Ebran, ADG: Automated generation and evaluation of many-body diagrams I. Bogoliubov many-body perturbation theory, Comput. Phys. Commun. 240 (2019) 202.  doi:10.1016/j.cpc.2018.11.023.
  * (18) P. Demol, M. Frosini, A. Tichai, V. Somà, T. Duguet, Bogoliubov many-body perturbation theory under constraint, Annals Phys. 424 (2021) 168358.  doi:10.1016/j.aop.2020.168358.
  * (19) A. Tichai, R. Roth, T. Duguet, Many-body perturbation theories for finite nuclei, Front. Phys. 8 (2020) 164.  doi:10.3389/fphy.2020.00164.
  * (20) T. Duguet, B. Bally, A. Tichai, Zero-pairing limit of Hartree-Fock-Bogoliubov reference states, Phys. Rev. C 102 (5) (2020) 054320.  arXiv:2006.02871, doi:10.1103/PhysRevC.102.054320.
  * (21) B. Ladóczki, M. Uejima, S. L. Ten-no, Third-order epstein–nesbet perturbative correction to the initiator approximation of configuration space quantum monte carlo, The Journal of Chemical Physics 153 (11) (2020) 114112.  arXiv:https://doi.org/10.1063/5.0022101, doi:10.1063/5.0022101.
URL https://doi.org/10.1063/5.0022101

  * (22) K. Hebeler, S. K. Bogner, R. J. Furnstahl, A. Nogga, A. Schwenk, Improved nuclear matter calculations from chiral low-momentum interactions, Phys. Rev. C 83 (2011) 031301.  doi:10.1103/PhysRevC.83.031301.
URL https://link.aps.org/doi/10.1103/PhysRevC.83.031301

  * (23) A. Nogga, S. K. Bogner, A. Schwenk, Low-momentum interaction in few-nucleon systems, Phys. Rev. C 70 (2004) 061002.  doi:10.1103/PhysRevC.70.061002.
URL https://link.aps.org/doi/10.1103/PhysRevC.70.061002

  * (24) E. Gebrerufael, K. Vobig, H. Hergert, R. Roth, Ab initio description of open-shell nuclei: Merging no-core shell model and in-medium similarity renormalization group, Phys. Rev. Lett. 118 (2017) 152503.  doi:10.1103/PhysRevLett.118.152503.
URL https://link.aps.org/doi/10.1103/PhysRevLett.118.152503

  * (25) J. M. Yao, J. Engel, L. J. Wang, C. F. Jiao, H. Hergert, Generator-coordinate reference states for spectra and $0\nu\beta\beta$ decay in the in-medium similarity renormalization group, Phys. Rev. C 98 (5) (2018) 054311.  doi:10.1103/PhysRevC.98.054311.
  * (26) J. M. Yao, B. Bally, J. Engel, R. Wirth, T. R. Rodríguez, H. Hergert, Ab initio treatment of collective correlations and the neutrinoless double beta decay of ${}^{48}\mathrm{Ca}$, Phys. Rev. Lett. 124 (2020) 232501.  doi:10.1103/PhysRevLett.124.232501.
URL https://link.aps.org/doi/10.1103/PhysRevLett.124.232501

  * (27) B. Bally, B. Avez, M. Bender, P.-H. Heenen, Beyond mean-field calculations for odd-mass nuclei, Phys. Rev. Lett. 113 (2014) 162501\.  doi:10.1103/PhysRevLett.113.162501.
URL https://link.aps.org/doi/10.1103/PhysRevLett.113.162501

  * (28) M. Borrajo, T. R. Rodríguez, J. Luis Egido, Symmetry conserving configuration mixing method with cranked states, Physics Letters B 746 (2015) 341–346.  doi:https://doi.org/10.1016/j.physletb.2015.05.030.
URL https://www.sciencedirect.com/science/article/pii/S0370269315003676

  * (29) H. Hergert, A guided tour of ab initio nuclear many-body theory, Frontiers in Physics 8 (2020) 379\.  doi:10.3389/fphy.2020.00379.
URL https://www.frontiersin.org/article/10.3389/fphy.2020.00379

  * (30) P. Businger, G. H. Golub, Linear Least Squares Solutions by Housholder Transformations, Springer Berlin Heidelberg, Berlin, Heidelberg, 1971, pp. 111–118.  doi:10.1007/978-3-662-39778-7_8.
URL https://doi.org/10.1007/978-3-662-39778-7_8

  * (31) G. W. Stewart, The qlp approximation to the singular value decomposition, SIAM Journal on Scientific Computing 20 (4) (1999) 1336–1348.  arXiv:https://doi.org/10.1137/S1064827597319519, doi:10.1137/S1064827597319519.
URL https://doi.org/10.1137/S1064827597319519

  * (32) C. C. Paige, M. A. Saunders, Solution of sparse indefinite systems of linear equations, SIAM Journal on Numerical Analysis 12 (4) (1975) 617–629.  arXiv:https://doi.org/10.1137/0712047, doi:10.1137/0712047.
URL https://doi.org/10.1137/0712047

  * (33) S.-C. T. Choi, C. C. Paige, M. A. Saunders, Minres-qlp: A krylov subspace method for indefinite or singular symmetric systems, SIAM Journal on Scientific Computing 33 (4) (2011) 1810–1836.  doi:10.1137/100787921.
URL http://dx.doi.org/10.1137/100787921

  * (34) A. J. Wathen, Preconditioning, Acta Numerica 24 (2015) 329–376.  doi:10.1017/S0962492915000021.
  * (35) O. Livne, G. Golub, Scaling by binormalization, Numerical Algorithms 35 (2004) 97–120.  doi:10.1023/B:NUMA.0000016606.32820.69.
  * (36) A. Bradley, Algorithms for the equilibration of matrices and their application to limited-memory quasi-newton methods (2010).
  * (37) H. G. A. Burton, A. J. W. Thom, J. Chem. Theory Comput. 16 (4) (2020) 5586.
  * (38) N. Forsberg, P. Åke Malmqvist, Multiconfiguration perturbation theory with imaginary level shift, Chemical Physics Letters 274 (1) (1997) 196–204.  doi:https://doi.org/10.1016/S0009-2614(97)00669-6.
URL https://www.sciencedirect.com/science/article/pii/S0009261497006696

  * (39) Sou-Cheng T. Choi, Minimal Residual Methods for Complex Symmetric, Skew Symmetric, and Skew Hermitian Systems, arXiv e-prints (2013) arXiv:1304.6782arXiv:1304.6782.
  * (40) Sou-Cheng, Choi, CS-MINRES-QLP, version 1” [Matlab/GNU-Octave Software] (2013).
  * (41) A. Tichai, J. Ripoche, T. Duguet, Pre-processing the nuclear many-body problem, The European Physical Journal A 55 (6) (Jun 2019).  doi:10.1140/epja/i2019-12758-6.
URL http://dx.doi.org/10.1140/epja/i2019-12758-6

  * (42) Y. Garniron, A. Scemama, E. Giner, M. Caffarel, P.-F. Loos, Selected configuration interaction dressed by perturbation, The Journal of Chemical Physics 149 (6) (2018) 064103\.  doi:10.1063/1.5044503.
URL http://dx.doi.org/10.1063/1.5044503

  * (43) A. Porro, V. Somà, A. Tichai, T. Duguet, Importance truncation in non-perturbative many-body techniques (2021).  arXiv:2103.14544.