Pre-processing the nuclear many-body problem: Importance truncation versus tensor factorization techniques

TL;DR

This paper compares importance truncation and tensor factorization as pre-processing methods to reduce the complexity of nuclear many-body calculations, showing importance truncation's potential for high-accuracy, low-cost computations.

Contribution

It introduces and benchmarks importance truncation against tensor factorization for nuclear many-body problem pre-processing, highlighting importance truncation's practical advantages.

Findings

Importance truncation shows great potential for reducing computational costs.

Tensor factorization requires significant numerical development for large spaces.

Both methods perform well in small model space benchmarks.

Abstract

The solution of the nuclear A-body problem encounters severe limitations from the size of many-body operators. These limitations are typically related to both the (iterative) storing of the associated tensors and to the computational time related to their multiple contractions in the calculation of various quantities of interest. However, not all the degrees of freedom encapsulated into these tensors equally contribute to the description of many-body observables. Identifying systematic and dominating patterns, a relevant objective is to achieve an \emph{a priori} reduction to the most relevant degrees of freedom via a pre-processing of the A-body problem. The present paper is dedicated to the analysis of two different paradigms to do so. The factorization of tensors in terms of lower-rank ones, whose know-how has been recently transferred to the realm of nuclear structure, is compared to a reduction of the tensors' index size based on an importance truncation. While the objective is to eventually utilize these pre-processing tools in the context of non-perturbative many-body methods, benchmark calculations are presently performed within the frame of perturbation theory. More specifically, we employ the recently introduced Bogoliubov many-body perturbation theory that is systematically applicable to open-shell nuclei displaying strong correlations. This extended perturbation theory serves as a jumpstart for non-perturbative Bogoliubov coupled cluster and Gorkov self-consistent Green's function theories. Results obtained in "small" model spaces are equally encouraging for tensor factorization and importance truncation techniques. While the former requires significant numerical developments to be applied in large model spaces, the latter is presently applied in this context and demonstrates great potential to enable high-accuracy calculations at a much reduced computational cost.

Tables

Table 1. Table 1: Number of entries and associated memory of J 𝐽 J -coupled mode-4 and mode-6 tensors expressed in the U ​ ( 1 ) 𝑈 1 U(1) -breaking quasi-particle basis as a function of the truncation parameter e max subscript 𝑒 max e_{\text{max}} of the initial one-body spherical HO basis. Two storage schemes are employed (see text). Storage estimates assume double precision for all tensor entries. The quoted numbers correspond to the complete tensor, i.e. they sum contributions from all possible J 𝐽 J ( J 12 , J 45 , J subscript 𝐽 12 subscript 𝐽 45 𝐽 J_{12},J_{45},J ) blocks.

	Index size		${}_{}{}^J\tilde{T}_{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4}^{40}$				${}_{}{}^{J_{12}{J}_{45}J}\tilde{T}_{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4{\tilde{k}}_5{\tilde{k}}_6}^{60}$
$e_{\text{max}}$	$N$	$\tilde{N}$	Naive		Optimal		Naive		Optimal
2	40	12	$1.2\cdot {10}^5$	$1.0$ Mb	$6.5\cdot {10}^2$	$5$ kb	$4.3\cdot {10}^8$	$3.1$ Gb	$2.2\cdot {10}^5$	$0.1$ Mb
4	140	30	$8.1\cdot {10}^6$	$61.8$ Mb	$1.9\cdot {10}^4$	$0.1$ Mb	$7.0\cdot {10}^{11}$	$5.0$ Tb	$7.3\cdot {10}^6$	$54.4$ Mb
6	336	56	$1.4\cdot {10}^8$	$1.0$ Gb	$2.4\cdot {10}^5$	$1.8$ Mb	$1.6\cdot {10}^{14}$	$1.1\cdot {10}^3$ Tb	$5.0\cdot {10}^8$	$3.8$ Gb
8	660	90	$1.2\cdot {10}^9$	$8.8$ Gb	$1.7\cdot {10}^6$	$13.0$ Mb	$1.5\cdot {10}^{16}$	$1.1\cdot {10}^6$ Tb	$1.4\cdot {10}^{10}$	$103.1$ Gb
10	1140	132	$6.7\cdot {10}^9$	$49.8$ Gb	$8.7\cdot {10}^6$	$66.0$Mb	$6.9\cdot {10}^{17}$	$5.0\cdot {10}^7$ Tb	$2.1\cdot {10}^{11}$	$1.5$ Tb
12	1820	182	$2.9\cdot {10}^{10}$	$212.5$ Gb	$3.4\cdot {10}^7$	$261.5$ Mb	$9.2\cdot {10}^{18}$	$6.7\cdot {10}^8$ Tb	$2.1\cdot {10}^{12}$	$15.2$ Tb

Table 2. Table 2: Number of different angular-momentum blocks of J 𝐽 J -coupled mode-4 and mode-6 tensors in the two storage formats.

	${}_{}{}^J\tilde{T}_{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4}^{40}$		${}_{}{}^{J_{12}{J}_{45}J}\tilde{T}_{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4{\tilde{k}}_5{\tilde{k}}_6}^{60}$
$e_{\text{max}}$	Naive	Optimal	Naive	Optimal
2	6	6	288	75
4	10	10	1400	289
6	14	14	3920	727
8	18	18	8424	1469
10	22	22	15488	2595
12	26	26	25688	4185

Table 3. Table 3: BMBPT contributions to the binding energy of 18 O as a function of e max subscript 𝑒 max e_{\text{max}} ; i.e. first-order (HFB) contribution, second-order correction and fourth-order correction from the leading triple BCC amplitude. Because the original matrix elements of the three-nucleon force in Eq. ( 1 ) are limited by an a priori e 3 ​ max = 14 subscript 𝑒 3 max 14 e_{3\text{max}}=14 truncation, numbers marked by † are slightly approximate compared to the calculation that would consistently use e 3 ​ max = 3 ​ e max subscript 𝑒 3 max 3 subscript 𝑒 max e_{3\text{max}}=3\,e_{\text{max}} . Furthermore, numbers marked by ∗ have been obtained by imposing an initial e 3 ​ max = 22 subscript 𝑒 3 max 22 e_{3\text{max}}=22 reduction on the leading triple BCC amplitude tensor.

$e_{\text{max}}$	$E_{\text{HFB}}$ $\left[\text{MeV}\right]$	$\Delta {\Omega }_0^{\left(2\right)}$ $\left[\text{MeV}\right]$	$\Delta {\Omega }_0^{\left[4_{\text{T}}\right]}$ $\left[\text{MeV}\right]$
2	-100.481	-10.328	1.228
4	-105.994	-25.815	2.493
6	-107.289†	-29.316†	2.961†
8	-107.588†	-29.834†	3.022†
10	-107.802†	-29.875†	2.962†∗
12	-108.051†	-29.783†	2.956†∗

Table 4. Table 4: Computational scaling of state-of-the-art expansion many-body methods. The displayed numbers correspond to optimized contraction patterns obtained via the introduction of intermediates in MBPT and CC or via the use of the natural basis in IMSRG.

Truncation scheme	Scaling	(B)MBPT(n)
(B)MBPT(2)	$N^4$	2
(B)MBPT(3)	$N^6$	3
(B)MBPT(4)	$N^7$	4
(B)CCSD	$N^6$	3
(B)CCSDT	$N^8$	4
IMSRG(2)	$N^6$	3
ADC(3)	$N^6$	3

Equations

$$H_{\text{nucl}}$$

$$\equiv \frac{1}{(1!{)}^2}\sum _{pq}{t}_{pq}{c}_p^{†}{c}_q$$

$$+\frac{1}{(2!{)}^2}\sum _{pqrs}{\bar{v}}_{pqrs}{c}_p^{†}{c}_q^{†}{c}_s{c}_r$$

$$+\frac{1}{(3!{)}^2}\sum _{pqrstu}{\bar{w}}_{pqrstu}{c}_p^{†}{c}_q^{†}{c}_r^{†}{c}_u{c}_t{c}_s$$

$$+...,$$

$$E_0$$

$${\beta }_k$$

$${\beta }_k^{†}$$

$$\Omega \equiv H-\lambda A.$$

$$\Omega$$

$$+\stackrel{{\Omega }^{[2]}}{\overbrace{{\Omega }^{20}+{\Omega }^{11}+{\Omega }^{02}}}$$

$$+\stackrel{{\Omega }^{[4]}}{\overbrace{{\Omega }^{40}+{\Omega }^{31}+{\Omega }^{22}+{\Omega }^{13}+{\Omega }^{04}}}$$

$$+\stackrel{{\Omega }^{[6]}}{\overbrace{{\Omega }^{60}+{\Omega }^{51}+{\Omega }^{42}+{\Omega }^{33}+{\Omega }^{24}+{\Omega }^{15}+{\Omega }^{06}}},$$

$${\Omega }^{31}$$

$${\Omega }_{k_1\dots k_i{k}_{i+1}\dots k_{i+j}}^{ij}$$

$${\Omega }_{k_1{k}_2}^{20}$$

$${\Omega }_{k_1{k}_2}^{11}$$

$$|\Psi ⟩\equiv e^{\mathcal{T}}|\Phi ⟩,$$

$${\mathcal{T}}_1$$

$${\mathcal{T}}_2$$

$${\mathcal{T}}_3$$

$${\mathcal{E}}_0$$

$${\mathcal{M}}_{k_1\dots k_p}\equiv ⟨{\Phi }^{k_1\dots k_p}|\bar{\Omega }|\Phi ⟩=0,$$

$$|{\Phi }^{k_1\dots k_p}⟩\equiv {\beta }_{k_1}^{†}\dots {\beta }_{k_p}^{†}|\Phi ⟩.$$

$$t_{k_1{k}_2}^{20(1)}$$

$$t_{k_1{k}_2{k}_3{k}_4}^{40(1)}$$

$$t_{k_1{k}_2{k}_3{k}_4{k}_5{k}_6}^{60(1)}$$

$$E_{k_1{k}_2{k}_3{k}_4\dots }$$

$$k\equiv (n_k,l_k,j_k,m_k,t_k),$$

$${{}^J\tilde{\Omega }}_{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4}^{ij}\equiv ⟨{\tilde{k}}_1{\tilde{k}}_2(J)|{\tilde{\Omega }}^{ij}|{\tilde{k}}_3{\tilde{k}}_4(J)⟩,$$

$$\tilde{k}\equiv (n_k,l_k,j_k,t_k),$$

$$\frac{1}{4!}\sum _J{\hat{J}}^2\sum _{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4}{}^J{\tilde{\Omega }}_{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4}^{04}{}^J{\tilde{t}}_{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4}^{40},$$

$$\Delta {\Omega }_0^{(2)}=\frac{1}{4!}\sum _J{\hat{J}}^2\sum _{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4}{}^J{\tilde{\Omega }}_{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4}^{04}{}^J{\tilde{t}}_{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4}^{40(1)},$$

