In-medium k-body reduction of n-body operators
Mikael Frosini, Thomas Duguet, Benjamin Bally, Yann Beaujeault-Taudière, Jean-Paul Ebran, Vittorio Som?
NNoname manuscript No. (will be inserted by the editor)
In-medium k -body reduction of n -body operators A flexible symmetry-conserving approach based on the sole one-body density matrix
M. Frosini , T. Duguet , B. Bally , Y. Beaujeault-Taudi`ere , J.-P.Ebran , V. Som`a IRFU, CEA, Universit´e Paris-Saclay, 91191 Gif-sur-Yvette, France KU Leuven, Department of Physics and Astronomy, Instituut voor Kern- en Stralingsfysica, 3001 Leuven, Belgium Departamento de F´ısica Te´orica, Universidad Aut´onoma de Madrid, 28049 Madrid, Spain CEA, DAM, DIF, 91297 Arpajon, France Universit´e Paris-Saclay, CEA, Laboratoire Mati`ere en Conditions Extrˆemes, 91680, Bruy`eres-le-Chˆatel, FranceReceived: February 23, 2021 / Revised version: date
Abstract
The computational cost of ab initio nuclearstructure calculations is rendered particularly acute bythe presence of (at least) three-nucleon interactions.This feature becomes especially critical now that many-body methods aim at extending their reach beyondmid-mass nuclei. Consequently, state-of-the-art ab initiocalculations are typically performed while approximat-ing three-nucleon interactions in terms of effective, i.e.system-dependent, zero-, one- and two-nucleon opera-tors. While straightforward in doubly closed-shell nu-clei, existing approximation methods based on normal-ordering techniques involve either two- and three-bodydensity matrices or a symmetry-breaking one-body den-sity matrix in open-shell systems. In order to avoid suchcomplications, a simple, flexible, universal and accurateapproximation technique involving the convolution ofthe initial operator with a sole symmetry-invariant one-body matrix is presently formulated and tested numeri-cally. Employed with a low-resolution Hamiltonian, thenovel approximation method is shown to induce errorsbelow 2 −
3% across a large range of nuclei, observablesand many-body methods.
Dealing fully with three-, possibly four-, nucleon inter-actions is non-trivial but tractable in a self-consistentmean-field Hartree-Fock (HF) or Hartree-Fock Bogoli-ubov (HFB) calculation. However, it becomes extremelycumbersome, if not impossible, beyond a certain nuclearmass when solving the A -body Schr¨odinger equation togood-enough accuracy beyond the mean field. Conse-quently, ab initio calculations of mid-mass nuclei aretypically performed on the basis of the so-called normal- ordered two-body (NO2B) approximation that capturesdominant effects of three-nucleon interactions while ef-fectively working with two-nucleon operators [1,2]. Inlarge-scale no-core shell-model calculations, the errorinduced by the NO2B approximation of the Hamilto-nian was estimated to be of the order of 1-3% up to theoxygen region .The NO2B approximation was originally designed bynormal ordering the Hamiltonian with respect to a Slaterdeterminant through standard Wick’s theorem [5]. Theprocedure involved the contraction of the three-body op-erator with the one-body density matrix of that product-state Slater determinant. The approximate Hamiltonianresulting from the NO2B approximation was consis-tently employed in many-body methods applicable toclosed-shell nuclei expanding the exact solution withrespect to a symmetry-conserving , e.g. J Π = 0 + , Slaterdeterminant. In this context, the approximate NO2BHamiltonian naturally displays the same symmetries asthe original one. However, the naive extension of theNO2B approximation to methods applicable to open- For low-resolution Hamiltonians obtained via, e.g., the ap-plication of similarity renormalization group (SRG) transfor-mations [3], the efficiency of the NO2B approximation canbe understood on the basis of phase-space arguments in thecalculation of homogeneous infinite nuclear matter [4]. In par-ticular, the analysis of Ref. [4] makes clear that the quality ofthe approximation can only improve as the density (mass) ofmatter (nuclei) increases. In the present work, a symmetry-conserving state representsa state whose associated one-body density matrix is symmetry-invariant, i.e. belongs to the trivial irreducible representationof the symmetry group of the Hamiltonian. While for theSU(2) group it makes necessary for the many-body state itselfto be symmetry invariant, i.e. to be a J = 0 state, for the U(1)group this condition is automatically satisfied for the normal one-body density matrix. a r X i v : . [ nu c l - t h ] F e b shell systems via the use of symmetry-breaking referencestates poses a difficulty in that respect. Indeed, ignoringthe normal-ordered three-body term delivers in sucha situation an approximate operator that itself explic-itly breaks the corresponding symmetry(ies) of the fullHamiltonian. This feature is unwelcome as it lacks thetransparency of restricting the symmetry breaking to ap-proximations of the many-body state, especially in viewof the eventual restoration of the symmetry(ies).Within the frame of many-body methods breaking U(1)symmetry associated with particle-number conservation[6,7,8,9,10] via the use of Bogoliubov reference states,a particle-number-conserving normal-ordered k -body(PNOkB) approximation of an arbitrary n -body opera-tor was recently formulated and validated numerically[11]. Using the PNOkB approximation, ab initio calcula-tions on singly open-shell nuclei based on U(1)-breakingand restored formalisms can thus be safely performed.As one currently becomes interested in methods (fur-ther) breaking SU(2) symmetry associated with angular-momentum conservation [12,13] to describe doubly open-shell nuclei, the symmetry-conserving NOkB approxima-tion should be extended to this symmetry group, whichhappens to be neither easy nor transparent.The difficulty is bypassed from the outset when de-scribing open-shell systems through a so-called multi-reference method based on an explicitly correlated andsymmetry-conserving reference state, e.g. in the multi-reference in-medium similarity renormalization groupmethod [14,15]. In this context, it is natural to approx-imate the three-body interaction through its normal-ordering with respect to the correlated reference state onthe basis of Kutzelnigg-Mukherjee’s Wick theorem [16].The benefit however comes with the prize of havingto contract the three-body operator not only with theone-body, but also with the two-body and three-bodydensity matrices of the correlated state. A somewhatsimilar situation occurs within self-consistent Green’sfunction (SCGF) theory that can be formulated in termsof effective k -body vertices obtained by contracting ini-tial n -body operators ( n ≥ k ) with fully correlated( n − k )-body density matrices [17].In conclusion, several approaches exist to produce so-called effective , i.e. nucleus-dependent, interactions. Theaim is to eventually discard the effective operator(s)of highest n -body character(s) whose contribution to,e.g., ground-state energies is (are) expected to be muchsmaller than for the original operator(s) carrying thesame n -body character(s). Such a procedure always Once again, this property can be justified for low-scale Hamil-tonians on the basis of phase-space arguments [4]. involves a contraction of the original operator(s) witha (set of) density matrix (matrices) reflecting (i) thesymmetries and (ii) the correlations of the many-bodystate it (they) originates from and that is typically thereference state or the fully correlated state at play inthe many-body method of interest.In this context, the present work introduces a novelmethod to build a set of effective k -body interactions inview of approximating the initial Hamiltonian. While theHamiltonian is indeed our primary target, the procedurecan in principle be applied to any observable. Our goalis thus to formally justify and test numerically a novelapproximation method that1. only invokes contractions with a one-body densitymatrix,2. uses a symmetry-invariant one-body density matrix,3. is flexible regarding the many-body state used tocompute that one-body density matrix,4. re-expresses the approximate Hamiltonian in normal-ordered form with respect to the particle vacuum.The benefits are that1. the method does not involve l -body density matriceswith l > k -body terms always possessesthe same symmetry group as the original one,3. the method does not necessarily have to employ theone-body density matrix associated with the many-body (reference) state at play in the method usedto solve Schr¨odinger’s equation,4. the resulting Hamiltonian is explicitly expressed inthe original single-particle basis such that it cannaturally be employed as the starting point of anymany-body method.Points (3) and (4) underline the fact that the approx-imation of the Hamiltonian and the resolution of theSchr¨odinger equation, although not unrelated, consti-tute two different problems and do not necessarily haveto be dealt with on the basis of the same many-bodyscheme.Per se, the method is applicable independently of theclosed- or (doubly) open-shell character of the system aswell as of the ground or excited nature of the targetedstate. Still, point (2) implies that only one-body densitiesderiving from a J Π = 0 + state can be employed inthe approximation procedure, which obviously impliesthat the employed one-body density matrix does not necessarily derive from the targeted state/nucleus. Withexcited states of even-even nuclei in mind, one can mostnaturally approximate the Hamiltonian through theuse of a one-body density matrix associated with theground state. With odd-even or odd-odd systems inmind, one can employ the symmetry-invariant densitymatrix associated with a fake odd system describedin terms of, e.g., a statistical mixture [18,19]. In thepresent paper, the focus is on even-even systems.The paper is organized as follows. Sec. 2 is dedicatedto the formulation of the method and its relation toexisting ones. After explicating in Sec. 3 the hierarchyof one-body density matrices and many-body methodspresently employed to test the approximation method,the corresponding numerical results are presented inSec. 4. While conclusions are provided in Sec. 5, severalappendices complement the paper with useful technicaldetails. An arbitrary particle-number-conserving operator O canbe written as O ≡ N (cid:88) n =0 O nn , (1)where each n -body component reads in an arbitrary ba-sis { c † a , c a } of the one-body Hilbert space H as O nn ≡ n ! 1 n ! (cid:88) a ··· a n b ··· b n o a ··· a n b ··· b n A a ··· a n b ··· b n , (2)where A a ··· a n b ··· b n ≡ c † a · · · c † a n c b n · · · c b (3)denotes a string of n one-particle creation and n one-particle annihilation operators such that ( A a ··· a n b ··· b n ) † = A b ··· b n a ··· a n . This string is in normal order with respect tothe particle vacuum | (cid:105) , i.e. N ( A a ··· a n b ··· b n ) = A a ··· a n b ··· b n , (4)where N ( . . . ) denotes the normal ordering with respectto | (cid:105) .In Eq. (2), the n -body matrix elements { o a ··· a n b ··· b n } consti-tute a mode-2 n tensor denoted as o ( n ) , i.e. a data arraycarrying 2 n indices associated with the n ( n ) particle creation (annihilation) operators they multiply. The n -body matrix elements are fully anti-symmetric underthe exchange of any pair of upper or lower indices, i.e. o a ··· a n b ··· b n = (cid:15) ( σ u ) (cid:15) ( σ l ) o σ u ( a ··· a n ) σ l ( b ··· b n ) , (5)where (cid:15) ( σ u ) ( (cid:15) ( σ l )) refers to the signature of the per-mutation σ u ( . . . ) ( σ l ( . . . )) of the n upper (lower) in-dices. The l -body density matrix associated with a many-bodystate | Θ (cid:105) constitutes a mode-2 l tensor defined through (cid:104) ρ ( l ) Θ (cid:105) b ··· b l a ··· a l ≡ (cid:104) Θ | A a ··· a l b ··· b l | Θ (cid:105)(cid:104) Θ | Θ (cid:105) . (6)In the following, the superscripts l and Θ are omitted for l = 1 and when dealing with a generic density matrix,respectively. The elements of ρ ( l ) Θ inherit from A a ··· a l b ··· b l a full anti-symmetry under the exchange of any pair ofupper or lower indices along with a hermitian character,i.e. (cid:104) ρ ( l ) Θ (cid:105) b ··· b l a ··· a l = (cid:18)(cid:104) ρ ( l ) Θ (cid:105) a ··· a l b ··· b l (cid:19) ∗ . (7)Given two density matrices ρ ( l ) Θ and ρ ( k ) Ψ , their tensorproduct ρ ( l ) Θ ⊗ ( k ) Ψ ≡ ρ ( l ) Θ ⊗ ρ ( k ) Ψ (8)defines a direct-product ( l + k )-body density matrixthrough the mode-2( l + k ) tensor whose elements aregiven by (cid:104) ρ ( l ) Θ ⊗ ( k ) Ψ (cid:105) b ··· b l + k a ··· a l + k ≡ (cid:104) ρ ( l ) Θ (cid:105) b ··· b l a ··· a l (cid:104) ρ ( k ) Ψ (cid:105) b l +1 ··· b l + k a l +1 ··· a l + k , (9)and display the hermitian property characterized inEq. (7). Because of the direct-product character of ρ ( l ) Θ ⊗ ( k ) Ψ , its elements are only partially anti-symmetrized,i.e. they are anti-symmetric under the exchange of anypair of the first l (or last k ) upper or lower indices.In case one considers the m -fold tensor product of thesame l -body density matrix ρ ( l ) Θ , the notation can befurther simplified according to ρ ⊗ ( ml ) Θ ≡ ρ ( l ) Θ ⊗ . . . ⊗ ρ ( l ) Θ . In particular, the m -fold tensor product of the Conventionally, Eq. (6) is consistently extended to l = 0 via ρ (0) Θ ≡ generic one-body density matrix ρ defines a mode-2 m tensor whose elements are (cid:104) ρ ⊗ ( m ) (cid:105) b ··· b m a ··· a m ≡ ρ b a · · · ρ b m a m . (10)Because of its pure direct-product character, the ele-ments of ρ ⊗ ( m ) display no property under the exchangeof any pair of upper or lower indices but inherit thehermitian property characterized by Eq. (7). In the following, the extent to which two one-bodydensity matrices ρ and ρ (cid:48) deviate from one another willneed to be characterized. The distance d ( ρ, ρ (cid:48) ) ≡ (cid:107) ρ − ρ (cid:48) (cid:107) , (11)provides such a diagnostic, with (cid:107) . (cid:107) the Frobenius norm reading for an arbitrary mode- n tensor T as (cid:107) T (cid:107) ≡ (cid:115) (cid:88) i ...i n T i ...i n T ∗ i ...i n , (12)where the superscript denotes elementwise complex con-jugation. The convolution of the mode-2 n tensor o ( n ) associatedwith a n -body operator O nn with the mode-2 m ten-sor ( m ≤ n ) defining a m -body density matrix ρ ( m ) generates the mode-2( n − m ) tensor o ( n ) · ρ ( m ) with ele-ments (cid:104) o ( n ) · ρ ( m ) (cid:105) a ··· a n − m b ··· b n − m ≡ (cid:88) a n − m +1 , ··· ,a n b n − m +1 , ··· ,b n o a ··· a n b ··· b n (cid:104) ρ ( m ) (cid:105) b n − m +1 ...b n a n − m +1 ...a n . (13)The tensor o ( n ) · ρ ( m ) is obviously a pure number whenever m = n and nothing but the initial tensor o ( n ) whenever m = 0.Given two density matrices ρ ( l ) Θ and ρ ( k ) Ψ , it is straight-forward to check that the convolution is such that thefollowing identity holds (cid:16) o ( n ) · ρ ( m ) Θ (cid:17) · ρ ( l ) Ψ = (cid:16) o ( n ) · ρ ( l ) Ψ (cid:17) · ρ ( m ) Θ = o ( n ) · (cid:16) ρ ( m ) Θ ⊗ ρ ( l ) Ψ (cid:17) . (14) 2.2 Standard NOkB approximation Let us consider a symmetry-conserving product state | Φ (cid:105) , i.e. a J Π = 0 + Slater determinant. Standard Wick’stheorem [5] with respect to | Φ (cid:105) entails four elementarycontractions c † a c † b − : c † a c † b : = 0 , (15a) c † a c b − : c † a c b : = ρ Φba , (15b) c a c † b − : c a c † b : = δ ab − ρ Φab , (15c) c a c b − : c a c b : = 0 , (15d)where : . . . : denotes the normal ordering with respectto | Φ (cid:105) .Applying Wick’s theorem, the operator O defined inEq. (1) is rewritten as O = N (cid:88) k =0 O kk [ ρ Φ ] , (16)where O kk [ ρ Φ ] is a k -body operator in normal-orderedform with respect to | Φ (cid:105) O kk [ ρ Φ ] ≡ k ! 1 k ! (cid:88) a ··· a k b ··· b k o a ··· a k b ··· b k [ ρ Φ ] : A a ··· a k b ··· b k : . (17)Considering O nn ( n ≤ N ) and k ≤ n , there are( n − k )! (cid:18) nn − k (cid:19)(cid:18) nn − k (cid:19) (18)ways to perform ( n − k ) non-zero contractions. Conse-quently, the matrix elements of O kk [ ρ Φ ] are related tothose defining the original contributions to O through o a ··· a k b ··· b k [ ρ Φ ] = N (cid:88) n = k n − k )! (cid:104) o ( n ) · ρ Φ ⊗ ( n − k ) (cid:105) a ··· a k b ··· b k . (19) The normal-ordered k -body (NOkB) approximation O NOkB [ ρ Φ ] to the operator O proceeds by truncatingthe sum in Eq. (16) to the desired maximum value k .While the original operator is obviously independent of ρ Φ , O NOkB [ ρ Φ ] does acquire such a dependence as soonas k < N . In the present paper, a k -body operator and the tensorrepresenting it are written with a standard font (bold font),e.g. O kk ( O kk ), if the operator is in normal order with respectto the particle vacuum | (cid:105) (a many-body state | Φ (cid:105) ). For example, the standard NO2B approximation consistsof ignoring beyond normal-ordered 2-body terms todefine the approximate Hamiltonian as [1,2] H NO2B [ ρ Φ ] ≡ H [ ρ Φ ] + H [ ρ Φ ] + H [ ρ Φ ] . (20)Generalizing the approach to a U(1)-breaking productstate, i.e. a Bogoliubov reference state, standard Wick’stheorem gives rise to non-zero anomalous contractions(Eqs. (15a) and (15d)) such that the truncation pro-cedure generates a particle-number-breaking operator.A different truncation scheme was thus formulated todesign a particle-number conserving normal-ordered k -body (PNOkB) approximation in Ref. [11]. A similarproblem arises when using a SU(2) non-invariant Slaterdeterminant, i.e. whenever | Φ (cid:105) is not a J Π = 0 + state.Indeed, the standard NOkB approximation delivers anoperator that is not rotationally invariant in such a case.Rather than extending the tedious approach designed inRef. [11] for the U(1) case, a novel method is proposedin Sec. 2.3 that avoids such complications from the out-set by involving the one-body density matrix stemmingfrom a symmetry-conserving many-body state. Starting from Eq. (16), it is interesting to re-express theoperator back into a normal-ordered form with respectto the particle vacuum (Eq. (2)). Doing so requires toapply Wick’s theorem backward, i.e. with respect to | (cid:105) .To do so, the only required non-zero contraction is givenby: c † a c b : − N (: c † a c b :) = : c † a c b : − c † a c b = − ρ Φba , (21)which is nothing but the opposite of the elementarycontraction at play in the first step. The original n -body part of O is obtained back in terms of the variouscontributions entering Eq. (16) such that the connectionbetween their matrix elements is given by o a ··· a n b ··· b n = N (cid:88) l = n ( − l − n ( l − n )! (cid:104) o ( l ) [ ρ Φ ] · ρ Φ ⊗ ( l − n ) (cid:105) a ··· a n b ··· b n . (22)Truncating Eq. (16) according to the NOkB approxima-tion and inserting the result into Eq. (22) delivers thematrix elements of the approximate n -body part of O innormal order with respect to the particle vacuum.2.3 Generalized k -body approximationThe standard NOkB approximation relies on standardWick’s theorem and is thus strictly defined with respect to a symmetry-conserving many-body product state. Be-cause this restriction is too severe in open-shell systems,a generalization of the procedure is now envisioned suchthat the involved one-body density matrix can originatefrom a more general many-body state. Given the operator O and the one-body density matrix ρ associated with an arbitrary J Π = 0 + state, one first defines the set of anti-symmetrized matrix elements o a ··· a k b ··· b k [ ρ ] ≡ N (cid:88) n = k n − k )! (cid:104) o ( n ) · ρ ⊗ ( n − k ) (cid:105) a ··· a k b ··· b k , (23)in strict analogy with Eq. (19) but relaxing the neces-sity for the density matrix to originate from a Slaterdeterminant .The key point of the present development relates tothe fact that, independently of the nature of ρ , theinverse operation embodied by Eq. (22) remains valid inthe present context and recovers the original operator’smatrix elements, i.e. o a ··· a n b ··· b n = N (cid:88) l = n ( − l − n ( l − n )! (cid:104) o ( l ) [ ρ ] · ρ ⊗ ( l − n ) (cid:105) a ··· a n b ··· b n . (24)This identity is proven in App. A. One can thus concludethat the combined operations embodied by Eqs. (19)and (22) are actually valid outside the reach of standardWick’s theorem. Indeed, the two steps are utterly generaloperations, i.e. tensor products that are inverse fromone another, holding independently of the nature of ρ (i.e. whether it stems from a Slater determinant or not).Standard Wick’s theorem is recovered as a particularcase of the general tensor identities (23)-(24), namelywhen the one-body density matrix does originate froma Slater determinant.Thus, the purpose of Eqs. (23) and (24) is to start fromthe set of tensors defining each n -body contribution tothe original operator O in Eqs. (1)-(2) and to recover itafter having gone through an intermediate set defined in The matrix elements introduced in Eq. (23) are not ob-tained through a set of algebraic operations on the originaloperator but via a straight convolution of tensors. Still, thepresent procedure could be formulated within the frame ofthe quasi -normal ordering of Ref. [20], which is itself an exten-sion of Kutzelnigg and Mukherjee’s universal normal-orderinginvolving the sole one-body density matrix. In this context, itbecomes possible to associate an actual quasi-normal-orderedoperator to the tensor o ( k ) [ ρ ]. However, given that such aquasi-normal-ordered operator is of no use in the presentcontext, there is no need to invoke it. strict analogy with the tensors generated via the single-reference normal ordering. While there is no benefit inapplying the two-step procedure per se , it ensures thatthe original operator is exactly recovered when doing so.Based on this property, the method provides a usefulway to produce nucleus-dependent approximations tothe operator through the truncation of the intermediateset of tensors. In close analogy with the NOkB approximation, the k -body approximation of O is now introduced. First,the set of tensors defined through Eq. (23) is truncatedaccording to ¯o ( l ) [ ρ ] ≡ o ( l ) [ ρ ] for l ≤ k , (25a) ¯o ( l ) [ ρ ] ≡ l > k . (25b)Second, inserting Eq. (25) into Eq. (24) generates theset of tensors ¯ o ( n ) [ ρ ] defining the k -body approximationof O in normal order with respect to the particle vacuumaccording to O kB [ ρ ] ≡ k (cid:88) n =0 ¯ O nn [ ρ ] , (26)where the truncation of the sum naturally derives fromEq. (25). While the original operator O is independentof ρ , O kB [ ρ ] does acquire such a dependence as a resultof the truncation characterized by Eq. (25).While O kB can be built on the basis of an arbitrarily cor-related (symmetry-conserving) state, it does not requirethe use of ρ ( l ) with l >
1. As a result, the procedure issignificantly simpler than the one associated with theapplication of Kutzelnigg-Mukherjee’s Wick theorem orthe one at play in SCGF theory. The practicality of theapproach also relates to the fact that the effective Hamil-tonian is expressed in normal-ordered form with respectto the particle vacuum in the working single-particlebasis. As a result, H kB [ ρ ] can be straightforwardly usedin place of H as the input to any many-body methodof interest. One is typically interested in the 2-body approximation H B [ ρ ] of an initial Hamiltonian containing a 3-body interaction H ≡ T + V + W ≡ (cid:88) a b t a b A a b + 1(2!) (cid:88) a a b b v a a b b A a a b b + 1(3!) (cid:88) a a a b b b w a a a b b b A a a a b b b , (27)where t a b denotes matrix elements of the kinetic en-ergy whereas v a a b b and w a a a b b b denote anti-symmetricmatrix elements of two- and three-body interactions,respectively.Setting O −→ ,O −→ T ,O −→ V ,O −→ W ,
Eq. (25) gives for k = 2 ¯h (0) [ ρ ] ≡ t (1) · ρ + 12! v (2) · ρ ⊗ (2) + 13! w (3) · ρ ⊗ (3) , (28a) ¯h (1) [ ρ ] ≡ t (1) + v (2) · ρ + 12! w (3) · ρ ⊗ (2) , (28b) ¯h (2) [ ρ ] ≡ v (2) + w (3) · ρ , (28c) ¯h (3) [ ρ ] ≡ . (28d)Except for the key fact that ρ does not necessarily relateto a Slater determinant, Eq. (28) is formally identicalto Eq. (20) defining H NO2B . Inserting Eq. (28) intoEq. (24), one eventually obtains the three tensors¯ h (0) [ ρ ] ≡ w (3) · ρ ⊗ (3) , (29a)¯ h (1) [ ρ ] ≡ t (1) − w (3) · ρ ⊗ (2) , (29b)¯ h (2) [ ρ ] ≡ v (2) + w (3) · ρ , (29c)defining the normal-ordered contributions to H B [ ρ ]with respect to the particle vacuum, i.e. H B [ ρ ] =¯ h (0) [ ρ ]+ 1(1!) (cid:88) a b ¯ h a b [ ρ ] A a b + 1(2!) (cid:88) a a b b ¯ h a a b b [ ρ ] A a a b b . (30) In addition to the fact that, by construction, H B [ ρ ] doesnot contain a three-body operator, its structure differsfrom the original operator expressed in normal orderwith respect to the particle vacuum (Eq. (27)) by thefact that it incorporates the pure number ¯ h (0) [ ρ ]. Equations (29) define a set of nucleus-dependent 0-, 1-and 2-body operators entering H B [ ρ ]. As in the NO2Bapproximation, the inclusion of a large part of W intothese effective operators, while treating T and V exactly,gives a clear argument that omitting h (3) [ ρ ] leads tosmall, hopefully small enough, errors. Still, one is leftwith the question of the optimal character of the one-body density matrix to be employed for a given systemand many-body approximation.In the hypothesis that exact eigenstates of H in the A -body Hilbert space H A are known, one may expectthat employing the one-body density matrix of the exactground-state is optimal to reproduce the ground-stateenergy . In fact, this intuition is not correct.From a formal viewpoint, it would be interesting tofind the optimal one-body density matrix to be used in H B [ ρ ] to best reproduce, e.g., the energy associatedwith the (approximate) ground state | Ψ (cid:105) of the fullHamiltonian H . This would however not be of practicaluse. Consequently, the numerical results displayed inSec. 4 rely on testing a set of trial one-body densitymatrices while obtaining the solution to the A -bodyproblem via various approximation methods. As will beconcluded, the results are very robust with respect tothe employed one-body density matrix as long as thelatter respects a minimal set of properties.3.1 Many-body methodsThe many-body methods presently used to solve the A -body Schr¨odinger equation for a collection of doublyclosed, singly open-shell and doubly open-shell even-evennuclei (to be specified later on) are1. axially deformed Hartree-Fock-Bogoliubov (dHFB)theory [21,13], This reasoning has of course a chance to be correct only ifthe target state | Ψ (cid:105) is a J Π = 0 + state. If not, the optimaldensity matrix cannot be equal to ρ Ψ for symmetry reasonsas already briefly discussed in the introduction. One may further think that the density matrix associatedwith a statistical symmetry-conserving average of a set of exactlow-lying states is optimal to best reproduce the low-lyingspectroscopy.
2. the particle-number- and angular-momentum-projectedHFB (PHFB) method [21,22,13] based on dHFBstates,3. the projected generator coordinate method (PGCM)[21,13] mixing PHFB states along the axial quadrupolemoment of the underlying dHFB states,4. quasi-particle random phase approximation for axi-ally deformed and superfluid nuclei (dQRPA) in thefinite amplitude method (FAM) formulation [23,24],5. deformed Bogoliubov many-body perturbation the-ory at third order (dBMBPT(3)) [25,26,27,13].Deformed HFB theory constitutes the mean-field base-line that can capture the bulk of static correlations inopen-shell nuclei through the spontaneous breaking ofU(1) and SU(2) symmetries. Based on it, PHFB, PGCMand dQRPA on the one hand and dBMBPT on the otherhand, provide systematic beyond-mean-field extensionswhose aim is to capture many-body correlations. WhilePHFB, PGCM and dQRPA do so via the additionof static correlations associated with the restorationof broken symmetries and the fluctuation of shapes,dBMBPT targets dynamical correlations through the re-summation of elementary, i.e. quasi-particle, excitations.Former approaches are well suited to the descriptionof spectroscopy whereas the latter naturally addressesabsolute binding energies and associated ground-stateobservables.The string of dHFB, PHFB, PGCM and dQRPA calcula-tions can presently be performed with the full inclusionof three-body forces, i.e. employing a realistic nuclearHamiltonian H without any form of approximation.Such ab initio calculations are the first of their kind [13,24] and allow us to benchmark the approximation of H by H B [ ρ ] on the basis of non-trivial many-body Whereas BMBPT has already been applied quite systemati-cally to semi-magic spherical nuclei [9,10], it is the first timeit is performed on top of a deformed Boboliubov state [13] inview of describing doubly open-shell nuclei. While helpful to discuss the performance of a many-bodymethod, the distinction between dynamical and static cor-relation effects involve a fuzzy boundary, prominently dis-played in the dQRPA case. Namely, the dQRPA equationscan be derived within different frames, e.g. as a harmoniclimit of the GCM equations [28] or via the linearization oftime-dependent HFB equations [29,30]. Depending on theseviewpoints, dQRPA either falls in the category of post-HFBextensions grasping static correlations (associated with fluc-tuation of shapes), or in the category of beyond-mean-fieldapproaches aiming at capturing dynamical correlations (interms of 2-quasi-particle excitations). In the present work, wemake the arbitrary choice to categorize dQRPA among theformer class of approaches. methods . While it can be envisioned to do so in thefuture [26], BMBPT is however not implemented yetwith full three-body interactions. Deformed BMBPTcalculations are thus presently performed with H B [ ρ ]for various approximations to ρ and compared to thosedone earlier [9,10] on the basis of the PNO2B approxi-mation [11].3.2 Trial one-body density matricesEmploying the many-body schemes introduced above,the goal is to approximate H by H B [ ρ ] with ρ computedfrom various J Π = 0 + trial states , i.e.1. spherical harmonic oscillator Slater determinant ( ρ sHOSD ),2. spherical HF(B) state ( ρ sHF(B) ),3. PHFB state ( ρ PHFB ),4. PGCM ground-state ( ρ PGCM ),5. standard spherical MBPT ground-state ( ρ sMBPT ).In the numerical results discussed in Sec. 4, ρ sHF(B) , ρ PHFB and ρ PGCM are extracted from the correspond cal-culations performed with the full H . Contrarily, ρ sMBPT is obtained from a calculation perform with the PNO2Bapproximation whereas ρ sHOSD does not require any apriori calculation.The two options ρ sHF(B) and ρ sMBPT originate from symmetry-restricted HFB and BMBPT calculations, i.e. spherical
HF(B) ensures the J Π = 0 + character of the PHFB and PGCM calculations based on realistic chiralHamiltonians have been performed recently for the first timebut at the price of approximating three-body operators [31,15]. The exact treatment of W in realistic PGCM calculationstypically increases the CPU time by three orders of magnitudecompared to using H B [ ρ ] [13]. Because correlations captured by QRPA do not feedbackinto the ground-state, there is no non-trivial one-body densitymatrix ρ sQRPA associated with the spherical QRPA solutionto be used in the construction of H B [ ρ ]. In open-shell nuclei, the invariant density matrix is obtainedvia the use of the equal filling approximation. This approachcan be justified on the basis of a specific statistical mixture ofsHO Slater determinants carrying the appropriate number ofparticles [19] or on the basis of a specific linear combination ofsHO Slater determinants carrying different number of particlessuch that the linear combination has the correct number ofparticles on average [32]. Standard spherical MBPT denotes many-body perturbationtheory based on a spherical Slater determinant reference staterather than on a particle-number-breaking Bogoliubov refer-ence state. The former is automatically obtained from thelatter in closed-shell nuclei where the dHFB reference statereduces to the spherical HF Slater determinant. state whereas standard spherical
MBPT ensures particle-number conservation . In the latter case, the restrictionimplies that the use of ρ sMBPT is limited to doublyclosed-shell nuclei.While the expression of ρ sHF(B) is textbook material [21],it is not the case for ρ PHFB and ρ PGCM . Consequently,the derivation of the corresponding expressions are pro-vided in App. B. For the sake of generality and futureuse [13], the derivation is actually performed for a moregeneral quantity than presently needed, i.e. App. Bprovides the expression of the transition one-body den-sity matrix between two arbitrary initial ( J Π i i ) andfinal ( J Π f f ) PGCM states. The result of present interestis then obtained by setting J Π i i = J Π f f = 0 + . Whilethe expression for ρ sMBPT is known material [33,34],the expression of ρ BMBPT it presently derives from isnot. Consequently, the derivation of ρ BMBPT is pro-vided in App. C for the sake of completeness and futureuse. O, Ca),2. singly open-shell ( O),3. doubly open-shell ( Ne, Ne, , Mg, , Ar),systems. The goal is to cover oblate, spherical andprolate representatives among which some nuclei are soft and others are hard with respect to axial deforma-tion . While it is not a problem to compute ρ from a particle-number-breaking state carrying the correct particle numberon average as in sHFB, it happens that ensuring the correctaverage particle-number requires a non-trivial procedure inBMBPT beyond HFB [27]. For simplicity, we thus presentlylimit ourselves to nuclei for which dBMBPT automaticallyreduces to standard spherical MBPT. Some of these nuclei, e.g. Mg, display a triaxial minimumif allowed to. Still, present calculations are restricted to axialsymmetry. e max ≡ max(2 n + (cid:96) ), where n and (cid:96) denotethe principal quantum number and the orbital angularmomentum, respectively. The value of the harmonicoscillator frequency (cid:126) ω is further needed to fully charac-terize the working basis. Except if specified otherwise,all calculations are presently performed with e max = 8and (cid:126) ω = 20 MeV. While these values do not permit togenerate fully converged calculations of all the nucleilisted above, the conclusions drawn at the end of thepaper are independent of them.When representing a n -body operator, the natural trun-cation of the tensor-product basis of the n -body Hilbertspace is set by e n max ≡ ne max . One and two-body op-erators are thus represented using e = e max and e = 2 e max , respectively. However, e = 8 , , (cid:28) e max ) is used to represent the three-nucleon inter-action given that employing 3 e max is largely beyondtoday’s capacities. This truncation will play a key roleregarding the quality of the approximation associatedwith H B [ ρ ] in medium-mass and/or neutron-rich nu-clei.The chiral effective field theory Hamiltonian H presentlyemployed combines a two-nucleon interaction at next-to-next-to-next-to-leading order (N3LO) [35,36] with aN LO three-nucleon interaction [37]. It is then evolved toa lower momentum scale λ srg via SRG transformations.While by default results obtained for λ srg = 1 .
88 fm − are discussed, λ srg = 2 .
23 fm − will also be used forcomparison.4.3 Ground-state binding energy Let us first discuss the use of H B [ ρ ] at the mean-field,i.e. dHFB, level. Fig. 1 displays the error (in %) of thecorresponding dHFB ground-state energies compared tothe reference values obtained from the full H . Resultsare provided for ρ = ρ sHOSD ( ), ρ sHF(B) ( ), ρ PHFB ( )and ρ PGCM ( ), as well as for ρ = ρ sMBPT ( ) wheneverapplicable, i.e. in doubly closed-shell nuclei. Fig. 1: Error (in %) of dHFB ground-state energiesobtained with H B [ ρ ] for the various test one-bodydensity matrices. The error corresponding to the useof ρ sHOSD for Mg amounts to 2.6% and lies outsidethe figure. Calculations are performed with e max = 8, e = 12 and λ srg = 1 .
88 fm − .One first observes that H B [ ρ ] perform well for thefive test one-body density matrices although a notabledegradation is visible for ρ = ρ sHOSD . As can be inferredfrom Tab. 1, the weaker performance of ρ = ρ sHOSD issystematic but especially pronounced as the mass and/orthe isospin-asymmetry of the system increases. Whilethe same trend is at play for ρ sHF(B) , ρ PHFB and ρ PGCM ,the error systematically remains below 0 .
2% for thesethree density matrices, with the exception of Mg whoseerror lies around 0 . ρ sHOSD , the average errorover the set is significantly larger (1 . .
6% in Mg (not shown inthe figure).From a general standpoint, it is not surprising that theerror due to the use of H B [ ρ ] is small at the mean-field level. To best appreciate this feature, let us focuson doubly closed-shell O and Ca for which dHFBreduces to sHF. As shown in Tab. 1, the error is below0 .
05% in these two nuclei for all test one-body densitymatrices but ρ sHOSD . As explained in App. E, the errorwould even be strictly zero in such a situation if ρ wereequated to the variational sHF one-body density matrixthroughout the sHF iterations based on H B [ ρ ]. Thisprocedure would be equivalent to working within theNO2B approximation, which is indeed exact at the sHFlevel, i.e. the NO2B approximation of the Hamiltonianonly impacts post-sHF methods by construction. Thefact that one rather takes ρ to be a fixed, e.g. ρ sHF ≤
30 Mass >
30 Neutron-rich All ρ sHOSD ρ sHF(B) ρ PHFB ρ PGCM ρ sMBPT Table 1: Average difference (in %) between ground-state dHFB energies computed with H B [ ρ ] and H for differentsub-categories in the test panel and the various test one-body density matrices. The neutron-rich subcategoryencompasses Ne, Mg and Ar. See Eq. (81a) for the definition of the cost function. Calculations are performedwith e max = 8, e = 12 and λ srg = 1 .
88 fm − .obtained from the full H , a priori determined one-bodydensity matrix to build H B [ ρ ] induces a marginal errorin sHF calculations.While the error remains below 0 .
05% in O and Cafor appropriate density matrices, the distinctly worseresult obtained in Ca for ρ = ρ sHOSD underlines thefact that obtaining a very accurate description is notautomatic even in this optimal situation, i.e. it is crucialthat the test density matrix contains relevant physicalinformation. Having said that, the results obtained withthe other four test density matrices are so similar thatno clear characteristic can be easily identified as far asthe optimal choice is concerned. Neither the consistencywith the employed many-body method nor the degreeof correlations encoded in the one-body density matrixseem to constitute a decisive feature. For example, ρ sHFB performs as well as the more advanced ρ sMBPT thatincorporates dynamical correlations beyond the meanfield, as can be seen in Tab. 1. We will come backrepeatedly to this question throughout the followingsections.While spherical doubly closed-shell nuclei are partic-ularly amenable to a very accurate description, it ispertinent to investigate the dependence of the approxi-mation on the axial quadrupole deformation of the HFBstate. All nuclei in the set but , O and Ca are dou-bly open-shell systems and thus spontaneously breakrotational symmetry at the dHFB level, Ne and Mgdisplaying the largest deformation of all.The upper panel of Fig. 2 displays the dHFB totalenergy curve (TEC) calculated in Ne from the full Because sHFB reduces to sHF in doubly closed-shell nuclei,notice that ρ sHFB = ρ PHFB in this case such that both testdensity matrices give identical results by construction.
Fig. 2: Upper panel: dHFB total energy curve of Ne asa function of the axial quadrupole deformation computedwith the full H . Lower panel: Error (in %) in the totalenergy curve when using H B [ ρ ] with the various testone-body density matrices. Calculations are performedwith e max = 8, e = 10 and λ srg = 1 .
88 fm − . H as a function of the axial quadrupole deformation .This nucleus is significantly deformed, the minimum The dimentionless axial quadrupole moment used in thefigures is defined as β ≡ π · . A (cid:104) Φ | r Y | Φ (cid:105)(cid:104) Φ | Φ (cid:105) . of the TEC being located at β = 0 .
45. As visiblefrom the lower panel, the error induced by H B [ ρ ] isessentially zero at sphericity , except for ρ = ρ sHOSD where it is equal to 0 . β = 0 . β = − .
4) on theprolate (oblate) side for ρ = ρ sHFB , ρ PHFB or ρ PGCM .For ρ = ρ sHOSD , the error is about twice as large alongthe TEC.There exists a trend along the TEC, the results obtainedwith ρ sHFB degrading slightly faster with the deforma-tion than those obtained with ρ PGCM and ρ PHFB . Thetrend is however not quantitatively significant as can beinferred from the systematic error over open-shell nucleiprovided in Tab. 1. Eventually, it is remarkable thatall three one-body density matrices give excellent andessentially equivalent results up to large deformations,especially given the fact that ρ sHFB does not encodeany information about deformation properties of Ne.This is a first indication of the robustness of the in-medium 2-body reduction method of 3-body interactionoperators.For orientation, it is interesting to analyze how thevarious components of H B [ ρ ] (Eqs. (29)-(30)) and H (Eq. (27)) contribute to the dHFB energy. Figure 3 de-composes the dHFB total energy accordingly in O and Ne. Results are provided for a schematic model spaceand ρ = ρ sHFB . Focusing first on O and making thehypothesis that the sHF density matrix is the same inboth calculations , the inspection of Eq. (29) makesclear that (a) the 0-body part of H B [ ρ ] is strictly equalto the sHF contribution originating from the three-bodyinteraction in H and that (b) the energy contributionassociated with the 1- and 2-body parts of H B [ ρ ] orig-inating from the three-body interaction exactly cancelout. These features are indeed observed in the upperpanel of Fig. 3 such that the total sHF energies areidentical in both calculations. While this formal analysisdoes not hold for dHFB in general, the results displayedin the lower panel demonstrate that it remains valid inpractice in a well-deformed nucleus such as Ne, whicheventually elucidates the high-quality results obtainedabove over a large set of nuclei. Contrarily to O and Ca, sHFB does not reduce to sHFat sphericity in Ne because neutrons and protons remainsuperfluid. Consequently, the error due to the use of H B [ ρ ]cannot be made strictly equal to zero by any optimization ofthe test one-body density matrix. Note that the edge of the displayed TEC lies 10 MeV(27 MeV) above the minimum on the prolate (oblate) side. This hypothesis is very well validated in practice, even moreso in a small model space such as the one employed in thepresent calculation.
Fig. 3: Contributions of the various components of H B [ ρ ] and H to the dHFB energy. Upper panel: O.Lower panel: Ne. Calculations are performed with e max = 6 and e = 6. While it is satisfying that the error induced by H B [ ρ ]is negligible at the mean-field, i.e. dHFB, level, it isto some extent expected and surely not sufficient toclaim victory. The performance of H B [ ρ ] must thus betested in beyond mean-field methods where the accuratecompensation observed above between the terms of H and those of H B [ ρ ] is not guaranteed to hold.While such a test must be carried out for various ab ini-tio methods, the present section focuses on ground-stateenergies obtained from dBMBPT that resums dynam-ical correlations in a perturbative fashion on top of a(possibly deformed and superfluid) HFB state. Presentcalculations are performed at the BMBPT(3) level thatis known to reproduce essentially exact results based onSRG-evolved interactions to better than 2% in oxygenisotopes and those computed from non-perturbative ex- Fig. 4: Difference of dBMBPT(3) ground-state energies(in %) obtained with H B [ ρ ] and within the PNO2Bapproximation of H for the various test one-body densitymatrices. Calculations are performed with e max = 8, e = 12 and λ srg = 1 .
88 fm − .pansion methods to better than 2% in semi-magic nucleiup to the nickel region [9,10].Although envisioned in the future, BMBPT calculationswith explicit three-nucleon forces are not available yet.Consequently, calculations with H B [ ρ ] are presentlybenchmarked against those obtained using the PNO2Bapproximation [11], which is the approximation em-ployed so far in all published BMBPT calculations ofsemi-magic nuclei [9,10]. In doubly closed-shell nuclei,the PNO2B approximation reduces to NO2B that hasitself being benchmarked against the use of full three-body interactions and shown to offer a typical 1 − O [1].Deformed BMBPT(3) binding energy differences (in%) are displayed in Fig. 4. Results produced withinboth approximations agree to better than 0 .
3% over thewhole set of considered nuclei, except for ρ = ρ sHOSD where the difference increases up to 2%. Just as forthe dHFB results discussed above, the use of a one-body density matrix encoding either static or dynamicalcorrelations beyond the mean-field does not have a sig-nificant impact on the quality of H B [ ρ ] such that theresults are essentially equivalent to those obtained with ρ = ρ sHF(B) . This can be confirmed quantitatively byinspecting the numbers reported in Tab. 2. Interestingly,the results also show that the average deviation is in-dependent of the closed- or open-shell character of thenuclei under consideration whereas it slightly increases with the mass even though the deviation remains tinyin all cases.These remarkable results indicate that the in-mediuminteraction and PNO2B approximation methods areequivalent as far as quantitative ab initio dBMBPT cal-culations of mid-mass nuclei are concerned. Given theearlier benchkmarking of the NO2B in doubly-closedshell nuclei, the presently developed in-medium approxi-mation method is well validated in fully-correlated bind-ing energy calculations.4.4 PHFB absolute energies and radiiIn the following, we wish to go beyond ground-stateenergies and test the in-medium approximation methodon spectroscopic properties. In order to do so, PHFB,PGCM and dQRPA calculations will be employed. Whilethese techniques resum static correlations associatedwith the restoration of broken symmetries and the fluc-tuation of shapes, they do not account for dynamicalcorrelations. As a result, whereas relative energies andspectroscopic quantities can be well converged and mean-ingful, absolute energies are not realistic, i.e. they arefar from converged ab initio values. Still, it is useful tofirst investigate how these absolute energies differ whencomputed from H and H B [ ρ ].In this section we thus analyse total (ground- andexcited-state) energies obtained at the PHFB level. Inaddition, corresponding ground-state matter radii arepresented. In doing so, the dependence of the results onnumerical parameters such as e , e max and λ srg isalso investigated. Upper panels of Fig. 5 display binding energies of thelowest-lying J Π = 0 + , + , + states obtained via PHFBcalculations with e = 8 and 12 (at fixed e max = 8).The energy of each state obtained from H ( ) iscompared to those generated from H B [ ρ ] with ρ = ρ sHOSD ( ), ρ sHF(B) ( ), ρ PHFB ( ) and ρ PGCM ( ),as well as with ρ = ρ sMBPT ( ) whenever applicable.Energies are shifted by the dHFB value obtained fromthe full H for the corresponding system such that allnuclei can be displayed on the same figure.Reference energies are well reproduced in light nuclei forall test one-body density matrices and both values of e , i.e. absolute deviations remain below 1 MeV until Mg. Increasing the mass and/or isospin asymmetryrenders the approximation more and more sensitive tothe value of e . Given that ab initio calculations are ≤
30 Mass >
30 Neutron-rich All ρ sHOSD ρ sHF(B) ρ PHFB ρ PGCM ρ sMBPT Table 2: Average difference (in %) of ground-state dBMBPT(3) energies obtained with H B [ ρ ] and within thePNO2B approximation of H for different sub-categories in the test panel and the various test one-body densitymatrices. The neutron-rich subcategory encompasses Ne, Mg and Ar. See Eq. (81c) for details on the costfunction. Calculations are performed with e max = 8, e = 12 and λ srg = 1 .
88 fm − .Fig. 5: Results of PHFB calculations with H and H B [ ρ ] for several test one-body density matrices ρ . Left andright panels display results obtained for e = 8 and 12, respectively, at fixed e max = 8. Upper panel: absoluteenergies of lowest J Π = 0 + , + , + states to which the dHFB energy obtained from H in each nucleus is subtracted.Lower panel: ground-state root-mean-square matter radii. Calculations are performed with λ srg = 1 .
88 fm − .known to be increasingly more sensitive to e withthe mass and isospin-asymmetry of the system [38], itis not surprising that any approximation of the three-nucleon displays the same feature. Going from e = 8,through e = 10 (not shown) and to e = 12, aclear convergence of the results is observed, although notquite yet for the heaviest and most neutron-rich nuclei of the panel. Eventually, converged results display asimilar error in medium-mass nuclei to the one obtainedin lighter systems, except for ρ = ρ sHOSD . Below, onlyresults obtained for the largest reachable value of e (typically 12 but not always) are shown. The lowest panels of Fig. 5 display ground-state root-mean-square matter radii (results are similar for excited-states radii). The conclusions are the same as for theenergies. Eventually, radii are extremely well reproducedfor all nuclei, states and test density matrices, with theexception of ρ = ρ sHOSD for which a slight underestima-tion is visible in the heaviest systems.Focusing on the right panels of Fig. 5, one does noticethat the situation regarding the performance of the testone-body density matrices is qualitatively and quanti-tatively similar to the one encountered in dHFB anddBMBPT(3) calculations. As soon as the results areconverged with respect to e , PHFB energies andradii obtained with H B [ ρ ] reproduce the reference re-sults equally well with all employed one-body densitymatrices but ρ sHOSD , i.e. it seems necessary (comparedto ρ sHOSD ) and sufficient (compared to ρ PHFB , ρ PGCM and ρ sMBPT ) to employ a test one-body density matrixencoding the information of the spherical mean-field, i.e. ρ sHF(B) .The above analysis is put in more quantitative terms viathe computation of systematic errors. Correspondingresults are shown in Tab. 3. By construction, PHFBresults are identical to sHF ones in doubly closed-shellnuclei given that the sole 0 + ground-state has beenconsidered for these nuclei and given that the projec-tions on particular number and angular momentumare superfluous for a sHF state. In the other nucleiwhere the projections typically add few MeV of corre-lations energy to the ground state, the average errorover J Π = 0 + , + , + , + PHFB states is essentially thesame as for dHFB ground-state energies, independentlyof the test one-body density matrix. e max Figure 6 probes the dependence of the results on thevalue of e max at fixed e . First, one notices that radiiare insensitive to e max and are perfectly reproduced. Sec-ond, no change is visible in the PHFB energies of Newhen going from e max = 8 to e max = 10. In the moreneutron-rich Ne isotope, there exists a slight change ofapproximate PHFB energies. While the agreement withthe reference results are still quantitatively good, theenergies degrade slightly when going from e max = 8 to e max = 10. The slight evolution away from the referenceresults relates in fact to the lack of convergence of theresults with respect to e discussed earlier. In thepresent case, e had to be set to 10 in order to beable to perform PHFB calculations with the explicit3-body interaction at e max = 10. Fig. 6: Same as Fig. 5 for Ne (left) and Ne (right)for e max = 8 and 10 at fixed e = 10.While pushing the calculations to large values of e would probably improve the agreement further, the over-all conclusion is that the high quality of H B [ ρ ] dependsonly mildly on e max as long as the reference calcula-tions themselves are converged enough. Although notshown, the same convergence behavior as a function of e max is at play in the BMBPT(3) ground-state energiesreported on in Sec. 4.3.2. λ srg Figure 7 probes the dependence of the results on thesmoothness of the Hamiltonian. The SRG Hamiltonianat λ srg = 2 .
23 fm − is less evolved than the one at λ srg = 1 .
88 fm − and produces spectra that are slightlyless compressed. Still, no significant dependence onthe hardness of the Hamiltonian is observed as far asthe quality of the results obtained with H B [ ρ ] is con-cerned. ≤
30 Mass >
30 Neutron-rich All ρ sHOSD ρ sHF(B) ρ PHFB ρ PGCM ρ sMBPT Table 3: Average error (in %) on absolute PHFB energies of low-lying J Π = 0 + , + , + and 6 + states for differentsub-categories in the test panel and the various test one-body density matrices. Calculations are performed with e = 8, e = 12 and λ srg = 1 .
88 fm − . See Eq. (81b) for details on the cost function.Fig. 7: Same as Fig. 6 but for two values of the SRGparameter λ srg . Calculations are performed with e max =8 and e = 12.Although not shown, the same convergence behavior asa function of λ srg is at play in the BMBPT(3) ground-state energies reported in Sec. 4.3.2. Fig. 8: Low-lying PHFB excitation spectra of doublyopen-shell nuclei. Reference results calculated from H are compared to those computed from H B [ ρ ] usingthe various one-body test density matrices. Calculationsare performed with e max = 8, e = 12 and λ srg =1 .
88 fm − .4.5 SpectroscopyHaving analyzed absolute PHFB energies and radii, weare now in position to investigate spectroscopic observ-ables. Low-lying PHFB excitation spectra of doubly open-shellnuclei computed from H and H B [ ρ ] are compared inFig. 8. Being based on the minimum of the dHFB TEC,these spectra describe the low-lying part of the ground-state rotational band.Reference results are well reproduced for all one-bodytest densities but ρ sHOSD for which a degrading arises Fig. 9: Same as the bottom panel of Fig. 2 for J Π =0 + , + , + PHFB states.with increasing mass. Even in nuclei for which absolutePHFB energies were not converged yet with respect to e (e.g. Mg and Ne), energy differences are fullyconsistent with the reference values.As for quantitative measures, systematic results are re-ported on in Tab. 4. Reference excitation energies arereproduced to better than 2% throughout the wholepanel for ρ sHF(B) , ρ PHFB and ρ PGCM , which amountsto making errors of the order of a few tens of keVs.This is obviously negligible compared to other sourcesof uncertainties in state-of-the-art ab initio calculations.While this outcome further demonstrates the robustnessof the approximation method, the 10% average error ob-tained for ρ sHOSD underlines the fact that the employedone-body density matrix must be realistic enough todeliver high accuracy results. Given that the purposeof ab initio PHFB (and PGCM below) calculations isto access excitation energies and not absolute ones, onecan be fully satisfied with the performances of H B [ ρ ]in the present context. While PHFB calculations already provide a good testwhenever the system is rigid with respect to collectivevariables, the PGCM opens the way to the wider classof so-called soft nuclei. More generally, it permits toinclude static correlations induced by shape fluctuationsand to access associated vibrational excitations.Presently, PGCM calculations of Ne and Ne alongthe axial quadrupole coordinate are performed. In orderto obtain a first indication of the performance of H B [ ρ ],Fig. 9 extends the study performed at the dHFB levelin Sec. 4.3.1 by displaying the error obtained for the Fig. 10: Low-lying part of the PGCM ground-state rota-tional band of Ne. Reference results calculated from H are compared to those computed from H B [ ρ ] using var-ious test one-body density matrices. Each energy levelis displayed along with the magnetic dipole (below) andelectric quadrupole (above) moments of the associatedstate. B(E2) transitions strengths are displayed usingred arrows. Calculations are performed with e max = 8, e = 10 and λ srg = 1 .
88 fm − .TEC of J Π = 0 + , + , + PHFB energies for the varioustest one-body density matrices. The J Π projected TECconstitutes the diagonal part of the Hamiltonian matrixat play in the Hill-Wheeler-Griffin secular equation ofthe PGCM calculation. The errors obtained along theprojected TECs are strictly similar to those displayedin Fig. 2 at the dHFB level. This result gives confidenceregarding the quality of the results that can be expectedat the PGCM level.Reference and approximate low-lying PGCM excitationenergies of the ground-state rotational band and as-sociated electromagnetic observables are compared for Ne and Ne in Figs. 10 and 11, respectively. Due tonumerical limitations, only three-body matrix elementsup to e = 10 could be included in the full calcu-lation, hence hindering the convergence in Ne. Still,building on the results reported in Fig. 9 an excellentagreement emerges in both nuclei for PGCM energiesand electromagnetic observables, even more so in Newhere sub-percent accuracy (see Tab. 5) is achieved. Asbefore, a decent but less optimal reproduction of the ref- ≤
30 Mass >
30 Neutron-rich Total ρ sHOSD ρ sHF(B) ρ PHFB ρ PGCM
Table 4: Average error (in %) on PHFB low-lying excitation energies computed from H B [ ρ ] for various sub-categories of nuclei and test one-body density matrices. See Eq. (82) for details on the cost function. Calculationsare performed with e max = 8, e = 12 and λ srg = 1 .
88 fm − . Nucleus Ne NeQuantity Spectrum Observables Spectrum Observables ρ sHOSD ρ sHF(B) ρ PHFB ρ PGCM
Table 5: Average error (in %) on PGCM excitation energies and spectroscopic observables computed from H B [ ρ ] in Ne and Ne for various test one-body density matrices. See Eq. (81d) for details on the cost function. Calculationsare performed with e max = 8, e = 10 and λ srg = 1 .
88 fm − .Fig. 11: Same as Fig. 10 for Ne.erence results is obtained for ρ = ρ sHOSD . The excellentresults obtained for electromagnetic observables testifythe stability of the PGCM wave-functions themselveswith respect to the in-medium approximation of thethree-nucleon interaction. QRPA is a method of choice to study excited states ofboth individual and collective characters, with energiesranging from a few MeV to tens of MeV. The method isappropriate as long as the excited states bear a strongresemblance with the ground-state owing to the har-monic approximation at the heart of the QRPA. In thiscontext, the performance of H B [ ρ ] can be assessed bylooking at, e.g., electromagnetic strength functions. Fig-ures 12 and 13 display the electric isovector dipole ( E H and H B [ ρ sHFB ], for O and Ne respectively. Similar results are obtainedfor the other test one-body density matrices and arereported in Table 6.In O, where dQRPA reduces to sRPA built on top ofa sHF Slater determinant, the difference between thestrength functions are hardly noticeable. As in sHF cal-culations discussed earlier on, this relates to the fact thatthe sRPA error would actually be strictly zero if ρ en-tering H B [ ρ ] were taken as the variational sHF densitymatrix obtained from that approximate Hamiltonian.Because ρ sHF coming from the calculation performedwith H slightly differs from the variational one obtainedfrom H B [ ρ sHF ], the error is not strictly zero but re-mains very tiny. As seen from the bottom panel of Fig. 12the relative error of the E ,
50] MeV.This error essentially relates to the horizontal position of Fig. 12: Electric isovector dipole strength of O as afunction of the excitation energy (upper panel). Refer-ence results calculated from H are compared to thosecomputed from H B [ ρ ] using ρ sHFB . The relative devia-tion from the strength computed with the exact Hamil-tonian is displayed as a function of the excitation energyin the lower panel. Calculations are performed with e max = 8, e = 10 and λ srg = 1 .
88 fm − .Fig. 13: Same as Fig. 12 for Ne. the individual dQRPA modes that, as testified in Tab. 6,are shifted by a tiny amount, i.e. 0.05% on average.Computing the observable total photo-emission crosssection by integrating the differential photo-emissioncross section deduced from the E ρ = ρ sHFB is optimal in thepresent case.Moving to Fig. 13, the dQRPA dipole strength of Neobtained from H B [ ρ sHFB ] at the minimum of the dHFBTEC is also visually very close to the reference one,although slightly deteriorated compared to the O case.Looking at the bottom panel, the situation suddenlyappears less favorable with a relative error at fixedexcitation energies that can reach nearly 30% on theleft side of the giant resonance. This error relates to thefact that approximating H by H B [ ρ ] slightly affectsdHFB quasi-particle energies, inducing in turn a smallshift (1% on average) in the position of the dQRPAeigenmodes. The steep slope of the E ρ = ρ sHFB is notapparent as all test one-body density matrices deliversimilar results.The above results validates the quality and robustnessof H B [ ρ ] in the dQRPA context. Although not shownfor brevity, essentially identical results hold for othermultipolarities of the one-body transition operator. Fur-thermore, the dependence of the results on e is alongthe same line as the one discussed in Sec. 4.4.1.4.6 Optimal one-body density matrixThe in-medium approximation of three-body interac-tions proposed in the present work appears to be veryrobust with respect to the employed symmetry-invariantone-body density matrix. All dHFB, dBMBPT, PHFBPGCM and dQRPA results presented above are of equal(excellent) quality for ρ = ρ sHF(B) , ρ PHFB and ρ PGCM but are systematically deteriorated for the more simplis-tic choice ρ = ρ sHOSD .In this context, it is of interest to better assess thisrobustness and possibly characterize the optimal one-body density matrix to be used in the design of H B [ ρ ].For this purpose, trial (symmetry-invariant) one-bodydensity matrices { ρ sRd } are generated by means of therandom sampling described in App. D. To evaluate the O NeQuantity Excitation energy Total photo-emission cross section Excitation energy Total photo-emission cross section ρ sHOSD ρ sHF(B) ρ PHFB ρ PGCM
Table 6: Average relative error (in %) on dQRPA excitation energies and on the total photo-emission cross sectioncomputed from H B [ ρ ] in O and Ne for various test one-body density matrices. Calculations are performedwith e max = 8, e = 10 and λ srg = 1 .
88 fm − .Fig. 14: Ground-state energy error ∆E BΨ [ ρ ] associated with the use of H B [ ρ ] as a function of the distance (cid:107) ρ − ρ Ψ (cid:107) between the test one-body density matrix ρ and the actual ground-state one ρ Ψ in log-log scale. Left panel: sHF( J Π = 0 + ) solution for O. Right panel: PHFB ( J Π = 0 + ) solution for Ne. Data points are for randomly-sampledand physical test one-body density matrices. The color scale characterizes one-body density matrices’ von Neumannentropy. The dashed-dotted lines denote the cubic envelop extracted from the left panel and reported on the rightpanel. Calculations are performed for e max = 6 and e = 6.corresponding approximation H B [ ρ sRd ], the ground-state energy error ∆E BΨ [ ρ ] ≡ (cid:104) Ψ | H B [ ρ ] | Ψ (cid:105)(cid:104) Ψ | Ψ (cid:105) − (cid:104) Ψ | H | Ψ (cid:105)(cid:104) Ψ | Ψ (cid:105) (31)is considered; see App. D.1 for the working expressionand a related discussion. The error function ∆E BΨ [ ρ sRd ]computed for a large set of randomly generated matricesis shown in Fig. 14 as a function of the distance (cid:107) ρ sRd − ρ Ψ (cid:107) between the trial one-body density matrix and the ground-state one ρ Ψ in the many-body calculation ofinterest. The data points corresponding to the physicalone-body density matrices ( ρ sHOSD , ρ sHF(B) , ρ PHFB and ρ PGCM ) are also displayed to better make sense of theresults obtained so far. In addition to the distance to ρ Ψ , each trial one-body density matrix is characterizedby its von Neumann entropy S [ ρ ] ≡ − Tr ( ρ ln ρ ) , (32) which, in the eigenbasis of ρ with eigenvalues { r a } readsas Shannon’s entropy of information theory S [ ρ ] ≡ − (cid:88) a r a ln r a . (33)In the present context, the size of the entropy essentiallycharacterizes how much the many-body state ρ differsfrom a Slater determinant for which S [ ρ ] = 0, i.e. it is ameasure of many-body correlations.Results for O computed in a small model space at thesHF level are shown in the left panel of Fig. 14 whilethe right panel displays results for Ne computed atthe PHFB level. The sHF calculation of O illustratesthe situation encountered for an uncorrelated state, i.e.the many-body solution | Ψ (cid:105) is nothing but a symmetry-conserving Slater determinant. In this particular case,the ground-state energy error (see Eq. (74)) takes thesimple form ∆E B sHF [ ρ ] = 13! w (3) · (cid:0) ρ − ρ sHF (cid:1) ⊗ (3) (34)and is thus minimal, actually null, for ρ = ρ sHF . The factthat the optimal one-body density matrix is nothingbut the one of the many-body state under scrutinyis confirmed numerically in the left panel of Fig. 14.In absence of genuine correlations, one expects fromEq. (34) that the sampled errors are bounded by acubic envelope in the variable (cid:107) ρ − ρ sHF (cid:107) , which indeedappears clearly in the numerical results. The coefficient(0 .
4) of that cubic envelope extracted from the datais a measure of the employed three-body interactionstrength in the utilized model space.Besides the null error delivered by ρ = ρ sHF , the errorsassociated with the physical one-body density matrices ρ sHOSD and ρ sMBPT are provided on the figure. Com-pared to the full range of sampled one-body densitymatrices ρ sHOSD and ρ sMBPT are rather close to ρ sHF .This is particularly true of ρ sMBPT , which is a sign of theweakly-correlated character of O when eventually go-ing beyond the mean-field. Given the cubic upper-bound,such a proximity between the two density matrices im-plies a tiny error on the energy obtained for ρ = ρ sMBPT .In spite of originating from a Slater determinant andthus sharing the same null entropy as ρ sHF , ρ sHOSD isabout 3 times more distant from it than ρ sMBPT . Inagreement with the cubic law governing the error, plus Given that ρ sHF relates to a Slater determinant with 16particles, the maximum distance is reached for densities as-sociated with Slater determinants obtained by promotingthe 16 particles from hole states into particle states, i.e.Max ρ (cid:107) ρ − ρ sHF (cid:107) = √ ≈ .
7, which is indeed the maximumvalue visible on the left panel of Fig. 14. being located closer to the envelope, the associated er-ror is about 170 times larger. Given the softness of theemployed three-body interaction, ρ sHOSD still providesa small absolute error in the end. Eventually, the sam-pling provides a fair understanding that, as long as thetest one-body density is not too distant from ρ sHF , itsdetailed properties do not matter much and the error isbound to be small.Compared to the previous case, the right panel of Fig. 14allows one to appreciate the qualitatively different situ-ation encountered for a genuinely-correlated state. In-deed, the error ∆E B PHFB [ ρ ] behaves now differently as afunction of the distance (cid:107) ρ − ρ PHFB (cid:107) . As visible fromEq. (74), ∆E B PHFB [ ρ ] contains non-zero constant andlinear terms in addition to the cubic term encounteredin Eq. (34).The constant term delivers the error ∆E B PHFB [ ρ PHFB ] as-sociated with the actual ground-state density, i.e. whensetting ρ = ρ PHFB . The fact that this error is differentfrom zero is a fingerprint of the fact the PHFB ground-state wave-function carries (a least) genuine three-bodycorrelations. The value of the corresponding error ad-ditionally depends on the size of the three-body in-teraction convoluted with the irreducible three-bodydensity matrix (see Eq. (74)). As analyzed in Ref. [4]in connection with the NO2B approximation, a low-scale Hamiltonian makes the energy contribution fromthe residual three-body interaction small. For ( e max =6; e = 6) this error is ∆E B PHFB [ ρ PHFB ] ≈ . e max = 8; e = 12) is 0 . . ρ = ρ PHFB , one can lowerthe error such that a minimum Min ρ ∆E B PHFB [ ρ ] is foundfor (cid:107) ρ − ρ PHFB (cid:107) = few 10 − with a value several timessmaller than ∆E B PHFB [ ρ PHFB ]. Passed the minimum theerror typically increases and is eventually dominatedby the cubic terms at large distances such that thecubic envelope extracted from the left panel becomeseffective.The physical density matrices ρ sHOSD and ρ sHFB arefound right passed the minimum such that their erroris small and in fact similar to the one found at theorigin. The profile of the error as a function of thedistance (cid:107) ρ − ρ PHFB (cid:107) rationalizes the fact that smallerrors can be found over a substantial range of densitymatrices to which the various physical one-body densitymatrices one may typically access all belong. This feature Because of the log-log scale employed, the point at zerodistance associated with ρ PHFB is artificially placed on theleft border of the figure.1 provides practitioners with a significant flexibility as faras the choice of the employed one-body density matrixis concerned. Passed that appropriate interval the errorrapidly increases with the distance, as testified by theuse of ρ sHOSD sitting on the edge of it, such that onemay not be too cavalier either regarding the choice of ρ .4.7 Lessons and perspectives The above results demonstrate the usefulness of theproposed in-medium reduction method of three-body in-teraction operators in nuclear ab initio calculations. Thefact that the method relies on the sole use of a one-bodydensity matrix gives much credit to the simplicity of themethod. Furthermore, the high-quality approximationwas shown to be robust with respect to the employedone-body density matrix, which gives much credit tothe flexibility of the method.These conclusions have been validated for nuclei withclosed and open-shell characters, i.e. displaying weakand strong correlations, for light and mid-mass systemsas well as for stable and exotic isotopes. While convinc-ingly substantiated via the use of the several many-bodymethods and for a large class of observables, a furthervalidation of the quality of the approximation on thebasis of ab initio methods built on different paradigmsare of course welcome in the future.
The independence of the results with respect to a largeclass of one-body density matrices is of prime importancefor practical applications in the future, especially giventhat ab initio calculations aspire to move up the nuclearchart towards heavy, doubly open-shell nuclei. Specif-ically, the high-quality results obtained for ρ = ρ sHFB allow one to build H B [ ρ ] at the sole cost of runningfirst a spherical HFB calculation with full three-bodyforces, thus bypassing the need to run any deformedHFB code followed by projections, which would alreadybe too costly with explicit three-body forces in heavynuclei requiring large values of e . Eventually, theenvisioned working algorithm is1. run a spherical HFB calculation with three-nucleonforces to extract ρ sHFB ,2. build H B [ ρ sHFB ],3. run the many-body method of interest with the two-body Hamiltonian H B [ ρ sHFB ], such that even in (heavy) open-shell nuclei – no two-body density matrix has to be extracted, – no genuine open-shell calculation with an explicitthree-nucleon operator has to be performed. While the method presented in this article relies on theuse of symmetry-invariant one-body densities, which canbe generated only starting from J Π = 0 + states (or a su-perposition of such states), it can also be easily used toconstruct effective k -body interactions for odd-even andodd-odd nuclei. For this purpose, one can employ theone-body density generated in a mean-field calculationof a spherical Bogoliubov vacuum constrained to haveodd-even or odd-odd numbers of particles on average.In the case of odd-even systems, it was demonstrated inRef. [18] that such a vacuum represents a good approx-imation to the true mean-field solution obtained withan odd number parity wave function.Of course, the accuracy of the k -body interactions con-structed in this way will have to be properly checkedbut there is no reason to believe that they would per-form particularly worse than those generated to describeeven-even nuclei. In the longer-term future, the in-medium reduction pro-cedure proposed in the present work can be tested todeal with, e.g., four-body interactions [39,40,41] and/orthree-body nuclear currents [42] at a reduced computa-tional cost.
The present work introduced a novel method to approx-imate n -body operators in terms of k -body ones with k < n . This is highly pertinent to overcome the steepincrease of the computational cost of many-body calcu-lations due to the presence of three-nucleon interactions,especially as ab initio calculations aspire to move toheavier nuclei than presently possible.The main advantages are that the method is accurate,universal, simple and flexible. The universality of themethod not only relates to its applicability to all nuclei,independently of their closed or open-shell character, butalso to its independence with respect to the many-bodymethod eventually used to solve Schr¨odinger’s equation.The simplicity of the method relates to the fact that itrequires the convolution of the, e.g., three-body operator with a sole symmetry-invariant one-body density ma-trix, even in open-shell nuclei. This it at variance withexisting methods that either convolute the three-bodyoperator with one-, two- and three-body density matri-ces in open-shell systems or with a symmetry-breakingone-body density matrix, thus leading to an approximateoperator that explicitly breaks symmetries of the initialHamiltonian. Eventually, the flexibility of the methodrelates to the possibility to use various one-body densitymatrices as an entry. As a matter of fact, the functionalform of the error due to the use of the approximateHamiltonian could be exploited to explain why accurateresults can be obtained for a rather large class of one-body density matrices. Such a flexibility can be exploitedto use a (not too) simple density matrix in practicalcalculations, e.g. the density matrix extracted from aspherical Hartree-Fock-Bogoliubov calculation.Extensive numerical results have demonstrated the highaccuracy of the approach over a large panel of nuclei andobservables. The approximation method is thus readyto be employed in routine ab initio calculations in thefuture. Furthermore, the in-medium reduction procedureis ready to be tested on four-nucleon interactions and/orthree-body nuclear currents in order to deal with themat a reduced computational cost. Acknowledgements
We would like to thank A. Tichai forhis help with the spherical solver and R. Roth for providingus with the interaction matrix element used in the numeri-cal simulations. Calculations were performed by using HPCresources from GENCI-TGCC (Contract No. A009057392).This project has received funding from the European Union’sHorizon 2020 research and innovation programme under theMarie Sk(cid:32)lodowska-Curie grant agreement No. 839847.
A Inverse tensor transformations
Given the operator O , i.e. the set of tensors { o ( n ) ; n =0 , . . . , N } , and the one-body density matrix ρ , a secondset of matrices is introduced through o ( k ) [ ρ ] ≡ N (cid:88) n = k n − k )! o ( n ) · ρ ⊗ ( n − k ) (35)with k = 0 , . . . , N . In a second step, the third set oftensors is defined via˜ o ( n ) [ ρ ] ≡ N (cid:88) l = n ( − l − n ( l − n )! o ( l ) [ ρ ] · ρ ⊗ ( l − n ) (36)with n = 0 , . . . , N . The third set of tensors is now shown to be the same asthe initial one. Noting that o ( l ) [ ρ ] · ρ ⊗ ( l − n ) = N (cid:88) k = l k − l )! (cid:16) o ( k ) · ρ ⊗ ( k − l ) (cid:17) · ρ ⊗ ( l − n ) = N (cid:88) k = l k − l )! o ( k ) · ρ ⊗ ( k − n ) , (37)one obtains˜ o ( n ) [ ρ ] = N (cid:88) l = n N (cid:88) k = l ( − l − n ( l − n )! 1( k − l )! o ( k ) · ρ ⊗ ( k − n ) = N (cid:88) k = n (cid:34) k (cid:88) l = n ( − l − n l − n )!( k − l )! (cid:35) o ( k ) · ρ ⊗ ( k − n ) = N (cid:88) k = n (cid:34) k (cid:88) l = n ( − l − n (cid:18) k − nl − n (cid:19)(cid:35) k − n )! o ( k ) · ρ ⊗ ( k − n ) = N (cid:88) k = n δ kn k − n )! o ( k ) · ρ ⊗ ( k − n ) = o ( n ) , (38)such that each original tensor o ( n ) , and thus the fulloriginal operator O , is recovered through the two-stepprocedure. B PGCM transition density matrix
A workable expression for the transition one-body den-sity matrix between two PGCM states is obtained andeventually reduced to the particular case of presentinterest.B.1 Inputs
B.1.1 Bogoliubov state
The deformed Bogoliubov state | Φ ( q ) (cid:105) characterized bythe collective deformation q is a vacuum for the set ofquasi-particle operators [21] β † k ( q ) ≡ (cid:88) a U ka ( q ) c † a + V ka ( q ) c a , (39a) β k ( q ) ≡ (cid:88) a U ∗ ka ( q ) c a + V ∗ ak ( q ) c † a . (39b) Equation (39) defines the unitary Bogoliubov transfor-mation that is inverted according to c † a = (cid:88) k U ∗ ka ( q ) β † k + V ka ( q ) β k , (40a) c a = (cid:88) k U ka ( q ) β k + V ∗ ak ( q ) β † k . (40b) B.1.2 PGCM state
Given a set of Bogoliubov states {| Φ ( q ) (cid:105)} differing bythe value of the collective deformation parameter q , aPGCM state reads as | Ψ σµ (cid:105) ≡ (cid:90) d qf σµ ( q ) P σ | Φ ( q ) (cid:105) . (41)While µ denotes a principal quantum number, σ ≡ (JM Π NZ) collects the set of symmetry quantum num-bers labelling the many-body state, i.e. the angularmomentum J and its projection M, the parity Π = ± P σ ≡ P JM P N P Z P Π (42)collects the projectors on good symmetry quantum num-bers P JM ≡ (cid:88) K g K ( q ) P JMK ≡ (cid:88) K g K ( q ) 2 J + 116 π (cid:90) [0 , π ] × [0 ,π ] × [0 , π ] d Ω D J ∗ MK ( Ω ) R (cid:126)J ( Ω ) , (43a) P N ≡ π (cid:90) π d ϕ N e iϕ n N R N ( ϕ n ) , (43b) P Z ≡ π (cid:90) π d ϕ Z e iϕ p Z R Z ( ϕ p ) , (43c) P Π ≡ (cid:88) ϕ π =0 ,π e i (1 − Π ) ϕ π Π ( ϕ π ) , (43d)such that | Ψ σµ (cid:105) is an eigenstate of J , J z , N , Z and Π ( π ), where the latter denotes the parity operator. Theunknown coefficients ˜ f σKµq ≡ f σµ ( q ) g K ( q ) are typicallyobtained by solving Hill-Wheeler-Griffin’s equation [21].In Eq. (43), Ω ≡ ( α, β, γ ), ϕ π and ϕ n ( ϕ p ) denote Euler, The covariant indices notation used in this document isextended to U and V matrices, such that creator (annihilator)indices are lowered (raised) using complex conjugation, e.g. U ∗ ka ( q ) = ( U ka ( q )) ∗ . By definition the coefficients ˜ f σKµq are such that the set ofPGCM states {| Ψ σµ (cid:105) ; µ = 1 , , . . . } emerging from a calculationare ortho-normalized. parity and neutron- (proton-) gauge angles, respectively.The rotation operators are given by R (cid:126)J ( Ω ) ≡ e − ıαJ z e − ıβJ y e − ıγJ z , (44a) R N ( ϕ n ) ≡ e − iϕ n N , (44b) R Z ( ϕ p ) ≡ e − iϕ p Z , (44c) Π ( ϕ π ) ≡ e − iϕ π F , (44d)with the one-body operator F ≡ (cid:88) ab f ab c † a c b (45)defined through its matrix elements [43] f ab ≡
12 (1 − π a ) δ ab , (46)where π a denotes the parity of one-body basis states thatare presently assumed to carry a good parity. Last butnot least, D JMK ( Ω ) ≡ (cid:104) JM | R (cid:126)J ( Ω ) | JK (cid:105) defines WignerD-matrices, where | JM (cid:105) denotes a generic eigenstate of J and J z . B.1.3 Off-diagonal one-body density matrix
Given a state | Φ ( q ) (cid:105) defined through the Bogoliubovtransformation ( U ( q ) , V ( q )) and the common index θ ≡ ( Ω, ϕ n , ϕ p , ϕ π ) (47)encompassing all rotation angles, the state obtainedthrough multiple rotations | Φ ( q, θ ) (cid:105) ≡ R (cid:126)J ( Ω ) R N ( ϕ n ) R Z ( ϕ p ) Π ( ϕ π ) | Φ ( q ) (cid:105) , (48)is also a Bogoliubov state whose Bogoliubov transforma-tion ( U ( q, θ ) , V ( q, θ )) can be obtained from ( U ( q ) , V ( q ))and from the characteristics of the rotation operators [21,44].A crucial quantity in terms of which the final resultswill be expressed is the so-called off-diagonal one-bodydensity matrix ρ ( q (cid:48) ; q, θ ) ba ≡ (cid:104) Φ ( q (cid:48) ) | c † a c b | Φ ( q, θ ) (cid:105)(cid:104) Φ ( q (cid:48) ) | Φ ( q, θ ) (cid:105) , (49)which involves two different Bogoliubov states and, assuch, can be computed explicitly from the sole knowledgeof ( U ( q (cid:48) ) , V ( q (cid:48) )) and ( U ( q, θ ) , V ( q, θ )) [21,44]. B.2 DefinitionConsidering two PGCM states, the one-body transition density matrix can now be defined through ρ σ f σ i µ f µ i ba ≡(cid:104) Ψ σ f µ f | c † a c b | Ψ σ i µ i (cid:105) = (cid:90) d q f d q i f σ f ∗ µ f ( q f ) f σ i µ i ( q i ) × (cid:104) Φ ( q f ) | P σ f † c † a c b P σ i | Φ ( q i ) (cid:105) = (cid:90) d q f d q i (cid:88) K f K i ˜ f σ f K f ∗ µ f q f ˜ f σ i K i µ i q i δ ( Π f π a π b ) Π i × δ N f N i δ Z f Z i ρ J f M f K f σ i K i q f q i ba , (50)where ρ J f M f K f σ i K i q f q i ba ≡ (cid:104) Φ ( q f ) | P J f † M f K f c † a c b P J i M i K i P N i P Z i P Π i | Φ ( q i ) (cid:105) , (51)and where the action of P N f P Z f P Π f was easily re-solved.B.3 Simplified expressions B.3.1 Expanding the projectors
In order to evaluate this matrix element, the left angular-momentum projector is expanded according to Eq. (43a)such that Eq. (51) is rewritten as ρ J f M f K f σ i K i q f q i ba = 2J f + 116 π (cid:90) d Ω (cid:88) M (cid:48) D J f M f K f ( Ω ) D J i ∗ M i M (cid:48) ( Ω ) × (cid:104) Φ ( q f ) | c † a [ Ω ] c b [ Ω ] P J i M (cid:48) K i P N i P Z i P Π i | Φ ( q i ) (cid:105) , (52)where the rotated creation and annihilation operatorsare defined as c † a [ Ω ] ≡ R (cid:126)J ( Ω ) † c † a R (cid:126)J ( Ω ) , (53a) c b [ Ω ] ≡ R (cid:126)J ( Ω ) † c b R (cid:126)J ( Ω ) , (53b)and where the identity R (cid:126)J ( Ω ) † P J i M i K i = (cid:88) M (cid:48) D J i ∗ M i M (cid:48) ( Ω ) P J i M (cid:48) K i , (54)has been used. In the present derivation, the density matrix is restricted tobe diagonal in the isospin quantum number, i.e. single-particlestates a and b carry the same isospin projection quantumnumber. B.3.2 Spherical one-body basis
In case one-body basis states carry spherical indices ( a ≡ n a , j a , m a , π a , q a ≡ α a , j a , m a ), the operators c † α a j a m a and ( − m a − j a c α a j a − m a transform like the m tha compo-nent of a rank- j a spherical tensor under the action of SU (2), which leads to c † α a j a m a [ Ω ] ≡ (cid:88) m D j a ∗ m a m ( Ω ) c † α a j a m , (55a) c α b j b m b [ Ω ] ≡ (cid:88) m D j b m b m ( Ω ) c α b j b m . (55b)Consequently, the transition density matrix can be sim-plified as ρ J f M f K f σ i K i q f q i ba = 2J f + 116 π (cid:88) mm (cid:48) (cid:88) M (cid:48) × (cid:18)(cid:90) d ΩD J f M f K f ( Ω ) D J i ∗ M i M (cid:48) ( Ω ) D j b m b m ( Ω ) D j a ∗ m a m (cid:48) ( Ω ) (cid:19) × (cid:104) Φ ( q f ) | c † α a j a m (cid:48) c α b j b m P J i M (cid:48) K i P N i P Z i P Π i | Φ ( q i ) (cid:105) . (56)The integral over Wigner-D matrices is performed an-alytically and generates a sum over Clebsch-Gordancoefficients and 3j-symbols according to116 π (cid:90) d ΩD J f M f K f ( Ω ) D J i ∗ M i M (cid:48) ( Ω ) D j b m b m ( Ω ) D j a ∗ m a m (cid:48) ( Ω )= min( J i + J f ,j a + j b ) (cid:88) λ =max( | J i − J f | , | j a − j b | ) ( − M f − K f + J f − J i + j b − j a × (cid:18) J f J i λ − M f M i M f − M i (cid:19) (cid:18) j b j a λm b − m a M f − M i (cid:19) × ( − m (cid:48) − m a C λ ( m − m (cid:48) ) J f − K f J i ( m − m (cid:48) + K f ) C λ ( m − m (cid:48) ) j b mj a − m (cid:48) . (57)The remaining matrix element in Eq. (56) is easily ob-tained in terms of the off-diagonal one-body densitymatrix defined through Eq. (49) (cid:104) Φ ( q f ) | c † α a j a m (cid:48) c α b j b m P J i M (cid:48) K i P N i P Z i P Π i | Φ ( q i ) (cid:105) = 2J i + 116 π π ) (cid:90) d θD J i ∗ M (cid:48) K i ( Ω ) e iϕ n N i e iϕ p Z i × e i (1 − Π ) ϕ π ρ ( q f ; q i , θ ) α b j b mα a j a m (cid:48) (cid:104) Φ ( q f ) | Φ ( q i , θ ) (cid:105) , (58)knowing that the overlap (cid:104) Φ ( q f ) | Φ ( q i , θ ) (cid:105) between twoarbitrary non-orthogonal Bogoliubov states can be com-puted in several ways [45,46]. B.3.3 Special case of present interest
One is presently interested in the one-body density ma-trix of a J Π = 0 + state. In the above set of equations,it corresponds to setting J i = J f = 0, M i = M f = 0 and Π i = +1. In this case, the triangular inequalitiesencoded in the 3j-symbols impose that λ = 0 , (59a) m a = m b , (59b) j a = j b , (59c) m = m (cid:48) , (59d)such that Eq. (57) becomes116 π (cid:90) d ΩD ( Ω ) D ∗ ( Ω ) D j b m b m ( Ω ) D j a ∗ m a m (cid:48) ( Ω )= δ m a m b δ j a j b δ mm (cid:48) j a + 1 . (60)The fact that the initial and final states are the sameand thus carry the same parity further requires that π a = π b . Eventually, Eq. (56) reduces to ρ + N i Z i q f q i ba ≡ δ j a j b δ m a m b δ π a π b π π ) × (cid:90) d θe iϕ n N i e iϕ p Z i (cid:104) Φ ( q f ) | Φ ( q i , θ ) (cid:105)× j a + 1 j a (cid:88) m = − j a ρ ( q f ; q i , θ ) α b j b mα a j a m . (61)The diagonal character of the one-body density matrixin ( j, m ) and its independence on m is made clear inEq. (61) and ends the derivation. C BMBPT transition density matrix
The BMBPT(2) transition one-body density matrixis presently derived within the frame of a so-called expectation-value many-body scheme rather than withina projective one, i.e. it is computed directly throughEq. (6) for l = 1 and | Θ (cid:105) ≡ | Ψ BMBPT(2) (cid:105) . The derivationcan actually be performed within the larger frame of theBogoliubov configuration-interaction (BCI) formalismsuch that BCI-like expansion coefficients are eventuallyobtained within the frame of BMBPT(2) [47].The present applications are eventually limited to stan-dard MBPT, i.e. calculations are restricted to doublyclosed-shell nuclei for which BMBPT reduces to MBPTon top a J Π = 0 + Slater determinant. C.1 Bogoliubov algebraIn methods based on a Bogoliubov vacuum | Φ ( q ) (cid:105) , thegrand potential , Ω ≡ H − µA , (62)must be used rather than the Hamiltonian to control theaverage particle-number in the system. The many-bodyalgebra is more conveniently worked out by normal-ordering all operators with respect to | Φ ( q ) (cid:105) and byexpressing them in terms of quasi-particle operators.Limiting oneself to two-body operators , the grandpotential is thus written as Ω = Ω ( q ) + Ω ( q ) + Ω ( q ) + Ω ( q )+ Ω ( q ) + Ω ( q ) + Ω ( q ) + Ω ( q ) + Ω ( q ) , where Ω ij ( q ) denotes the normal-ordered componentinvolving i ( j ) quasi-particle creation (annihilation) op-erators, e.g., Ω ( q ) ≡ (cid:88) k k k k Ω k k k k ( q ) β † k ( q ) β † k ( q ) β † k ( q ) β k ( q ) , where the matrix elements are anti-symmetric with re-spect to the exchange of any pair of upper or lowerindices. For more details about the normal orderingprocedure, see Refs. [25,9,26,27,10,13].C.2 BCI stateIn BCI, many-body states are written as a CI-like ex-pansion on top of the (deformed) Bogoliubov vacuum.Presently truncated to single and double excitations,the BCISD ansatz reads as | Ψ ( q ) (cid:105) ≡ (cid:32) (cid:88) k k C k k ( q ) β † k ( q ) β † k ( q )+ (cid:88) k k k k C k k k k ( q ) β † k ( q ) β † k ( q ) β † k ( q ) β † k ( q ) (cid:33) × | Φ ( q ) (cid:105) , (63)where the unknown coefficients, anti-symmetric with re-spect to the exchange of any pair of upper indices, can beobtained by diagonalization of Ω or via BMBPT. In practice two separate Lagrange multipliers µ N and µ Z are introduced to account for neutron and proton chemicalpotentials, such that the average neutron and proton numbersare conserved individually. This is the case when working with two-body forces only orwithin the PNO2B or the presently developed 2B approxima-tion.6
C.3 Expression in quasi-particle space
C.3.1 Definition
Considering two different BCISD states | Ψ i ( q ) (cid:105) and | Ψ f ( q ) (cid:105) , the four transition one-body density matricesdefined in terms of quasi-particle operators are given by ρ fik k ( q ) ≡ (cid:104) Ψ f ( q ) | β † k ( q ) β k ( q ) | Ψ i ( q ) (cid:105) (cid:112) (cid:104) Ψ f ( q ) | Ψ f ( q ) (cid:105)(cid:104) Ψ i ( q ) | Ψ i ( q ) (cid:105) , (64a) κ fik k ( q ) ≡ (cid:104) Ψ f ( q ) | β k ( q ) β k ( q ) | Ψ i ( q ) (cid:105) (cid:112) (cid:104) Ψ f ( q ) | Ψ f ( q ) (cid:105)(cid:104) Ψ i ( q ) | Ψ i ( q ) (cid:105) , (64b) − κ ∗ fik k ( q ) ≡ (cid:104) Ψ f ( q ) | β † k ( q ) β † k ( q ) | Ψ i ( q ) (cid:105) (cid:112) (cid:104) Ψ f ( q ) | Ψ f ( q ) (cid:105)(cid:104) Ψ i ( q ) | Ψ i ( q ) (cid:105) , (64c) − σ ∗ fik k ( q ) ≡ (cid:104) Ψ f ( q ) | β k ( q ) β † k ( q ) | Ψ i ( q ) (cid:105) (cid:112) (cid:104) Ψ f ( q ) | Ψ f ( q ) (cid:105)(cid:104) Ψ i ( q ) | Ψ i ( q ) (cid:105) , (64d)among which the relations κ ∗ fik k ( q ) = (cid:16) κ if k k ( q ) (cid:17) ∗ , (65a) σ ∗ fik k ( q ) = (cid:16) ρ if k k ( q ) (cid:17) ∗ − , (65b)hold. C.3.2 Matrix elements
Starting from Eqs. (63)-(64) and applying Wick’s theo-rem, one obtains ρ fik k ( q ) = 1 (cid:112) (cid:104) Ψ f ( q ) | Ψ f ( q ) (cid:105)(cid:104) Ψ i ( q ) | Ψ i ( q ) (cid:105)× (cid:34) (cid:88) k C f ∗ k k ( q ) C ik k ( q )+ 112 (cid:88) k k k C f ∗ k k k k ( q ) C ik k k k ( q ) (cid:35) , and κ fik k ( q ) = 1 (cid:112) (cid:104) Ψ f ( q ) | Ψ f ( q ) (cid:105)(cid:104) Ψ i ( q ) | Ψ i ( q ) (cid:105)× (cid:20) C ik k ( q )+ 14 (cid:88) k k C f ∗ k k ( q ) C ik k k k ( q ) (cid:35) . The expressions of κ fi ∗ and σ ∗ fi are then deduced viaEq. (65) whereas the norm entering the denominators of the transition one-body density matrices reads, e.g.,as (cid:104) Ψ i ( q ) | Ψ i ( q ) (cid:105) ≡ (cid:88) k k C i ∗ k k ( q ) C ik k ( q )+ (cid:88) k k k k C i ∗ k k k k ( q ) C ik k k k ( q ) . C.4 Expression in one-particle spaceInserting the inverse Bogoliubov transformation (Eq. (40)),the normal one-body density matrix expressed in termsof particle operators is obtained under the form ρ fiba ( q ) ≡ (cid:104) Ψ f ( q ) | c † a c b | Ψ i ( q ) (cid:105) (cid:112) (cid:104) Ψ f ( q ) | Ψ f ( q ) (cid:105)(cid:104) Ψ i ( q ) | Ψ i ( q ) (cid:105) = (cid:88) k k (cid:104) U k b ( q ) ρ fik k ( q ) U ∗ k a ( q ) − V ∗ bk ( q ) σ ∗ fik k ( q ) V k a ( q ) − V ∗ bk ( q ) κ ∗ fik k ( q ) U ∗ k a ( q )+ U k b ( q ) κ fik k ( q ) V k a ( q ) (cid:105) . (66)C.5 One-body density matrixIn the present work, one is interested in the case where | Ψ i ( q ) (cid:105) = | Ψ f ( q ) (cid:105) ≡ | Ψ ( q ) (cid:105) such that Eq. (66) reducesto ρ Ψ ba ( q ) = (cid:88) k k (cid:2) V ∗ bk ( q ) V k a ( q )+ U k b ( q ) ρ Ψ k k ( q ) U ∗ k a ( q ) − (cid:16) V k b ( q ) ρ Ψ k k ( q ) V ∗ ak ( q ) (cid:17) ∗ − (cid:16) V k b ( q ) κ Ψ k k ( q ) U ∗ k a ( q ) (cid:17) ∗ + U k b ( q ) κ Ψ k k ( q ) V k a ( q ) (cid:105) . (67)When time-reversal symmetry is preserved, the one-bodydensity matrix can be chosen to be real such that the final expression reads as ρ Ψ ba ( q ) = ρ Φba ( q )+ (cid:88) k k (cid:104) U k b ( q ) ρ Ψ k k ( q ) U ∗ k a ( q ) − V k b ( q ) ρ Ψ k k ( q ) V ∗ ak ( q ) − V k b ( q ) κ Ψ k k ( q ) U ∗ k a ( q )+ U k b ( q ) κ Ψ k k ( q ) V k a ( q ) (cid:35) , (68)where ρ Φ ( q ) denotes the one-body density matrix of thereference state | Φ ( q ) (cid:105) such that the additional termsrelate to BCISD corrections on top of it.C.6 BMBPT coefficientsThe coefficients of the BCISD state obtained at first-and second-order in BMBPT are now explicitly pro-vided [47]. Given that the present application is limitedto a J Π = 0 + Slater determinant reference state | Φ (cid:105) ,the resulting one-body density matrix ρ Ψ ba is actuallyproportional to δ j a j b and relates to a many-body statethat is an eigenvector of the particle-number opera-tor. C.6.1 Partitioning
To formulate BMBPT with respect to the Bogoliubovreference state | Φ ( q ) (cid:105) , the grand potential is split intoan unperturbed part Ω and a residual part Ω [25,9,26] Ω = Ω ( q ) + Ω ( q ) , (69)such that Ω ( q ) ≡ Ω ( q ) + ¯ Ω ( q ) ,Ω ( q ) ≡ Ω ( q ) + ˘ Ω ( q ) + Ω ( q )+ Ω ( q ) + Ω ( q ) + Ω ( q ) + Ω ( q ) + Ω ( q ) , with ˘ Ω ( q ) ≡ Ω ( q ) − ¯ Ω ( q ) and where the normal-ordered one-body part of Ω ( q ) is diagonal, i.e.,¯ Ω ( q ) ≡ (cid:88) k E k ( q ) β † k ( q ) β k ( q ) , with E k ( q ) > k . C.6.2 First-order correction
At first order in BMBPT, singles and doubles BCI-likecoefficients read as C (1) k k ( q ) ≡ − Ω k k ( q ) E k k ( q ) , (70a) C (1) k k k k ( q ) ≡ − Ω k k k k ( q ) E k k k k ( q ) , (70b)where E k k ... ( q ) ≡ E k ( q ) + E k ( q ) + . . . . C.6.3 Second-order correction
At second order in BMBPT, singles and doubles BCI-likecoefficients read C (2) k k ( q ) ≡ P ( k /k ) (cid:88) k k k C (1) k k k k ( q ) Ω k k k k ( q ) E k k ( q )+ 12 (cid:88) k k C (1) k k k k ( q ) Ω k k ( q ) E k k ( q )+ 12 (cid:88) k k C (1) k k ( q ) Ω k k k k ( q ) E k k ( q )+ P ( k /k ) (cid:88) k C (1) k k ( q ) ˘ Ω k k ( q ) E k k ( q ) , (71a) C (2) k k k k ( q ) ≡ P ( k k /k k ) (cid:88) k k C (1) k k k k ( q ) Ω k k k k ( q ) E k k k k ( q )+ P ( k /k k k ) (cid:88) k C (1) k k k k ( q ) ˘ Ω k k ( q ) E k k k k ( q )+ P ( k /k k k ) (cid:88) k C (1) k k ( q ) Ω k k k k ( q ) E k k k k ( q )+ P ( k k /k k ) (cid:104) C (1) k k ( q ) C (1) k k ( q ) (cid:105) , (71b)where anti-symmetrizing operators P ( · · · / · · · ) are ex-pressed in terms of the of the permutation operator as P ( k /k ) ≡ − P k k , (72a) P ( k /k k k ) ≡ − P k k − P k k − P k k , (72b) P ( k k /k k ) ≡ − P k k − P k k − P k k − P k k + P k k P k k . (72c) D Error-function sampling
D.1 Error functionThe error function introduced in Eq. (31) can be inter-preted as the first-order correction to the energy due tothe perturbation δH [ ρ ] ≡ H B [ ρ ] − H . Knowing thatfor a generic n -body operator O ( nn ) (cid:104) Ψ | O ( nn ) | Ψ (cid:105)(cid:104) Ψ | Ψ (cid:105) = 1 n ! 1 n ! (cid:88) a ··· a i b ··· b i o a ··· a n b ··· b n (cid:104) Ψ | A a ··· a n b ··· b n | Ψ (cid:105)(cid:104) Ψ | Ψ (cid:105) = (cid:18) n ! (cid:19) o ( n ) · ρ ( n ) Ψ , (73)Eq. (31) can be written as ∆E BΨ [ ρ ] = − (cid:18) (cid:19) w (3) · ρ (3) Ψ + (cid:18) (cid:19) (cid:16) w (3) · ρ (2) Ψ (cid:17) · ρ − (cid:16) w (3) · ρ Ψ (cid:17) · ρ ⊗ (2) + 13! w (3) · ρ ⊗ (3) = − (cid:18) (cid:19) w (3) · λ (3) Ψ + (cid:18) (cid:19) (cid:16) w (3) · λ (2) Ψ (cid:17) · (cid:0) ρ − ρ Ψ (cid:1) + 13! w (3) · (cid:0) ρ − ρ Ψ (cid:1) ⊗ (3) (74)where λ ( n ) Ψ denotes the irreducible n -body density ma-trix (or cumulants) [48,49] that, for n ≥
2, encodesgenuine n -body correlations in | Ψ (cid:105) . Whenever | Ψ (cid:105) re-duces to a Slater determinant, one has λ ( n ) Ψ = 0 for n ≥ | Ψ (cid:105) ,2. can be written as a cubic polynomial in the vari-able ρ − ρ Ψ whenever involving irreducible densitymatrices of | Ψ (cid:105) ,3. is zero (and thus minimal in absolute value) for ρ = ρ Ψ whenever | Ψ (cid:105) contains at most genuine 2-body correlations, which is notably the case whenever | Ψ (cid:105) reduces to a Slater determinant,4. is in general non-zero for ρ = ρ Ψ and measuresin that case genuine 3-body correlations encodedinto λ (3) Ψ . While the error does not minimize for ρ = ρ Ψ in general, the fact that ρ = ρ Ψ is the optimalsolution whenever | Ψ (cid:105) contains at most genuine 2-body correlations indicates that the optimal ρ cannotbe very different from ρ Ψ whenever | Ψ (cid:105) is a weaklycorrelated J Π = 0 + state.D.2 Random one-body density matricesThe goal is to sample the error function ∆E BΨ [ ρ ] withinthe space of one-body density matrices { ρ } associatedwith J Π = 0 + states . Thus, a large set of densitymatrices { ρ sRd } is randomly generated in the sHO basis {| a (cid:105) ; a ≡ α a j a m a } under the constraints that ρ ba = ( ρ ab ) ∗ , (75a) ρ ba ≡ δ j a j b δ m a m b δ π a π b (cid:37) α b α a , (75b)Tr ρ = (1) · ρ = A , (75c)0 ≤ diag( ρ ) bb ≤ , ∀ b , (75d)where (1) denotes the identity operator on H andwhere diag( ρ ) gathers the eigenvalues. More specifically,the procedure works as follows1. choice of a reference one-body density matrix ρ ref ,2. diagonalization of ρ ref ρ ref ≡ L T diag(r) L , (76)where L denotes an orthogonal matrix.3. choice of two coefficients α d , α o characterizing theamplitude of the random perturbation to be per-formed next.4. sampling of a random perturbation δr of the diagonalmatrix elements of r verifying (cid:88) a δr a = 0 , (77a) r a + δr a ∈ [0 , , ∀ a , (77b) | δr a | ≤ α d , ∀ a . (77c)5. sampling of a random skew-symmetric matrix δl withall upper-diagonal coefficients chosen via a normaldistribution N (0 , Strickly speaking, and as the procedure detailed in Sec. 2.3makes clear, the one-body density matrix employed in theconstruction of H B [ ρ ] does not have to be actually related toa many-body state, i.e. it does not have to be N-representable .It is at least mandatory to use trial one-body density ma-trices carrying the fingerprints of the symmetry constraintsassociated with a true state in order for H B [ ρ ] to display ap-propriate symmetries, which translates into Eqs. (75a), (75b)and (75c) as far as hermiticity, angular momentum, parityand particle number are concerned.9
6. exponentiation of δl to obtain an orthogonal matrix δL ≡ exp [ α o δl ] . (78)7. computation of the random neighbour of ρ ref ρ sRd ≡ ( LδL ) T diag( r + δr ) LδL . (79)Although the sampling is not uniform, all densities withthe required properties can in principle be obtained viathis method.
E Spherical Hartree-Fock field
Given a test one-body density matrix ρ and consideringthat the many-body state of interest | Ψ (cid:105) is a Slaterdeterminant, the one-body Hamiltonian at play in theHF minimization problem based on H B [ ρ ] is, givenEq. (74), (cid:104) h HF(2 B ) [ ρ Ψ ; ρ ] (cid:105) ab ≡ δ (cid:104) Ψ | H B [ ρ ] | Ψ (cid:105) δ [ ρ Ψ ] ba (80)= δ (cid:104) Ψ | H | Ψ (cid:105) δ [ ρ Ψ ] ba + δ∆E BΨ [ ρ ] δ [ ρ Ψ ] ba = (cid:2) h HF [ ρ Ψ ] (cid:3) ab − (cid:104) w (3) · (cid:0) ρ Ψ − ρ (cid:1) ⊗ (2) (cid:105) ab . In Eq. (80), h HF [ ρ Ψ ] denotes the one-body HF Hamilto-nian obtained from the full H whose associated solutionis ρ sHF . Equation (80) allows one to appreciate the im-plications of using H B [ ρ ] at the sHF level, i.e. in themean-field calculation of a doubly closed-shell nucleussuch as O and Ca. One observes that – in general, the use of H B [ ρ ] generates an additionalterm on top of h HF [ ρ Ψ ], eventually leading to ρ Ψ (cid:54) = ρ sHF and ∆E BΨ [ ρ sHF ] (cid:54) = 0 at convergence, – even when using H B [ ρ sHF ], the additional term dif-fers from zero such that ρ Ψ (cid:54) = ρ sHF and ∆E BΨ [ ρ sHF ] (cid:54) =0 at convergence, – only if one were to set ρ = ρ Ψ in H B [ ρ ] throughoutthe iterative procedure, thus modifying the approxi-mate Hamiltonian along the way, would the correc-tion term vanish in Eq. (80) at convergence and thesHF solution based on H B [ ρ ] be the same as the oneobtained from H . This particular case is equivalentto constructing H B [ ρ ] through Wick’s theorem withrespect to the self-consistent sHF Slater determinantitself and is thus identical to the NO2B procedure,which indeed does not lead to any approximation atthe HF level. F Measure of the systematic deviations
In order to assess quantitatively the errors induced bythe approximation and compare the different effectiveinteractions, measures of the average deviation betweenthe results obtained with H B [ ρ ] and H are introducedfor each method, i.e. r HFB ≡ n data (cid:88) n nuclei (cid:12)(cid:12)(cid:12)(cid:12) E HFB [ ρ ] − E HFB E HFB (cid:12)(cid:12)(cid:12)(cid:12) , (81a) r PHFB ≡ n data (cid:88) n nuclei (cid:88) J (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) E J Π PHFB [ ρ ] − E J Π PHFB E J Π PHFB (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (81b) r BMBPT ≡ n data (cid:88) n nuclei (cid:12)(cid:12)(cid:12)(cid:12) E BMBPT [ ρ ] − E PNO2BBMBPT E PNO2BBMBPT (cid:12)(cid:12)(cid:12)(cid:12) , (81c) r PGCM ≡ n data (cid:88) n nuclei (cid:88) J (cid:88) O (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) O J Π PGCM [ ρ ] − O J Π PGCM O J Π PGCM (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (81d)where in practice states up to J = 6 are taken intoaccount when applicable and where O denotes any ob-servable computed within the PGCM formalism (energy,electric and magnetic moments, transitions and radii).In each case, n data denotes the number of terms in thesum(s).The deviation on PHFB excitation energies is given bythe formula r PHFB-S ≡ n data (cid:88) n nuclei (cid:88) J (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) δE J Π PHFB [ ρ ] − δE J Π PHFB δE J Π PHFB (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (82) References
1. R. Roth, S. Binder, K. Vobig, A. Calci, J. Langhammer,P. Navr´atil, Ab Initio Calculations of Medium-Mass Nu-clei with Normal-Ordered Chiral NN+3N Interactions,Physical Review Letters 109 (2012) 052501.2. E. Gebrerufael, A. Calci, R. Roth, Open-shell nucleiand excited states from multireference normal-orderedHamiltonians, Phys. Rev. C 93 (3) (2016) 031301. doi:10.1103/PhysRevC.93.031301 .3. 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 .4. A. Dyhdalo, S. K. Bogner, R. J. Furnstahl, Estimates andpower counting in uniform nuclear matter with softenedinteractions, Phys. Rev. C 96 (5) (2017) 054005. doi:10.1103/PhysRevC.96.054005 .5. G. C. Wick, The Evaluation of the Collision Matrix, Phys.Rev. 80 (1950) 268. doi:10.1103/PhysRev.80.268 .06. V. Som`a, T. Duguet, C. Barbieri, Ab initio self-consistentGorkov-Green’s function calculations of semi-magic nuclei:Formalism at second order with a two-nucleon interaction,Phys. Rev. C 84 (2011) 064317. doi:10.1103/PhysRevC.84.064317 .7. V. Som`a, C. Barbieri, T. Duguet, P. Navr´atil, Movingaway from singly-magic nuclei with Gorkov Green’s func-tion theory arXiv:2009.01829 .8. A. Signoracci, T. Duguet, G. Hagen, G. Jansen, Ab initioBogoliubov coupled cluster theory for open-shell nuclei,Phys. Rev. C 91 (6) (2015) 064320. doi:10.1103/PhysRevC.91.064320 .9. A. Tichai, P. Arthuis, T. Duguet, H. Hergert, V. Som`a,R. Roth, Bogoliubov Many-Body Perturbation Theoryfor Open-Shell Nuclei, Phys. Lett. B 786 (2018) 195. doi:10.1016/j.physletb.2018.09.044 .10. A. Tichai, R. Roth, T. Duguet, Many-body perturbationtheories for finite nuclei, Front. Phys. 8 (2020) 164. doi:10.3389/fphy.2020.00164 .11. J. Ripoche, A. Tichai, T. Duguet, Normal-orderedk-body approximation in particle-number-breakingtheories, The European Physical Journal A 56 (2). doi:10.1140/epja/s10050-020-00045-8 .URL http://dx.doi.org/10.1140/epja/s10050-020-00045-8
12. S. J. Novario, G. Hagen, G. R. Jansen, T. Pa-penbrock, Charge radii of exotic neon and mag-nesium isotopes, Phys. Rev. C 102 (2020) 051303. doi:10.1103/PhysRevC.102.051303 .URL https://link.aps.org/doi/10.1103/PhysRevC.102.051303
13. M. Frosini, T. Duguet, J.-P. Ebran, R. Roth, V. Som`a,A. Tichai, unpublished (2020).14. H. Hergert, S. K. Bogner, T. D. Morris, A. Schwenk,K. Tsukiyama, The In-Medium Similarity Renormaliza-tion Group: A Novel Ab Initio Method for Nuclei, Phys.Rept. 621 (2016) 165. doi:10.1016/j.physrep.2015.12.007 .15. J. Yao, B. Bally, J. Engel, R. Wirth, T. Rodr´ıguez, H. Herg-ert, Ab Initio Treatment of Collective Correlations andthe Neutrinoless Double Beta Decay of Ca, Phys. Rev.Lett. 124 (23) (2020) 232501. doi:10.1103/PhysRevLett.124.232501 .16. W. Kutzelnigg, D. Mukherjee, Normal order and extendedwick theorem for a multiconfiguration reference wave func-tion, J. Chem. Phys. 107 (1997) 432.17. A. Carbone, A. Cipollone, C. Barbieri, A. Rios, A. Polls,Self-consistent Green’s functions formalism with three-body interactions, Phys. Rev. C 88 (2013) 054326. doi:10.1103/PhysRevC.88.054326 .URL https://link.aps.org/doi/10.1103/PhysRevC.88.054326
18. T. Duguet, P. Bonche, P.-H. Heenen, J. Meyer,Pairing correlations. i. description of odd nuclei inmean-field theories, Phys. Rev. C 65 (2001) 014310. doi:10.1103/PhysRevC.65.014310 .URL https://link.aps.org/doi/10.1103/PhysRevC.65.014310
19. S. Perez-Martin, L. M. Robledo, Microscopic justificationof the equal filling approximation, Phys. Rev. C 78 (2008)014304. doi:10.1103/PhysRevC.78.014304 .URL https://link.aps.org/doi/10.1103/PhysRevC.78.014304
20. L. Kong, M. Nooijen, D. Mukherjee, J. Chem. Phys. 132(2010) 234107.21. P. Ring, P. Schuck, The Nuclear Many-Body Problem,Springer Verlag, New York, 1980. 22. B. Bally, M. Bender, Projection on particle number andangular momentum: Example of triaxial bogoliubovquasiparticle states, Phys. Rev. C 103 (2021) 024315. doi:10.1103/PhysRevC.103.024315 .URL https://link.aps.org/doi/10.1103/PhysRevC.103.024315
23. T. Nakatsukasa, T. Inakura, K. Yabana, Finite am-plitude method for the solution of the random-phaseapproximation, Phys. Rev. C 76 (2007) 024318. doi:10.1103/PhysRevC.76.024318 .URL https://link.aps.org/doi/10.1103/PhysRevC.76.024318
24. Y. Beaujeault-Taudi`ere, M. Frosini, J.-P. Ebran, V. Som`a,R. Roth, T. Duguet, unpublished (2021).25. T. Duguet, A. Signoracci, Symmetry broken and restoredcoupled-cluster theory. II. Global gauge symmetry andparticle number, J. Phys. G 44 (1) (2017) 015103, [Er-ratum: J.Phys.G 44, 049601 (2017)]. arXiv:1512.02878 , doi:10.1088/0954-3899/44/1/015103 .26. P. Arthuis, T. Duguet, A. Tichai, R.-D. Lasseri, J.-P.Ebran, ADG: Automated generation and evaluation ofmany-body diagrams I. Bogoliubov many-body pertur-bation theory, Comput. Phys. Commun. 240 (2019) 202. doi:10.1016/j.cpc.2018.11.023 .27. P. Demol, M. Frosini, A. Tichai, V. Som`a, T. Duguet,Bogoliubov many-body perturbation theory under con-straint, Annals Phys. 424 (2021) 168358. doi:10.1016/j.aop.2020.168358 .28. B. Jancovici, D. H. Schiff, Nucl. Phys. 58 (1964) 678.29. E. Khan, N. Sandulescu, N. Van Giai, M. Grasso, Two-neutron transfer in nuclei close to the drip line, Phys. Rev.C 69 (2004) 014314. doi:10.1103/PhysRevC.69.014314 .URL https://link.aps.org/doi/10.1103/PhysRevC.69.014314
30. B. Avez, C. Simenel, P. Chomaz, Pairing vibra-tions study with the time-dependent hartree-fock-bogoliubov theory, Phys. Rev. C 78 (2008) 044318. doi:10.1103/PhysRevC.78.044318 .URL https://link.aps.org/doi/10.1103/PhysRevC.78.044318
31. J. Yao, J. Engel, L. Wang, C. Jiao, H. Hergert, Generator-coordinate reference states for spectra and 0 νββ decayin the in-medium similarity renormalization group, Phys.Rev. C 98 (5) (2018) 054311. doi:10.1103/PhysRevC.98.054311 .32. T. Duguet, B. Bally, A. Tichai, Zero-pairing limit ofHartree-Fock-Bogoliubov reference states, Phys. Rev. C102 (5) (2020) 054320. arXiv:2006.02871 , doi:10.1103/PhysRevC.102.054320 .33. M. Strayer, W. Bassichis, A. Kerman, Phys. Rev. C 8(1973) 1269.34. J. Hoppe, A. Tichai, M. Heinz, K. Hebeler, A. Schwenk,Natural orbitals for many-body expansion methods, Phys.Rev. C 103 (1) (2021) 014321. doi:10.1103/PhysRevC.103.014321 .35. D. R. Entem, R. Machleidt, Accurate charge-dependentnucleon-nucleon potential at fourth order of chiralperturbation theory, Phys. Rev. C 68 (2003) 041001. doi:10.1103/PhysRevC.68.041001 .URL https://link.aps.org/doi/10.1103/PhysRevC.68.041001
36. R. Machleidt, D. R. Entem, Chiral effective field theoryand nuclear forces, Physics Reports 503 (1) (2011) 1 – 75. doi:https://doi.org/10.1016/j.physrep.2011.02.001 .URL doi:10.1007/s00601-007-0193-3 .URL http://dx.doi.org/10.1007/s00601-007-0193-3
38. V. Som`a, P. Navr´atil, F. Raimondi, C. Barbieri, T. Duguet,Novel chiral Hamiltonian and observables in light andmedium-mass nuclei, Phys. Rev. C 101 (1) (2020) 014318. doi:10.1103/PhysRevC.101.014318 .39. D. Rozpedzik, J. Golak, R. Skibinski, H. Witala,W. Glockle, E. Epelbaum, A. Nogga, H. Kamada, A Firstestimation of chiral four-nucleon force effects in He-4,Acta Phys. Polon. B 37 (2006) 2889–2904.40. T. Kr¨uger, I. Tews, K. Hebeler, A. Schwenk, Neutronmatter from chiral effective field theory interactions,Phys. Rev. C 88 (2013) 025802. doi:10.1103/PhysRevC.88.025802 .41. N. Kaiser, R. Milkus, Reducible chiral four-body interac-tions in nuclear matter, Eur. Phys. J. A 52 (1) (2016) 4. doi:10.1140/epja/i2016-16004-7 .42. H. Krebs, Nuclear Currents in Chiral Effective Field The-ory, Eur. Phys. J. A 56 (9) (2020) 234. doi:10.1140/epja/s10050-020-00230-9 .43. J. L. Egido, L. M. Robledo, Nucl. Phys. A524 (1991) 65.44. K. Hara, S. Iwasaki, Nucl. Phys. A332 (1979) 61.45. L. M. Robledo, The sign of the overlap of HFB wavefunctions, Phys. Rev. C79 (2009) 021302.46. B. Bally, T. Duguet, Norm overlap between many-bodystates: Uncorrelated overlap between arbitrary Bogoliubovproduct states, Phys. Rev. C 97 (2) (2018) 024304. doi:10.1103/PhysRevC.97.024304 .47. J. Ripoche, Projected bogoliubov many-body pertur-bation theory : Overcoming formal and technical chal-lenges, Ph.D. thesis, th`ese de doctorat dirig´ee par Duguet,Thomas Structure et r´eactions nucl´eaires Universit´e Paris-Saclay (ComUE) 2019 (2019).URL