$${}^J{\tilde{t}}_{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4}^{40(1)}=-\frac{{}^J{\tilde{\Omega }}_{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4}^{40}}{E_{{\tilde{k}}_1{\tilde{k}}_2{\tilde{k}}_3{\tilde{k}}_4}}.$$

$$\Vert T\Vert \equiv \sqrt{\sum _{i_1...i_n}{T}_{i_1...i_n}{T}_{i_1...i_n}^{\ast }}.$$

$$\Delta T\equiv \frac{\Vert T-\hat{T}\Vert }{\Vert T\Vert }.$$

$$T_{i_1{i}_2{i}_3{i}_4}=\sum _{\alpha \beta }{X}_{i_1\alpha }^1{X}_{i_2\alpha }^2{W}_{\alpha \beta }{X}_{i_3\beta }^3{X}_{i_4\beta }^4.$$

Full text

1 Introduction

Due to the growing amount of data at play in applied and fundamental sciences there exists a strong need for efficient data analysis tools to store and process more efficiently the underlying information. This statement, generally valid independently of the type of data or research area, applies acutely to the nuclear many-body problem that provides a particularly interesting challenge beyond the lightest nuclei.

Over the past decade, tremendous progress have been made to extend the reach of ab initio methods to larger mass numbers and better accuracy. In particular the development of many-body methods based on a systematic expansion of the exact solution around a conveniently chosen reference state has allowed an efficient description of medium-mass (semi-)magic nuclei. Examples of such methods are many-body perturbation theory (MBPT) Langhammer:2012jx ; Hu:2016txm ; Tichai:2017rqe ; Tichai:2018mll ; Arthuis:2018yoo ; Hu18arxiv , self-consistent Green’s function (SCGF) Dickhoff:2004xx ; Soma:2011aj ; Soma:2013xha ; Carbone:2013eqa ; Lapoux:2016exf ; Duguet:2016wwr ; Raimondi:2017kzi ; Raimondi:2018mtv , coupled-cluster (CC) Hagen:2013nca ; Signoracci:2014dia ; Morris:2017vxi or the in-medium similarity renormalization group (IM-SRG) Hergert:2015awm ; Hergert:2016iju ; Parzuchowski:2017wcq ; Morris:2017vxi approaches. This category of polynomially-scaling many-body methods is the focus of our attention.

The guiding idea behind these methods is to account for dynamical correlations via a particle-hole (quasi-particle) expansion that is approximated according to a chosen truncation scheme. In the simplest case, low-order MBPT corresponds to keeping the first few terms of an expansion in powers of the residual interaction. In more advanced non-perturbative frameworks, all-order resummations of MBPT contributions are typically accounted for by solving a non-linear set of equations, e.g., the amplitude equations in CC theory. To extend their reach to open-shell nuclei and capture noticeably challenging static correlations, these single-reference methods employ a symmetry-breaking reference state Soma:2011aj ; Signoracci:2014dia ; Tichai:2018mll and possibly restore the broken symmetry in a two-step approach Duguet:2014jja ; Duguet:2015yle ; Qiu:2018edx .

In all expansion methods, the working equations are eventually expressed as multiple contractions of two different sets of tensors whose indices relate to a given (truncated) basis of the one-body Hilbert space ${\cal H}_{1}$. The first set of tensors relates to the given of the many-body Hamiltonian

[TABLE]

where $T$ indicates the kinetic energy whereas $V$ and $W$ characterize two- and three-body interactions, respectively. Dots embody possible higher-body operators. The Hamiltonian is thus represented via a set of mode-n (i.e. n-index) tensors, with $n=2,4$ and $6$ in the above example. The storage cost of the Hamiltonian is thus dominated by the mode-$2k$ tensor defining the highest $k$-body interaction and scales as $N^{2k}$, where $N$ denotes the dimension of the underlying one-body basis used.

The second set of tensors at play is presently denoted as many-body tensors and is method specific. For example, in the Hartree-Fock- (HF) based CC or MBPT implemented on the basis of a Hamiltonian containing two-body interactions only, the ground-state binding energy can be written as

[TABLE]

where $t^{a}_{i}$ and $t^{ab}_{ij}$ denote so-called single and double cluster amplitudes. These many-body amplitudes respectively constitute mode-$2$ and mode-$4$ tensors that are contracted with the mode-$4$ tensor originating from the Hamiltonian to compute, e.g., the ground-state energy. Depending on the level at which MBPT or CC are implemented, the calculation of $t^{a}_{i}$ and $t^{ab}_{ij}$ may actually require to carry higher-mode tensors as intermediate objects whose storage scales as $N^{l}$ and whose computation scales as $N^{p}$. In the non-perturbative CC method, these many-body tensors need to be solved for and stored repeatedly. A similar identification of the relevant many-body tensors can be carried out for SCGF or IMSRG methods.

Eventually, the memory load and the computational cost of a given many-body implementation are respectively driven by the highest mode tensor carried in the calculation and by the complexity of the tensor network associated with the working equations. In non-perturbative approaches, this cost is significantly augmented by the iterative character of the method that requires a repeated computation and storage of possibly high-mode many-body tensors. It is to be noted that the computational load of a given many-body method is significantly more pronounced for nuclei than for, e.g., electronic systems at play in atomic physics or quantum chemistry. This is due (i) to the importance of (at least) three-nucleon interactions, (ii) to significant low-to-high momentum couplings encapsulated into the Hamiltonian tensors and (iii) to the dominance of open-shell nuclei displaying strong correlations. While the first point translates into the necessity to carry (at least) a mode-6 tensor from the outset, the second point renders necessary to use significantly larger one-body bases. The third point makes mandatory to design versatile many-body methods capable of consistently grasping dynamical and static correlations. Eventually, these three features currently forbid the application of polynomially-scaling methods to high accuracy (below 1 per cent error) as well as to nuclei with mass $A>100$ and/or with a doubly open-shell character.

Generically speaking, overcoming the curse of dimensionality associated with the exponential growth of the $A$-body Hilbert space is a long-standing problem that has given birth to numerous method- or system-specific strategies. In one way or another, the idea is to exploit an a priori knowledge based on physical considerations, educated guess or a simplified evaluation of the importance of selected configurations. This concerns for instance (a) recent advances in tensor network theory (TNT) Ver08 ; Sch11 ; Or14 where the localization of quantum entanglement over the $A$-body Hilbert space is exploited to describe low-dimensional systems, (b) the recent development of infinite-basis extrapolation techniques in nuclear physics FuHa12 ; Fu14 ; We15 where an educated guess on the long-distance behavior of the many-body wave function is exploited, or (c) the method at the core of the Monte Carlo shell model KoDe97 ; OtHo01 where the relevance of randomly generated configurations are tested ’on the fly’ via an appropriate importance measure.

The present work takes place in this context and wishes to test two systematic strategies appropriate to the many-body methods introduced above. The two types of pre-processing are

Tensor factorization (TF) techniques whose goal is to produce a low-rank decomposition of both the Hamiltonian and the many-body tensors at play in the method. The anticipated merit of these techniques is to effectively reduce the value $p$ of the $N^{p}$ storage and computational costs of the full-fledge implementation. 2. 2.

Importance truncation (IT) techniques whose goal is to discard entries in the many-body tensors at play in the many-body framework. This is done by estimating the importance of each tensor entry on the basis of a low-order, i.e. computationally cheap, estimate. The anticipated merit of this technique is to effectively reduce the value $N$ of the $N^{p}$ storage and computational costs of the full-fledge implementation.

Recently, TF techniques have shown to be a promising tool to lower the computational requirements of nuclear many-body calculations Tichai:2018eem dedicated to closed-shell nuclei. Two tensor decomposition formats have been applied to state-of-the-art nuclear Hamiltonians to obtain compressed tensors yielding accurate description of nuclear observables benchmarked within the well-tested second-order Hartree-Fock MBPT (HF-MBPT). The goal of the present paper is to extend this analysis to single-reference methods applicable to open-shell nuclei, i.e. to methods authorizing the spontaneous breaking of symmetries. More specifically, $U(1)$ global-gauge symmetry associated to particle-number conservation is presently authorized to break in order to tackle singly open-shell nuclei. This feature is characteristic of Bogoliubov many-body perturbation theory (BMBPT) Tichai:2018mll ; Arthuis:2018yoo , particle-number projected BMBPT (PBMBPT) Duguet:2015yle , Bogoliubov CC (BCC) theory Signoracci:2014dia ; Henderson:2014vka , particle-number projected BCC (PBCC) Duguet:2015yle ; Qiu:2018edx and Bogoliubov-based, i.e. Gorkov, SCGF (GSCGF) theory Soma:2011aj ; Soma:2013xha . The present challenge is thus to factorize (i) the Hamiltonian matrix elements expressed in a quasi-particle basis and (ii) the tensor networks associated to the many-body method and truncation level of interest. Regarding point (i), nothing guarantees that the performance of the TF observed in Ref. Tichai:2018eem extends to matrix elements represented in a symmetry-breaking, i.e. Bogoliubov quasi-particle, basis. Regarding point (ii), tests are presently performed using the simplest of methods, i.e. second-order BMBPT, which is thus tensor decomposed and benchmarked against its full-fledged implementation. Extensions to BCC, PBMBPT, PBCC and GSCGF are left to future works.

While TF numerically filters high-rank components of the many-body tensors, another option consists of discarding a large set of their entries on the basis of a robust but inexpensive a priori estimate on their importance. Iterative equations are then only solved for the dominant entries while an a posteriori correction may be considered for the omitted part of the many-body tensors. As for TF, IT introduces a systematic approximation error that needs to be characterized through benchmark calculations. So far, IT techniques have been mainly used in the context of the no-core shell-model (NCSM) that allows for brute-force solutions of the stationary Schrödinger equation up to mass number $A\approx 12$. As a remedy to the exponential growth of the $A$-body Hilbert space, the IT-NSCM based on an a priori estimate of the relevance of each many-body basis state has allowed to extend the application of diagonalization techniques to $A\approx 25$ in the last decade Roth:2007sv ; Roth09 ; Roth:2008qd . While it is also in connection with configuration interaction (CI) that IT techniques have been mostly developed in quantum chemistry buenker74a ; buenker75a ; illas91a , it can also be combined with Monte Carlo calculations giner13a . As for non-perturbative expansion techniques, e.g (B)CC, (G)SCGF and IMSRG methods, this has been seldom used. Very recently, an interesting work based on CC theory bearing great resemblance with the presently developed idea, has been realized in quantum chemistry deustua17a . In this scheme, a large part of the triple CC amplitude is omitted on the basis of a prior full configuration interaction quantum Monte Carlo (FCIQMC) booth09a ; booth10a or CC Monte Carlo (CCMC) thom10a ; spencer16a ; scott17a calculation.

The paper is organized as follows. Section 2 introduces Bogoliubov coupled cluster theory as the many-body formalism of reference employed throughout the paper to formally illustrate the potential merit of TF and IT techniques. Its perturbative reduction, i.e. Bogoliubov many-body perturbation theory, is further characterized both in view of performing numerical benchmarks and as the inexpensive auxiliary method providing the IT measure. Section 3 provides an overview of TF techniques before displaying the results of their application in small model spaces to quasi-particle matrix elements of the nuclear Hamiltonian. In Sec. 4, the IT concept is introduced and applied in small model spaces to be compared to TF techniques. Next, IT results in large model spaces appropriate for medium-mass systems are displayed and analyzed in details. Section 6 concludes with a set of envisioned extensions for future research.

2 Many-body formalism

The present paper wishes to exemplify the use of TF and IT techniques within the frame of non-perturbative Bogoliubov coupled-cluster theory Si15 . The latter thus serves as a baseline for the discussion although actual BCC calculations based on TF and IT are postponed to a future work. In the present paper, the perturbative reduction of BCC, i.e. BMBPT Tichai:2018mll ; Arthuis:2018yoo , is employed to perform benchmark calculations.

2.1 Hamiltonian tensors

Bogoliubov CC theory and BMBPT expand the exact many-body ground-state around a particle-number breaking Bogoliubov state $|\Phi\rangle$. The reference state $|\Phi\rangle$ is a vacuum for a complete set of quasi-particle operators obtained from particle ones through a unitary Bogoliubov transformation

[TABLE]

such that $\beta_{k}|\Phi\rangle=0$ for all $k$. The columns of the transformation matrices $(U,V)$ are typically obtained as eigenvectors of the HFB equation RiSc80 that also delivers the set of quasi-particle energies $\{E_{k}>0\}$ as eigenvalues. Since the Bogoliubov reference state is not an eigenstate of the particle-number111In practice the constraint has to be done for neutron and proton-number operators $N$ and $Z$, respectively by introducing two separate chemical potentials $\lambda_{N}$ and $\lambda_{Z}$. In our formalism $A$ stands for either one of them. operator $A$, the HFB equation is obtained while constraining the expectation value of $A$ in $|\Phi\rangle$ to match the physical particle number via the use of a Lagrange term, i.e. by minimizing the grand potential

[TABLE]

The many-body methods of interest are most naturally expressed in terms of normal-ordered contributions to the grand potential $\Omega$ with respect to the Bogoliubov reference state. Starting from the Hamiltonian expressed in the one-body basis (Eq. (1)), one obtains

[TABLE]

where $\Omega^{ij}$ denotes the normal-ordered component involving $i$ ($j$) quasi-particle creation (annihilation) operators, e.g.,

[TABLE]

It is clear from Eqs. (5) and (6) that the Hamiltonian collects now a set of mode-$(i\!+\!j)$ tensors $\Omega^{ij}_{k_{1}\ldots k_{i}k_{i+1}\ldots k_{i+j}}$ that are expressed in a symmetry breaking, i.e. quasi-particle, basis. These tensors display antisymmetry properties, i.e.

[TABLE]

where $\sigma(P)$ refers to the signature of the permutation $P$. The notation $P(\ldots|\ldots)$ denotes a separation into the $i$ quasiparticle creation operators and the $j$ quasiparticle annihilation operators such that permutations are only considered among members of the same group. Details on the normal-ordering procedure as well as expressions of the various contributions to $\Omega^{[0]}$, $\Omega^{[2]}$ and $\Omega^{[4]}$ in terms of the original matrix elements of $H$ and of the $(U,V)$ matrices can be found in Ref. Si15 . When the Bogoliubov reference state is chosen to solve the HFB equations, one has

[TABLE]

To circumvent the explicit treatment of three-body operators, i.e., mode-$6$ tensors, state-of-the-art many-body calculations employ the so-called normal-ordered two-body approximation (NO2B). In large-scale NCSM calculations, the error induced by the NO2B approximation to the Hamiltonian was estimated to be of the order of $1$-$3\%$ RoBi12 ; Geb16 up to the oxygen region. While straightforwardly defined in symmetry-conserving methods, the design of a symmetry-conserving approximation of the Hamiltonian in methods based on a symmetry-breaking reference state is non trivial ripoche19a . Eventually, the approximation leads to omitting $\Omega^{[6]}$ (but not only) such that the dominant contribution of the original three-body interaction is included into the retained terms $\Omega^{[0]}$, $\Omega^{[2]}$ and $\Omega^{[4]}$ via the normal-ordering procedure222While the terms $\Omega^{[0]}$, $\Omega^{[2]}$ and $\Omega^{[4]}$ are in fact also affected by the approximation ripoche19a , the original notation is kept for simplicity in the remainder of the present paper.. The present analysis is performed in this context and can later be extended to the use of $\Omega^{[6]}$.

2.2 Many-body tensors

Bogoliubov CC theory relies on the use of the ground-state wave-function ansatz Si15

[TABLE]

where the quasiparticle cluster operator $\mathcal{T}\equiv\mathcal{T}_{1}+\mathcal{T}_{2}+\mathcal{T}_{3}+\ldots$ is defined through

[TABLE]

etc. The BCC amplitudes $t^{m0}_{k_{1}\ldots k_{m}}$, which need to be determined, constitute the fully antisymmetric many-body tensors of present interest.

In the NO2B approximation, the BCC ground-state energy reads as

[TABLE]

where the similarity-transformed grand potential is defined through $\bar{\Omega}\equiv e^{-\mathcal{T}}\Omega e^{\mathcal{T}}$. Similarly to Eq. (2) for standard CC, the BCC correlation energy $\Delta\Omega_{0}\equiv{\cal E}_{0}-\Omega^{00}$ is an explicit function of the sole single and double BCC amplitudes. The determination of the BCC amplitudes defining the connected cluster operator $\mathcal{T}$ relies on solving a set of coupled non-linear equations given in a compact form by

[TABLE]

where even quasi-particle excitations of the vacuum are defined through

[TABLE]

The standard hierarchy of truncation schemes consists in solving the set of equations for $\{\mathcal{T}_{1},\ldots,\mathcal{T}_{r}\}$, i.e. solving Eq. (12) for $p\leq r$, while setting $\mathcal{T}_{q}=0$ for $q>r$. For example, BCC with singles and doubles (BCCSD) retains $\mathcal{T}_{1}$ and $\mathcal{T}_{2}$ while setting $\mathcal{T}_{q}=0$ for $q>2$. This corresponds to solving Eq. (12) in the subspace of Fock space spanned by two and four quasi-particle excitations and denoted as ${\cal F}^{\text{SD}}$ while ignoring $\mathcal{T}_{q}$ with $q>2$. The explicit tensor network associated with Eq. 12 and used to determine $t^{20}_{k_{1}k_{2}}$ and $t^{40}_{k_{1}k_{2}k_{3}k_{4}}$ within BCCSD can be found in Ref. Si15 . The next approximation level, coined as BCCSDT, consists of further considering triples, i.e. the mode-6 tensor $t^{60}_{k_{1}k_{2}k_{3}k_{4}k_{5}k_{6}}$, which impacts the determination of $\mathcal{T}_{1}$ and $\mathcal{T}_{2}$, and thus ultimately the energy. Thus, BCCSDT consists of solving Eq. (12) in the larger subspace of Fock space further spanned by six quasi-particle excitations and denoted as ${\cal F}^{\text{SDT}}$.

It is possible to move to a perturbative version of BCC, i.e. to employ BMBPT Tichai:2018mll ; Arthuis:2018yoo . It consists of bypassing the iterative solving of the amplitudes equations by relying on their perturbative approximations. In the NO2B approximation, cluster amplitudes are given to first order in perturbation by333Using a HFB reference state leads to $t^{20(1)}_{k_{1}k_{2}}=0$ given that $\Omega^{20}_{k_{1}k_{2}}=0$ in this case. The first non-zero contribution to $t^{20}_{k_{1}k_{2}}$ arises in this case at second order.

[TABLE]

where

[TABLE]

The two non-zero contributions in Eq. (14) correspond to the Hugenholtz diagrams displayed in Fig. 2 Si15 . Inserting these expressions in Eq. (2.2) while omitting the last term provides the second-order BMBPT correlation energy $\Delta\Omega_{0}^{(2)}$.

2.3 $J$-coupled scheme

Present calculations are performed while imposing spherical symmetry. One starts from the $N$ lowest one-body eigenstates of the spherical harmonic oscillator (HO) whose associated quantum numbers are

[TABLE]

where $n_{k}$ denotes the radial HO quantum number, $l_{k}$ the orbital angular-momentum quantum number, $j_{k}$ the total angular-momentum quantum number, $m_{k}$ the angular momentum projection and $t_{k}$ the isospin projection. Proceeding to the $J$-coupling of cross-coupled Tichai19unp grand potential and BCC amplitude matrix elements, one obtains the $J$-coupled form of the mode-4 tensors, e.g.

[TABLE]

where each index denotes now the reduced set of quantum numbers

[TABLE]

differing from $k$ by the removal of the magnetic quantum number $m_{k}$. Correspondingly, one initial mode-$n$ tensor444For mode-6 tensors the $J$-coupling is such that the first and second (fourth and fifth) individual angular momenta are coupled to a given $J_{12}$ ($J_{45}$), which is further coupled with the third (sixth) individual angular momentum to a good overall diagonal $J$. See App. A.4 for details. generates a set of $J$-coupled tensors associated with all possible values of the recoupled $(n/2)$-body angular momentum $J$.

Starting from the original tensor indices $\{k_{i}\}$ whose range $N$ is governed by the truncation of the HO one-body basis characterized by the maximum value of the quantum number $e_{\text{max}}=2n_{k}+l_{k}$, the range of the indices $\{\tilde{k}_{i}\}$ of the $J$-coupled tensors is reduced to $\tilde{N}$.

Correspondingly, the tensor networks at play in the many-body method of interest are processed via angular momentum coupling techniques, e.g. the second term in the BCC correlation energy (Eq. 2.2) is re-expressed as

[TABLE]

where $\hat{x}\equiv\sqrt{2x+1}$ VaMo88 . Similarly, the second-order BMBPT correlation energy based on a HFB reference state reduces to

[TABLE]

where the $J$-coupled form of the first-order approximation to the double amplitudes (Eq. 14b) reads as555Working with a spherically-restricted Bogoliubov state, HFB quasi-particle energies are $m_{k}$-independent, i.e., $E_{\tilde{k}}=E_{k}$.

[TABLE]

2.4 Dimensionalities

In order to anticipate the benefit of pre-processing the solving of the $A$-body Schrödinger equation, let us now briefly discuss typical dimensionalities and memory requirements. Table 1 provides the numbers of entries and associated memory of $J$-coupled mode-4 and mode-6 tensors expressed in the $U(1)$-breaking quasi-particle basis as a function of the truncation parameter $e_{\text{max}}$ of the initial spherical HO basis. To this truncation parameter correspond the basis dimension $N$ of the spherical HO basis and the reduced range $\tilde{N}$ of the indices actually labelling the $J$-coupled tensors. Numbers quoted in Tab. 1 correspond to the full mode-4 (mode-6) tensor, i.e. they sum the contributions from all $J$ ($J_{12},J_{45},J$) blocks generated through the angular momentum recoupling. For further orientation, Tab. 2 provides the number of those summed blocks as a function of $e_{\text{max}}$. The mode-4 (mode-6) tensor under consideration in Tab. 1 is representative of $\Omega^{40}$, $\Omega^{04}$ or ${\cal T}_{2}$ ($\Omega^{60}$, $\Omega^{06}$ or ${\cal T}_{3}$) but numbers would be similar for the other contributions to $\Omega^{[4]}$ ($\Omega^{[6]}$), e.g. $\Omega^{22}$ ($\Omega^{33}$).

Two storage schemes are considered. First, $J$-coupled tensors are naively stored with all $\{\tilde{k}_{i}\}$ indices running independently over the $\tilde{N}$ possible sets of reduced quantum numbers. Second, all symmetries of the $J$-coupled tensors, i.e. antisymmetry under the exchange of pair of indices, parity, isospin666The coupling to two-body isospin cannot be exploited for quasi-particle matrix elements. The only constraint left from isospin symmetry is that an even number of indices must carry a neutron or a proton label. and triangular inequalities associated with the (successive) recoupling of pairs of angular momenta, are exploited to avoid storing a large set of null entries. This defines the optimal storage scheme.

While proof-of-principle calculations using $e_{\text{max}}=4$ are performed in Secs. 3 and 4.4, converged calculations of mid-mass nuclei with $A\sim 40-80$ whose results are discussed in Sec. 4.5 typically require $e_{\text{max}}=12$ when using a Hamiltonian softened via a Similarity Renormalization Group (SRG) transformation as described below. Even higher values of $e_{\text{max}}$ are necessary to employ Hamiltonians that are not processed via SRG and/or to compute heavier nuclei. While today’s high performance computers may allow up to 3TB of RAM for a designated memory node and the use of MPI parallelisation may allow to go beyond that, as a rule of thumb, one may consider useful to work with tensors requiring less than 200GB of storage in production calculations777Non-perturbative methods typically require to store several copies of the same tensor produced through successive iterations.. Table 1 demonstrates that, while an optimal storage scheme makes the handling of $J$-coupled mode-4 tensors not problematic for $e_{\text{max}}=12$, the use of $J$-coupled mode-6 tensors, e.g. the handling of triples in BCC theory, is already a tremendous task. Everything becomes all the more challenging when working in $m$-scheme to, e.g., further authorize rotational symmetry to address doubly open-shell nuclei. The above numbers justify to investigate systematic methods to bypass the full handling of large tensors for a (hopefully) negligeable loss of accuracy.

2.5 Hamiltonian

The nuclear Hamiltonian (Eq. (1)) employed in this work has been derived within the frame of chiral effective field theory We90 ; We91 ; Ep09 . It combines a chiral two-nucleon interaction at next-to-next-to-next-to leading order with a cutoff of $\Lambda_{2N}=500\,\text{MeV}$ EnMa03 with a three-nucleon interaction at next-to-next-to leading order with a local regulator based on a cutoff of $\Lambda_{3N}=400\,\text{MeV}$ Na07 ; Roth:2011vt .

The Hamiltonian is further softened using a SRG transformation with a flow parameter $\alpha=0.08\,\text{fm}^{4}$ BoFu07 ; HeRo07 ; RoRe08 ; RoLa11 ; JuMa13 . This transformation induces many-nucleon forces that are included consistently up to the three-nucleon level, i.e., chiral and induced many-body forces beyond that level are neglected. SRG-evolved Hamiltonians have already been used in a number of medium-mass calculations and have been shown to be soft enough to be used meaningfully in MBPT Tichai:2016joa ; Tichai:2017rqe ; Hu16 and BMBPT Tichai:2018mll calculations.

3 Tensor factorization

In Ref. Tichai:2018eem , TF techniques have been used for the first time in the context of ab initio nuclear structure calculations. In such an approach, tensors are decomposed and approximated on the basis of a particular tensor format.

3.1 Basic considerations

A mode-$n$ tensor $T$ is a multivariate data array $T_{i_{1}...i_{n}}$ with index ranges $\{I_{1},...,I_{n}\}$ and can be seen as a higher-mode analogue of vectors and matrices. In order to quantify the distance between tensors, the Frobenius norm of a tensor $T$ is introduced as

[TABLE]

Given an approximation $\hat{T}$ of $T$, the relative error is thus defined as

[TABLE]

3.2 Tensor hypercontraction

The key idea behind a given approximation is the decomposition of a mode-$n$ tensor into a sum of products of mode-$k$ tensors with $k<n$. The particular topology of the decomposition is referred to as a tensor format. Over the past years various tensor formats have been successfully applied in atomic and solid-state physics as well as in quantum chemistry, e.g., canonical polyadic decomposition (CPD), resolution of identity (RI) or tensor hypercontraction (THC).

Tensor hypercontraction Ho12a ; Ho12b ; Schu17 consists in expanding a mode-4 tensor under the particular form

[TABLE]

Every tensor format comes with a tensor rank that is defined as the size of the auxiliary index used when truncating the decomposition. In the present case, the THC rank $r_{\text{THC}}$ is taken to be the same for both summations over $\alpha$ and $\beta$ in Eq. (24). Increasing the tensor rank of a given decomposition lowers the error such that a systematically improvable approximation is obtained. The THC decomposition involves five factor matrices: four factors $\{X^{i}\}$ whose size scales linearly with $r_{\text{THC}}$ and the core tensor $W$ that scales quadratically with $r_{\text{THC}}$.

The THC format is very flexible and encompasses the CPD or the purely separable ansatz as particular cases, i.e. setting $W_{\alpha\beta}=\delta_{\alpha\beta}$ the CPD format is recovered as

[TABLE]

whereas further assuming $r_{\text{THC}}=1$ provides the purely separable approximation

[TABLE]

Of importance is the data compression ratio

[TABLE]

which relates the initial amount of data in the optimal storage scheme to the compressed amount of data after the decomposition process and the truncation have been achieved. Whenever $R_{c}>1$, the compressed tensor $\hat{T}$ requires less storage than the original one.

3.3 Factorized grand potential

The algorithmic procedure to determine the THC factor matrices was laid out in Ref. Tichai:2018eem . A numerical analysis of each of the individual steps was provided through the application of THC to state-of-the-art 2N+3N nuclear Hamiltonians represented in various symmetry-conserving bases. The same procedure is presently employed to decompose each component888It is only necessary to investigate ${}^{J}{\Omega^{22}}$, ${}^{J}{\Omega^{31}}$ and ${}^{J}{\Omega^{40}}$ given that matrix elements of ${}^{J}{\Omega^{13}}$ and ${}^{J}{\Omega^{04}}$ are trivially related to those of ${}^{J}{\Omega^{31}}$ and ${}^{J}{\Omega^{40}}$ Si15 . of $\Omega^{[4]}$ in $J$-coupled scheme.

The first step of the THC decomposition involves a (truncated) singular value decomposition (SVD). Grouping the first and second pair of indices999Given that quasi-particle matrix elements do not necessarily have the same number of creator and annihilation operators, there is no a priori natural coupling between the first two and the last two labels. However, the step leading to cross-coupled matrix elements Tichai19unp does provide a natural coupling, which is indeed used here., the mode-4 tensors are rewritten as

[TABLE]

where $\tilde{K}\equiv(\tilde{k}_{1},\tilde{k}_{2})$ and $\tilde{K}^{\prime}\equiv(\tilde{k}_{3},\tilde{k}_{4})$ are two-body indices of range $\tilde{N}^{2}$. The SVD is applied to ${}^{J}\tilde{\Omega}^{40}_{\tilde{K}\tilde{K}^{\prime}},^{J}\tilde{\Omega}^{31}_{\tilde{K}\tilde{K}^{\prime}}$ and ${}^{J}\tilde{\Omega}^{22}_{\tilde{K}\tilde{K}^{\prime}}$ and the absolute sizes of the (ordered) sequence of singular values are displayed in Fig. 3. Calculations are performed in 18O using $e_{\text{max}}=4$. As visible in Tab. 1, the corresponding one-body index range is $N=140$ for the original $m$-scheme matrix elements and $\tilde{N}=30$ for the $J$-coupled matrix elements under present investigation.

Singular values behave differently for the three grand potential components. The most efficient truncation can be achieved for ${}^{J}\tilde{\Omega}^{40}$ for which only 20% of singular values need to be retained even for intermediate $J$ values. In the most optimal $J$ blocks, up to 95% of the singular values can be safely discarded. The situation is slightly worse for ${}^{J}\tilde{\Omega}^{31}$ and ${}^{J}\tilde{\Omega}^{22}$ that require keeping up to 35% of the singular values, although large $J$ values, again, authorize more severe truncations. As Fig. 4 testifies, these results directly correlate with the sparsity of each tensor ${{}^{J}}\tilde{\Omega}^{ij}_{\tilde{K}\tilde{K}^{\prime}}$, i.e. the smaller the number of initial non-zero entries, the faster the decrease of the singular values.

While not shown here, similar conclusions have been reached for other semi-magic nuclei. Furthermore, results are qualitatively and quantitatively similar to those obtained for matrix elements of the Hamiltonian represented in the symmetry-conserving HF basis in closed-shell nuclei Tichai:2018eem . As a matter of fact, singular values presently extracted for ${}^{J}\tilde{\Omega}^{40}$ decrease even faster. Inspecting the expression of $\Omega^{ij}_{k_{1}k_{2}k_{3}k_{4}}$ in terms of the Bogoliubov matrices $(U,V)$ Si15 , such a positive outcome could not be easily anticipated. Indeed, the grouping of quasiparticle indices in ${{}^{J}}\tilde{\Omega}^{ij}_{\tilde{K}\tilde{K}^{\prime}}$ relates to a mixture of so-called natural and un-natural grouping of single-particle indices in the original matrix elements Tichai:2018eem . This feature directly reflects the presence of pairing correlations. Knowing that the un-natural grouping of single-particle indices was shown to lead to very slowly decreasing singular values and inefficient truncations Tichai:2018eem , the results displayed in Fig. 3 constitute a non-trivial result that need to be further investigated in heavier nuclei that are more strongly paired than the oxygen isotopes.

Figure 5 shows the THC error as a function of the THC rank $r_{\text{THC}}\equiv\tilde{N}^{\alpha}$. The THC rank is expressed as a certain power $\alpha$ of the range of the tensor indices to emphasize the relation to the basis dimension. The behavior of the decomposition error as a function of the two-body angular-momentum $J$ is the same for the three grand potential components, i.e., for high values of $J$ the approximation error drops more rapidly than for lower values of $J$. For all $J$ values, a THC approximation error $\Delta^{J}{\tilde{\Omega}^{ij}}=10^{-1}$ can be obtained for $r_{\text{THC}}=\tilde{N}^{1.6}$ for ${}^{J}{\tilde{\Omega}^{40}}$ and for $r_{\text{THC}}=\tilde{N}^{1.8}$ for ${}^{J}{\tilde{\Omega}^{31,22}}$, thus, allowing for a robust extension of the THC ansatz to quasi-particle matrix elements.

To illustrate the memory gain obtained via the application of THC, the data compression factor $R_{C}$ is displayed in the left panel of Fig. 6 for 18,20O. The approximation on the overall $\Omega^{[4]}$ tensor is set by using an error threshold $\epsilon$ common to all $J$ channels. This corresponds to using $J$- and component-dependent THC ranks obtained by setting $\Delta^{J}{\tilde{\Omega}^{ij}}=\epsilon$ for all $(i,j)$ and $J$. While $R_{C}=1$ in the limit where no approximation is made on the grand potential tensors, the compression factor increases with decreasing accuracy, i.e. increasing $\epsilon$. The trend is monotonous and identical for both nuclei under investigation. Eventually, setting $\epsilon=10^{-1}$ authorizes to compress the data by about a factor of 10. While this result is achieved in a small model space, much larger compression factor (for the same error threshold) are envisioned for larger model spaces appropriate to realistic ab initio calculations.

3.4 Factorized many-body tensor and energy

Eventually, the goal is to decompose both Hamiltonian and many-body tensors at play, eventually leading to a tensor-factorized form of the many-body formalism Schu17 . This is envisioned to be truly beneficial in non-perturbative methods where the iterative equations can be solved for the low-rank factors in a specific tensor format rather than for the original many-body tensors. While this is our goal to do so for, e.g., BCC in the future, second-order BMBPT introduced in Secs. 2.2 and 2.3 is presently employed to gauge the propagated error to the ground-state energy.

The decomposition of the tensor network associated with $\Delta\Omega^{(2)}_{0}$ can be achieved by inserting the decomposition of ${}^{J}{\tilde{\Omega}^{04}_{\tilde{k}_{1}\tilde{k}_{2}\tilde{k}_{3}\tilde{k}_{4}}}$ and ${}^{J}{\tilde{t}^{40(1)}_{\tilde{k}_{1}\tilde{k}_{2}\tilde{k}_{3}\tilde{k}_{4}}}$ into Eq. (20). The THC decomposition of the grand potential tensor is already at hands under the form

[TABLE]

A decomposition101010An alternative would be to perform the numerical THC decomposition of ${}^{J}{\tilde{t}^{04(1)}_{\tilde{k}_{1}\tilde{k}_{2}\tilde{k}_{3}\tilde{k}_{4}}}$ itself rather than combining the decompositions of the tensors that compose it as is presently done. In an actual non-perturbative BCC calculation, an ansatz for the factorization of the cluster amplitudes would be made a priori. of ${}^{J}{\tilde{t}^{40(1)}_{\tilde{k}_{1}\tilde{k}_{2}\tilde{k}_{3}\tilde{k}_{4}}}$ (Eq. 21) is achieved by combining Eq. (29) with the inverse Laplace transform of the energy factor in Eq. 21

[TABLE]

where a numerical quadrature has been used. Equation (30) corresponds to a CPD whose rank is $r_{D}\equiv 2M+1$. Numerical values of $M$, $\omega_{s}$ and $t_{s}$ are tabulated in the literature BrHa05 such that high precision can be reached with very few grid points. In particular, the size of the integration mesh was shown to be independent of the system. The pairing gap associated with the spontaneous breaking of particle-number symmetry makes all energy denominators non-vanishing and, therefore, the quadrature well-defined111111If $U(1)$ symmetry does not spontaneously break, HFB quasi-particle energies reduce to HF single-particle energies with a significant shell gap near the Fermi energy, i.e. a significant HOMO-LUMO gap in the quantum-chemistry language. Consequently, the quadrature is also well defined in this particular limit.. Combining the above elements, the factorized form of the second-order BMBPT energy correction is obtained as

[TABLE]

where the intermediates

[TABLE]

have been introduced. The evaluation cost of these intermediates is $\mathcal{O}(r^{2}_{\text{THC}}Nr_{D})$. Defining new intermediates via

[TABLE]

one is eventually left with

[TABLE]

The computational scaling of the factorized $\Delta\Omega^{(2)}_{0}$ depends on the THC ranks $r_{\text{THC}}$ in each $J$ block, which themselves depend on the chosen approximation error $\epsilon$ on the grand potential. As discussed in Ref. Tichai:2018eem , the scaling is typically worse than the naive $\mathcal{O}(N^{4})$ scaling of the original second-order correction. However, the aim of the present work is not to derive a low-scaling approximation of the already low-cost BMBPT(2) but rather to benchmark the propagation of the TF error to nuclear observables. One has to move to many-body methods that are more expensive to begin with to generate a reduction of the numerical scaling.

To measure the impact of TF, the relative error on the second-order energy correction is introduced

[TABLE]

which goes to zero in the limit of an exact THC decomposition. Because of the highly-accurate Laplace transformed introduced above, the error presently propagates entirely from the approximation made on ${}^{J}{\tilde{\Omega}^{04}}$.

The right panel of Fig. 6 displays $\Delta(\Delta\Omega^{(2)}_{0})$ for 18,20O as a function of $\epsilon$. A global trend is visible such that lower values of $\epsilon$ yield lower $\Delta(\Delta\Omega^{(2)}_{0})$. One observes that the behavior is identical in both nuclei. Eventually, a THC approximation error of $\epsilon\approx 10^{-1}$ is sufficient to obtain $\Delta(\Delta\Omega^{(2)}_{0})$ below $1\%$. Therefore, even though the matrix elements are only approximated to an accuracy of $10^{-1}$, the precision on the observable of interest is more than one order of magnitude better, which is achieved by retaining one order of magnitude less entries than in the original tensor. Even a quite crude approximation on the matrix elements thus yields an accuracy that is good enough to perform precision studies.

4 Importance truncation

As an alternative to discarding high-rank components of the many-body tensors via a numerical factorization, a procedure to remove entries on the basis of an importance measure is now investigated.

4.1 Basic concept

The general idea is to discard irrelevant entries of the largest mode-$n$ tensors at play in the many-body framework of interest. This is done by a priori estimating the importance of each of its entries on the basis of a less costly many-body method than the envisioned one. In order to illustrate the concept, it is necessary to slightly reformulate the BCC method introduced in Sec. 2.2.

One considers so-called $P$ and $Q$ subspaces ${\cal F}^{(P)}$ and ${\cal F}^{(Q)}$ of Fock space ${\cal F}$ spanned by two selected sets of quasi-particle excitations of the Bogoliubov vacuum $|\Phi\rangle$ such that ${\cal F}^{(Q)}\subseteq({\cal F}^{(0)}+{\cal F}^{(P)})^{\perp}$, where ${\cal F}^{(0)}$ is the one-dimensional subspace spanned by $|\Phi\rangle$. Typically, ${\cal F}^{(P)}$ is spanned by a selected set of low-rank excitations $\{|\Phi^{K}\rangle\}$ of $|\Phi\rangle$. Considering the cluster operator ${\cal T}^{(P)}$ associated with those excitations and solving the set of P-space BCC amplitude equations

[TABLE]

where $\bar{\Omega}^{(P)}\equiv e^{-\mathcal{T}^{(P)}}\Omega e^{\mathcal{T}^{(P)}}$, one obtains a first approximation to the ground-state energy via

[TABLE]

Once this is done, a non-iterative correction associated with quasi-particle excitations in ${\cal F}^{(Q)}$ is computed such that the energy eventually reads a

[TABLE]

Several variants of (B)CC theories fit into the above formulation. For instance, setting ${\cal T}^{(P)}\equiv{\cal T}_{1}+\ldots{\cal T}_{r}$ and taking a null $Q$ space, standard BCC truncations with $\delta(P;Q)=0$ are recovered. Choosing the same $P$ space but computing $\delta(P;Q)$ in larger and larger $Q$ spaces via the method of moments, the Bogoliubov extension of the completely renormalized CC hierarchy piecuch05a ; piecuch06a ; shen12a is obtained. Using a more flexible definition of ${\cal F}^{(P)}$ while still computing $\delta(P;Q)$ via the method of moments, the Bogoliubov extension of the CC(P;Q) hierarchy shen12a ; shen12b ; shen12c ; bauman17a is additionally obtained.

In the present paper, we wish to exploit the above formulation in a spirit close to BCC(P;Q) except that the energy correction $\delta(P;Q)$ is meant to be computed in perturbation and not via the more advanced method of moments121212This choice is made for simplicity and using the method of moments is envisioned for future applications..

4.2 Importance measure and ${\cal F}^{(P)}$ selection

The general notion of IT has been developed in many different forms in the past. It already concerns the very initial truncation of the many-body basis or of the one-body basis used to represent the many-body tensors at play.

In a conventional NCSM framework for example, one builds the $A$-body configuration space by including all particle-hole excitations below $N_{\text{max}}\hbar\Omega$, where $N_{\text{max}}$ denotes the sum of individual harmonic-oscillator excitation quanta with respect to the reference configuration. It corresponds to using $A$-body unperturbed excitation energies as an importance measure for the $A$-body basis states.

In the present context, a similar truncation is at play when originally representing the $k$-body (mode-$2k$) operators (tensors) in the HO basis while employing a so-called "$e_{k\text{max}}$ truncation" to select the entries on the basis of the associated $k$-body HO energies. While one originally envisions to use a consistent scheme131313This consistency is necessary for the truncated object initially representing the $k$-body operator to be an actual mode-$2k$ tensor. for all tensors involved by using $e_{k\text{max}}=ke_{\text{max}}$, further truncations are typically envisioned for the largest $k$-body operators by relaxing the consistency between the various $k$-body sectors through the actual use of $e_{k\text{max}}<ke_{\text{max}}$.

Given the initial representation of the tensors, the key point relates to the further selection of ${\cal F}^{(P)}$ and ${\cal F}^{(Q)}$. Typically, the goal is to approximate a full BCCSD or BCCSDT calculation at a (much) reduced computational cost. As such, one chooses ${\cal F}^{(P)}+{\cal F}^{(Q)}={\cal F}^{\text{SD}}$ or ${\cal F}^{\text{SDT}}$ and exploit the flexibility in the partitioning between $P$ and $Q$ to only solve the non-iterative equations in an optimal subspace ${\cal F}^{(P)}$ given a targeted accuracy. Aiming at BCCSD (BCCSDT), this translates into the fact that only a small subset of the original $t^{40}_{k_{1}k_{2}k_{3}k_{4}}$ ($t^{60}_{k_{1}k_{2}k_{3}k_{4}k_{5}k_{6}}$) tensor entries141414Given that treating $\mathcal{T}_{1}$ in full is doable even for large $N$, single excitations are included in ${\cal F}^{(P)}$ by default. is carried along when solving Eq. (35) iteratively.

To actually select ${\cal F}^{(P)}$, i.e. the subset of $t^{40}_{k_{1}k_{2}k_{3}k_{4}}$ ($t^{60}_{k_{1}k_{2}k_{3}k_{4}k_{5}k_{6}}$) entries to be solved for, the simplest idea is to rely on a zeroth-order-like measure

[TABLE]

which corresponds to selecting the tensor entries according to the sum of unperturbed quasi-particle energies associated with its $2k$ indices. While computationally inexpensive, such a scheme does neither take into account properties of the full Hamiltonian. Anticipating the too crude character of this truncation, one is led to evaluating the importance of the entries on the basis of a more optimal measure. Generically speaking, the relevance of a given importance measure relies on the following considerations

the measure must provide a robust estimate of the importance of the tensor entries delivered by the targeted calculation, 2. 2.

the evaluation of the importance measure must be as computationally inexpensive as possible, i.e. it must offer a significant gain compared to the evaluation of the tensor entries in the targeted calculation.

In practice, a tradeoff between both requirements must be found. While the most precise estimate is obtained by evaluating the tensors of interest from a large-scale full configuration interaction calculation, it is of no practical interest since it already requires the full solution of the quantum many-body body problem. On the other hand a too simplistic estimate, e.g., through $\kappa^{(0)}$ introduced in Eq. (38), might not be related strongly enough to the full-fledged solution of interest.

In the present work, the lowest-order BMBPT counterpart of the tensors of interest is used as an importance measure. Alternatives, e.g. the early stages of FCIQMC or CCMC propagations deustua17a , can be envisioned if necessary. The first-order BMBPT estimates $\kappa^{(1)}$ of the double and triple amplitudes based a HFB vacuum and the NO2B approximation read as

[TABLE]

The merit of $\kappa^{(1)}(t^{40}_{k_{1}k_{2}k_{3}k_{4}})$ over $\kappa^{(0)}(t^{40}_{k_{1}k_{2}k_{3}k_{4}})$ will be exemplified in realistic calculations below.

While first-order BMBPT offers a $N^{4}$, i.e. a low cost, estimate of the double amplitudes, it does not do so for the triple amplitudes that are strictly zero at that order151515In presence of $\Omega^{[6]}$, the leading triple amplitudes do arise at first order with a $N^{6}$ complexity.. A non-zero estimate of the triple amplitudes requires to go to second order and reads as

[TABLE]

where the operator $P(k_{1}k_{2}k_{3}/k_{4}k_{5}k_{6})$ permutes members of the triplet $(k_{1},k_{2},k_{3})$ with those of the triplet $(k_{4},k_{5},k_{6})$, i.e. it generates 20 different terms. The Hugenholtz diagram corresponding to Eq. (40) is displayed in Fig. 7. Its evaluation is a costly non-iterative $N^{7}$ process. While constituting a state-of-the-art many-body calculation in itself, the evaluation of $t_{k_{1}k_{2}k_{3}k_{4}k_{5}k_{6}}^{60(2)}$ remains very advantageous compared to the even more challenging iterative $N^{8}$ BCCSDT calculation of the $\mathcal{T}_{3}$ tensor entries. Full-fledged $J$-scheme calculations of $t_{k_{1}k_{2}k_{3}k_{4}k_{5}k_{6}}^{60(2)}$ in large model spaces will be presented below to both set up the IT of a future BCCSDT calculation and to evaluate the corresponding $\delta(P;Q)$ perturbative correction to the energy. Performing the angular momentum recoupling of the mode-6 tensor $t_{k_{1}k_{2}k_{3}k_{4}k_{5}k_{6}}^{60(2)}$ is also challenging in itself, especially because the recoupling is different for the 20 terms generated by $P(k_{1}k_{2}k_{3}/k_{4}k_{5}k_{6})$ and spans several levels of complexity/computational cost. See App. A.4 for a detailed discussion of this point leading to the introduction of the $J$-coupled tensor

[TABLE]

whose second-order BMBPT estimate will be investigated below.

4.3 Importance-truncated tensor

Considering the $n$-tuple BCC amplitude tensor

[TABLE]

the corresponding importance-truncated tensor based on the IT measure $\kappa^{(p)}(t^{2n0}_{k_{1}\ldots k_{2n}})$ is obtained as

[TABLE]

where $\kappa^{(p)}_{\text{min}}$ defines the IT threshold. The original tensor is obtained in the limit of $\kappa^{(p)}_{\text{min}}\rightarrow 0$, i.e.,

[TABLE]

Additionally, for $\kappa^{(p)}_{\text{min}}<{\kappa^{(p)\prime}_{\text{min}}}$ the subset relation holds

[TABLE]

such that lowering $\kappa^{(p)}_{\text{min}}$ increases monotonically the number of entries in the truncated tensor.

4.4 Small-scale IT applications

The IT concept is first applied to the double BCC amplitude tensor $\mathcal{T}_{2}$ calculated for 18,20O in an $e_{\text{max}}=4$ one-body HO basis, i.e. for $N=140$ and $\tilde{N}=30$. This choice authorizes a direct comparison to the results previously obtained with TF techniques.

The simplest importance measure $\kappa^{(0)}$ is investigated first, with identical results for both nuclei. The left panel of Fig. 8 displays the compression ratio as a function of $\kappa^{(0)}_{\text{min}}$. For $\kappa^{(0)}_{\text{min}}\leq 5\cdot 10^{-2}$ the compression ratio is $R_{c}=1$ such that the IT-truncation is inactive. As expected, the compression ratio decreases monotonically when lowering the IT threshold. Reaching, e.g., a compression ratio $R_{c}=10$ corresponds to setting $\kappa^{(0)}_{\text{min}}\approx 10^{-2}$ and relates to an error $\Delta\mathcal{T}_{2}=0.5$ on the tensor. For comparison, the same compression ratio induced an error $\epsilon=10^{-1}$ on the Hamiltonian tensor $\Omega^{[4]}$ in the THC framework. The right panel of Fig. 8 displays the relative error on the second-order BMBPT correlation energy. For $R_{c}=10$, i.e. $\kappa^{(0)}_{\text{min}}\approx 10^{-2}$, a 50$\%$ error is generated, which compares very negatively to the sub-percent error obtained in the THC for a similar compression factor of the grand potential tensor.

As a next step, the left panel of Fig. 9 displays $R_{c}$ as a function of $\kappa^{(1)}_{\text{min}}$ for the same nuclei. Again, the compression ratio decreases monotonically when lowering the IT threshold. For $\kappa^{(1)}_{\text{min}}=10^{-3}$, a compression ratio $R_{c}\approx 10$ is obtained for both nuclei, which corresponds to an error $\Delta\mathcal{T}_{2}=10^{-2}$ on the tensor. The associated relative error on the second-order BMBPT correlation energy (right panel) is equal to $1\%$. Lowering the IT threshold to $\kappa^{(1)}_{\text{min}}=3\cdot 10^{-4}$ barely reduces $R_{c}$ while decreasing the error by another order of magnitude.

Figure 10 summarizes the results by correlating $\Delta(\Delta\Omega^{(2)}_{0})$ with $R_{c}$ for both TF and IT pre-processing methods. While operating in a very different fashion, TF based on THC and IT based on the first-order BMBPT estimation of $\mathcal{T}_{2}$ provide very similar performances. In particular, maintaining the relative error below $1\%$ can be achieved while compressing the data by one order of magnitude. While already satisfactory, Sec. 4.5 below illustrates how this result extends very favourably to larger model spaces and/or mode-6 tensors. It is also clear from Fig. 10 that employing a too naive IT measure, e.g. $\kappa^{(0)}$ based on unperturbed excitation energies of the associated configurations, is highly inefficient and must be discarded.

4.5 Large-scale IT applications

The conceptual simplicity of IT allows its implementation and application in large model spaces. Furthermore, at the price of a state-of-the-art many-body development, the IT pre-processing is also applied to the mode-6 triple BCC amplitude tensor ${\cal T}_{3}$ in large model spaces.

4.5.1 Double BCC amplitudes

Applying IT to $\mathcal{T}_{2}$, the goal is to solve BCCSD equations for the retained entries (defining ${\cal F}^{(P)}$) and correct for the omitted ones (defining ${\cal F}^{(Q)}$) in perturbation. While solving CCSD in large model spaces is doable in both $J$- and $m$-scheme, BCCSD is more demanding and actually challenging given that the $N^{4}$ storage cost invokes the full basis dimension $N$ that cannot be separated into particle and hole states. In this context, IT can either provide a significant leverage in $J$-scheme161616As seen below, the optimal storage cost of $\mathcal{T}_{2}$ in $J$-scheme is of the order of 200 Mb for $e_{\text{max}}=12$. or a mandatory leverage in $m$-scheme at the BCCSD level171717The storage cost of $\mathcal{T}_{2}$ in $m$-scheme is of the order of 3 Tb for $e_{\text{max}}=12$ Signoracci:2014dia .. While this will actually be investigated in a future publication, the present analysis is limited to characterizing the compression ratio obtained in $J$-scheme as a function of the IT threshold and to evaluating the error induced on the second-order BMBPT correlation energy by straightforwardly omitting the entries below the IT threshold.

The performance of IT techniques applied to the double BCC amplitude tensor $\mathcal{T}_{2}$ in a small model space has been characterized in Sec. 4.4. The calculations are now extended up to $e_{\text{max}}=12$, i.e. up to $N=1820$ and $\tilde{N}=182$. In order to appreciate the data compression achieved for a given accuracy, Fig. 11 displays the number of entries of the $J$-coupled $\mathcal{T}_{2}$ tensor as a function of the size of the one-body HO basis. This is done for the naive and optimal storage schemes of the initial tensor as well as for the IT tensor based on $\kappa^{(1)}$. The numbers displayed for the initial tensor correspond to those given in Tab. 1. As for the IT tensor, results are provided for several values of the relative error on the second-order BMBPT correlation energy in 18O, ranging arbitrarily from $0.01\%$ to $10\%$.

The first observation is that the naive and optimal storage schemes both exhibit exponential growth, however, with different rates. The naive scheme grows much more rapidly, finally requiring about three orders of magnitude more storage than the optimal one in $e_{\text{max}}=12$. The final storage requirement for the optimal scheme is about $0.5\,$Gb. While applying IT permitted to reduce the number of entries by one order of magnitude in an $e_{\text{max}}=4$ model space for a $1\%$ error on the correlation energy, the performance increases tremendously as a function of $e_{\text{max}}$. For $e_{\text{max}}=12$, the data is compressed by three orders of magnitude, six orders of magnitude compared to the naive storage scheme, for the same $1\%$ error on the correlation energy. Even for a very accurate $0.01\%$ error, one still obtains a data compression by more than two orders of magnitude.

In order to better characterize the effect of IT, Fig. 12 compares the number of entries per $J$ block before and after IT. The IT results are displayed for $\kappa^{(1)}_{\text{min}}=9\cdot 10^{-3}$ and $\kappa^{(1)}_{\text{min}}=10^{-5}$, which correspond to $1\%$ and $0.01\%$ relative errors on the second-order energy correction, respectively. While a significant compression is observed for all $J$ values, a clear trend emerges in average: the larger the angular momentum associated to the entries, the lesser their importance. While the initial entries extend up to $J=2e_{\text{max}}+1=25$, there remains no entry beyond $J=8$ in the very precise IT tensor corresponding to $\kappa^{(1)}_{\text{min}}=10^{-5}$. No entries beyond $J=6$ are necessary to reach a $1\%$ error on the second-order energy correction.

In addition to the storage of $\mathcal{T}_{2}$, the CPU runtime is the other critical component of a computational analysis. While being computationally very simple, the evaluation of the second-order energy correction in large model spaces provides a useful testground for a CPU analysis. In order to produce a transparent comparison, the parallelization of the code was turned off for the present analysis.

Figure 13 displays the speedup obtained by IT-BMBPT compared to the exact evaluation of the second-order correction in 18O. Two orders of magnitude in runtime are gained for very accurate calculations corresponding to $\kappa^{(1)}_{\text{min}}=10^{-5}$ and four orders of magnitude are gained at the $1\%$ error level. Whereas the goal is not to obtain a speedup for a low-cost theory such as BMBPT(2), these numbers are very encouraging in view of performing IT-BCCSD calculations in the future.

So far, results have been displayed for weakly-paired Oxygen isotopes181818In prior BMBPT calculations using the same Hamiltonian, oxygen isotopes were identified to be only weakly paired as testified by the small particle-number variance Tichai:2018mll .. We now wish to investigate how the IT technique behaves for heavier nuclei and when going from a closed-shell to an open-shell system . To do so, IT calculations are repeated for the unpaired doubly closed-shell 40Ca nucleus and for the singly open-shell 44Ca isotope that exhibits stronger pairing191919Note however that 44Ca computed with the same Hamiltonian via GSCGF theory at the ADC(2) truncation level, which is close to BMBPT(2) Tichai:2018mll , displays only about half of the experimental pairing gap measured via the three-point mass difference formula Barb18 . than 18,20O.

Figure 14 displays the relative error on the second-order BMBPT ground-state correlation energy against the data compression factor for 18O and 40,44Ca. Comparing first the two paired systems 18O and 44Ca, one observes a mass dependence while working with a fixed basis size ($e_{\text{max}}=12$) that is relatively larger to begin with for 18O. While the two curves follow the same overall trend, the data compression achieved for a given accuracy is larger in 18O, e.g. the $1\%$ accuracy achieved for $R_{c}=2\cdot 10^{3}$ in 18O is only reached for $R_{c}=4\cdot 10^{2}$ in 44Ca. Focusing next on 40Ca and 44Ca, the curves follow qualitatively two different trends. A significantly smaller error is generated in 40Ca than in 44Ca when truncating the smallest entries. Because 40Ca is not superfluid, the associated HFB vacuum reduces to a HF Slater determinant such that the matrix elements of ${}^{J}{\tilde{\Omega}^{40}}$, and thus of ${}^{J}{\tilde{t}^{40(1)}}$, display a particle-hole symmetry and thus more zero entries than for 44Ca to begin with. It eventually leads to a more efficient data compression for 40Ca. These characteristics are visible in Fig. 16. For example, the compression obtained for a $1\%$ error is twice as better in 40Ca than in 44Ca. For errors larger than $10\%$, the benefit associated with the particle-hole symmetry for high accuracy calculations disappears and both curves converge towards one another such that the closed- or open-shell character of the system becomes irrelevant.

4.5.2 Triple BCC amplitudes

Applying IT to $\mathcal{T}_{3}$, the goal is to solve BCCSDT equations for the retained entries (defining ${\cal F}^{(P)}$) and correct for the omitted ones (defining ${\cal F}^{(Q)}$) in perturbation. The feasibility of the approach directly depends on the reduction offered by the IT for a desired accuracy given that a full BCCSDT calculation is currently undoable in realistic model spaces, even in $J$-scheme.

The IT measure requires the perturbative evaluation of $\mathcal{T}_{3}$ as defined by Eq. (40) (see Eq. (A.14) for the $J$-coupled form), which is a $N^{7}$ process. Based on this estimate, the associated contribution to the fourth-order ground-state correlation energy

[TABLE]

which is necessarily positive, is used to gauge the propagated error. The evaluation of $\Delta\Omega^{[4_{T}]}_{0}$ is itself a $N^{6}$ process.

In order to appreciate the data compression achieved for a given accuracy, Fig. 17 displays the number of entries of the $J$-coupled $\mathcal{T}_{3}$ tensor as a function of the size of the one-body basis. This is done for the naive and optimal storage schemes of the initial tensor as well as for the IT tensor based on $\kappa^{(2)}$. The numbers displayed for the initial tensor correspond to those given in Tab. 1. For $e_{\text{max}}=10,12$, additional numbers are given for the optimal storage scheme that correspond to performing an additional reduction to $e_{3\text{max}}=22<3\cdot e_{\text{max}}$. This additional truncation, which still leads to working with an extremely large mode-6 tensor, is necessary to evaluate $t_{k_{1}k_{2}k_{3}k_{4}k_{5}k_{6}}^{60(2)}$. The IT results correspond to a $1\%$ relative error on the second-order BMBPT correlation energy in 18O and are based on the initially reduced tensor for $e_{\text{max}}=10,12$.

Once again, the naive and optimal storage schemes both exhibit exponential growth. The storage requirements in the optimal scheme are about $15\,$Tb in $e_{\text{max}}=12$ for $e_{3\text{max}}=36$ and $200\,$Gb for the reduced $e_{3\text{max}}=22$, which is unmanageable in a BCCSDT calculation. While the IT permits to reduce the number of entries by one order of magnitude at $e_{\text{max}}=4$ for a $1\%$ error on $\Delta\Omega^{[4_{T}]}_{0}$, the performance increases tremendously as a function of $e_{\text{max}}$. For $e_{\text{max}}=12$ and $e_{3\text{max}}=22$, the data is compressed by 2 orders of magnitude202020Starting from $e_{3\text{max}}=36$, the compression factor would undoubtedly be significantly larger than 2 orders of magnitude such that the number of remaining entries would probably not be much larger than in the present calculation. for the same $1\%$ error such that the IT tensor to handle requires less than $1\,$Gb. As illustrated in Tab. 3, it must be noted that $\Delta\Omega^{[4_{T}]}_{0}$ being one order of magnitude smaller than $\Delta\Omega^{[2]}_{0}$, requiring $1\%$ error on the former is equivalent to requiring $0.1\%$ error on the latter, i.e. an error of $30$ keV out of a total binding energy of about $135$ MeV in 18O. This level of error is largely negligible given the predictive power of current ab initio calculations.

In order to better characterize the effect of IT, Fig. 18 compares the number of entries per $J$ block before and after IT. The IT results are displayed for $\kappa^{(2)}_{\text{min}}=10^{-6}$ and $\kappa^{(2)}_{\text{min}}=6\cdot 10^{-6}$, which corresponds to $1\%$ and $5\%$ relative error on $\Delta\Omega^{[4_{T}]}_{0}$, respectively. While a significant compression is observed for all $J$ values, the same trend emerges as for $\mathcal{T}_{3}$, i.e. the larger the three-body angular momentum associated to the entries, the lesser their importance. While the initial entries extend up to $J=43/2$, there remains no entry beyond $J=19/2$ for $\kappa^{(2)}_{\text{min}}=10^{-6}$. No entries beyond $J=17/2$ are necessary to reach a $5\%$ error on $\Delta\Omega^{[4_{T}]}_{0}$.

4.6 Discussion

Let us now make three comments on IT to anticipate further developments

While the entire benefit of the IT cannot be obtained in this way, it does emerge a posteriori that a large part of its effects is to eliminate the entries entering large-$J$ blocks. A systematic study as a function of $e_{\text{max}}$ could allow us to understand quantitatively which two-body and three-body $J$ blocks are entirely irrelevant for a given accuracy on the correlation energy. Based on the results of such a study, one could avoid producing those $J$ blocks from the outset, i.e. design an a priori truncation with respect to $J$, thus leading to the need to run the IT pre-processing on much smaller tensors to begin with. This could help pushing calculations to yet higher-mode tensors, i.e. higher orders in the many-body expansion, and/or larger $e_{\text{max}}$. 2. 2.

The benefits of IT are expected to be even more pronounced when going to doubly open-shell nuclei that necessitate to work in an $m$-scheme basis. For $e_{\text{max}}=12$, a full-fledged BCCSD calculation would require over $3$ Tb of storage Si15 such that only techniques like IT could make such a calculation manageable. 3. 3.

Similar to many-body tensors in CC theory, IT can be straightforwardly applied to Bogoliubov configuration interaction truncated up to, e.g., 8 quasi-particle excitations yielding the IT-BCISDTQ approach. Doing so requires the derivation of the $\kappa^{(2)}$ estimate of the leading quadruple contributions to the BMBPT wave function. Applying Lanczos diagonalization, ground-state and low-lying excitation energies could be accessed, as well as other observables. Eventually, particle-number symmetry can be consistently restored by applying projection techniques. The lack of size-extensivity of truncated (B)CI calculations can be approximately cured via an a posteriori corrections of Davidson type.

5 Computational analysis of TF and IT

Typically, expansion methods admit polynomial scaling with respect to system size212121In practice, the cost is also driven by a prefactor that can be sufficiently different in two methods to obscure the classification associated to the scaling law over a certain interval of $N$ values.

[TABLE]

where $\alpha$ is a characteristic exponent of the theory222222Valence-space approaches involve the diagonalization of a dressed Hamiltonian in an active space of limited size. Consequently, such methods are excluded from the present analysis even though the dressing of the Hamiltonian itself requires only polynomial effort.. Table 4 displays the value of $\alpha$ for a collection of many-body methods routinely used in nuclear theory. Furthermore, the order of MBPT-completeness associated to each method/truncation is added, i.e., the order up to which all MBPT diagrams are included. While TF and IT aim at reducing the computational effort required to solve the many-body problem, they operate differently to achieve this goal.

First, TF potentially allows one to decrease the exponent to an effective value $\alpha^{\prime}$ by decomposing many-body tensors and finally only invoking contractions between low-rank operators. Assuming a mild scaling of $r_{\text{THC}}$ with $N$, the complexity of the tensor network might be drastically reduced, i.e. ideally, CC theory at arbitrary truncation level could be evaluated at $N^{4}$ if the THC ranks admit (nearly) linear scaling. Although it is unclear at this point if this formal counting is contaminated in practice by a large prefactor, it is indeed very promising.

Second, while IT does not act on the scaling exponent, it decreases the size of the one-body basis to an effective value $N^{\prime}$ by only keeping a small subset of the entries. Focusing first on the double BCC amplitude ${\cal T}_{2}$, Fig. 12 shows that for $e_{\text{max}}=12$ ($\tilde{N}\approx 182$), about $10^{5}$ tuples are retained after the IT truncation when targeting a very conservative $0.01\%$ error on the second-order correlation energy, which effectively corresponds to $\tilde{N}^{\prime}\approx 45$ in the optimal storage scheme, i.e. to an effective number of major HO shells $e_{\text{max}}^{\prime}\approx 5$. Moving to the triple BCC amplitude ${\cal T}_{3}$, Fig. 17 shows that less than $10^{8}$ tuples are retained after the IT truncation when targeting a conservative $1\%$ error on the fourth-order correlation energy. This yields a very similar effective one-body basis size $\tilde{N}^{\prime}\approx 45$, i.e. $e_{\text{max}}^{\prime}\approx 5$. In conclusion, although IT does not change the scaling of the theory, it allows one to perform high-accuracy non-perturbative calculations in a realistic model space characterized by, e.g., $e_{\text{max}}=12$ at an effective price corresponding to a much reduced model space characterized by $e_{\text{max}}^{\prime}\approx 5$. This is very promising.

6 Conclusions

In this work tensor-factorization and importance-truncation techniques are introduced as two different paradigms to pre-process the nuclear many-body problem with the goal to push ab initio calculations based on expansion methods to (i) higher accuracies, (ii) doubly open-shell nuclei and (iii) nuclei with $A>100$.

Bogoliubov many-body perturbation theory calculations of semi-magic nuclei are used in small model spaces to test the numerical implementation of tensor-factorization and importance-truncation techniques. Both methods yield very promising computational advantages with respect to storage requirements of many-body tensors while generating a small error on nuclear ground-state energies. Eventually, large-scale IT benchmarks based on the state-of-the-art computation of perturbative triple amplitudes in Bogoliubov coupled cluster theory confirm the great promises of this pre-processing method.

The next step consists of implementing THC and IT techniques in actual non-perturbative calculations, which seems particularly straightforward for IT. While these pre-processing method can be adapted to any method, our focus will be on Bogoliubov extensions of CC and truncated CI with the goal to address open-shell physics from simple single-reference methods.

Acknowledgements

We thank Robert Roth for providing us with nuclear NN and 3N matrix elements. This publication is based on work supported in part by the framework of the Espace de Structure et de réactions Nucléaires Théorique (ESNT) at CEA.

Appendix A Angular momentum coupling

In all applications performed in this work a spherical formulation of the underlying many-body theory is employed and, consequently, the evaluation of working expressions must be performed in an angular-momentum coupled scheme. The aim of this section is to discuss the implementational details of the evaluation of importance measures in such a spherical scheme requiring an angular-momentum reduction of the underlying many-body diagrams. While the basic notations and coupling symbols are presently introduced, the reader is referred to Ref. VaMo88 for an extensive treatment of angular-momentum theory.

A.1 Basics

A generic state of the one-body Hilbert space ${\cal H}_{1}$ is denoted as $|k\rangle$, where $k$ stands for a collective index

[TABLE]

with $n_{k}$ the radial quantum number, $l_{k}$ the orbital angular-momentum quantum number, $j_{k}$ the total angular-momentum quantum number with projection $m_{k}$ and $t_{k}$ the isospin projection distinguishing protons and neutrons. A reduced set of quantum numbers is introduced through

[TABLE]

where the angular-momentum projection is explicitly excluded. Effectively, Eq. (A.2) contains all quantum numbers necessary to describe a rotationally invariant system. Forming the tensor product of two one-body states, basis states of the two-body Hilbert space ${\cal H}_{2}$ are obtained in the uncoupled representation as

[TABLE]

In the coupled representation, the two total angular momenta of the states $|k_{1}\rangle$ and $|k_{2}\rangle$ are coupled to a total two-body angular momentum $J$ and projection232323While the two-body state does indeed depend on $M$, the label is omitted for brevity given that the reduced tensors eventually built in that basis are diagonal in $M$ and independent of it. $M$,

[TABLE]

where $C_{j_{k_{1}}m_{k_{1}}j_{k_{2}}m_{k_{2}}}^{JM}$ denotes the Clebsch-Gordan coefficient (CGC) mitigating the transformation from the uncoupled to the coupled basis. The inverse transformation is given by

[TABLE]

Applying the same rationale to the bra two-body states leads to defining the $J$-coupled representation of, e.g., the interaction matrix elements ${}^{J}v_{\tilde{k}_{1}\tilde{k}_{2}\tilde{k}_{3}\tilde{k}_{4}}$ via

[TABLE]

Due to rotational invariance, nuclear matrix elements are diagonal with respect to the two-body angular momenta of the bra and ket states and independent of their projection $M$.

Furthermore, angular-momentum coupling is extended to basis states of ${\cal H}_{3}$, i.e., the tensor product of three single-particle states

[TABLE]

In order to do so a coupling order needs to be fixed. In the subsequent derivations this is chosen to be

[TABLE]

where the intermediate two-body angular-momentum quantum number $J_{12}$ is integer and the three-body angular-momentum quantum number $J$ is half-integer.

While the initial choice of the coupling order for three-body states is arbitrary, a consistent treatment throughout the derivation is crucial. A proper selection of this coupling order may significantly simplify or complicate the resulting final expressions in the many-body framework.

A.2 Quasi-particle matrix elements

Since the building blocks, e.g. grand potential matrix elements, of a particle-number broken many-body formalism are defined in quasi-particle space, the angular-momentum coupling needs to be extended to these general objects. However, the different normal-ordered components admit different symmetry properties with respect to parity and $M$ conservation. In order to unify the treatment the notion of cross-coupled matrix elements is introduced yielding matrix elements with the same block structure for all $\Omega^{ij}$ components Tichai19unp . In the following all working equations are expressed in terms of such cross-coupled matrix elements indicated by a tilde on top of the corresponding symbols, e.g., $\tilde{\Omega}^{ij}$ instead of $\Omega^{ij}$.

The angular-momentum coupling of cross-coupled matrix elements of the normal-ordered grand-potential components reads as

[TABLE]

where the actual expression of angular-momentum coupled quantities ${{}^{J}}\tilde{\Omega}^{ij}$ will be detailed in a forthcoming publication. Note the appearance of an additional ($m$-dependent) phase factor for $\tilde{\Omega}^{40}$ and $\tilde{\Omega}^{31}$ compared to $\tilde{\Omega}^{22}$.

When performing the angular-momentum reduction of complex tensor networks the use of recoupling symbols is inevitable. Of particular importance is the Wigner 6j-symbol

[TABLE]

which arises naturally from the coupling of three angular momenta VaMo88 . In principle higher recoupling symbols like $9j$- and $12j$-symbols may also appear. However, due to limited cache size, it is convenient to precompute and store only $6j$-symbols. Thanks to angular-momentum identities, higher-order recoupling symbols can be reexpressed as products of $6j$-symbols that indeed constitute indeed the most complex recoupling symbols needed in our calculations.

The reduction of both $\tilde{t}^{40(1)}$ and $\tilde{t}^{60(2)}$ is now performed. The construction of $\tilde{t}^{20(1)}$ or $\tilde{t}^{20(2)}$ can be done analogously. However, the low number of entries of $t^{20}$ typically does not require an IT treatment.

A.3 Reduction of $\tilde{t}^{40(1)}$

The cross-coupled matrix elements of ${\cal T}_{2}^{(1)}$ read in $m$-scheme as

[TABLE]

and is thus trivially proportional to $\tilde{\Omega}^{40}$ such that no explicit permutation of external indices needs to be applied to ensure the full antisymmetry of $\tilde{t}^{40(1)}$. Inserting $J$-coupled matrix elements of $\tilde{\Omega}^{40}$ in Eq. A.11, one straighforwardly obtains those of interest. While the $m$-scheme expression requires $N^{4}$ evaluations, its $J$-coupled partner (A.11) only necessitates $J_{2\text{max}}\cdot\tilde{N}^{4}$ evaluations where $J_{2\text{max}}$ defines the number of channels of the two-body angular-momentum $J$. The number of $J$ channels is $J_{2\text{max}}=2e_{\text{max}}+2$ as long as no additional truncation on $e_{2\text{max}}\equiv e_{1}+e_{2}$ is employed.

A.4 Reduction of $\tilde{t}^{60(2)}$

The $m$-scheme cross-coupled matrix elements $\tilde{t}^{60(2)}_{k_{1}k_{2}k_{3}k_{4}k_{5}k_{6}}$ associated with the diagram displayed in Fig. 7 and whose algebraic form was given in Eq. (40) can be written in a compact form as

[TABLE]

where the operator $\mathcal{P}(k_{1}k_{2}k_{3}/k_{4}k_{5}k_{6})$ permutes indices from the first group with indices from the second group in all possible ways. The permutation operator generates $20$ different terms yielding a complex expression for $\tilde{t}^{60(2)}$.

The angular-momentum coupling of the cross-coupled matrix elements of the triple BCC amplitude reads as

[TABLE]

where in both bra and ket states the first two quasi-particle labels are coupled to intermediate two-body angular momenta $J_{12}$ and $J_{45}$, respectively. Whereas the angular-momentum coupling of the double BCC amplitude matrix elements is comparatively simple, performing it for the 20 terms defining the mode-6 tensor of present interest is much more challenging and error prone due to the increasing number of CGCs. Additionally, to the four CGCs arising from the external coupling there are four additional ones from inserting the transformations from (A.9). Overall this requires the angular-momentum reduction of a string of eight CGCs. In order to deal with this problem systematically, a graph-theory based tool for automatized angular-momentum algebra of many-body tensor networks has been used ripoche19b . Much larger strings of CGCs can actually be treated in this framework, thus paving the way to relax state-of-the-art many-body truncations and extend high-accuracy ab initio calculations of open-shell nuclei. For benchmarking purposes, an independent $m$-scheme code was implemented to verify the correctness of the results from the $J$-scheme implementation in small model spaces.

Omitting the harmless energy denominator that only depends on external indices, the $J$-coupled form of the 20 contributions is given by {strip}

[TABLE]

To the second-order triple amplitude is associated the fourth-order perturbative ground-state energy correction

[TABLE]

which is a $N^{6}$ process. The $J$-coupled form of (A.15) is given by

[TABLE]

and requires $J^{2}_{2\text{max}}J_{3\text{max}}\cdot\tilde{N}^{6}$ evaluations, where $J_{k\text{max}}$ denotes the number of $k$-body angular momentum channels. For the case $k=3$ this number of channels is given by $J_{3\text{max}}=3e_{\text{max}}+2$ as long as no additional truncation on $e_{3\text{max}}\equiv e_{1}+e_{2}+e_{3}$ is employed.

Due to working in quasi-particle basis, the evaluation of Eq. (A.14) and Eq. (A.16) is a computationally challenging task requiring orders of magnitudes more resources than its symmetry-conserving counterpart in closed-shell nuclei. Taking mid-mass Ca or Ni isotopes as an example, the runtime increases by a factor $10^{2}-10^{3}$ assuming an $e_{\text{max}}=12$ model space corresponding to $\tilde{N}=182$.

Finally, it is worth noting that a consistent inclusion of all permutations is mandatory. From (A.14) it is obvious that the evaluation of some of the terms is computationally simpler than others due to the appearance of additional recoupling symbols accompanied by angular-momentum summations. However, from an $m$-scheme perspective all terms are equally important and the difference in computational complexity is just an artefact of the chosen external coupling order. In actual applications, different contributions corresponding to recoupling patterns of different complexity (i.e. different number of $6j$-symbols) are equally important such that a consistent inclusion of all terms is necessary.