ADG: Automated generation and evaluation of many-body diagrams III. Bogoliubov in-medium similarity renormalization group formalism
EEur. Phys. J. A manuscript No. (will be inserted by the editor)
ADG: Automated generation and evaluation of many-body diagrams
III. Bogoliubov in-medium similarity renormalization group formalism
A. Tichai , P. Arthuis , H. Hergert , T. Duguet Technische Universität Darmstadt, Department of Physics, 64289 Darmstadt, Germany ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum für Schwerionenforschung GmbH, 64291 Darmstadt, Germany Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany NSCL / FRIB Laboratory and Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824-1321, USA IRFU, CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette, France KU Leuven, Instituut voor Kern- en Stralingsfysica, 3001 Leuven, BelgiumReceived: February 23, 2021 / Accepted: date
Abstract
The goal of the present paper is twofold. First,a novel expansion many-body method applicable to super-fluid open-shell nuclei, the so-called
Bogoliubov in-mediumsimilarity renormalization group (BIMSRG) theory, is for-mulated. This generalization of standard single-referenceIMSRG theory for closed-shell systems parallels the recentextensions of coupled cluster, self-consistent Green’s func-tion or many-body perturbation theory. Within the realm ofIMSRG theories, BIMSRG provides an interesting alternativeto the already existing multi-reference IMSRG (MR-IMSRG)method applicable to open-shell nuclei.The algebraic equations for low-order approximations,i.e., BIMSRG(1) and BIMSRG(2), can be derived manuallywithout much di ffi culty. However, such a methodology be-comes already impractical and error prone for the derivationof the BIMSRG(3) equations, which are eventually needed toreach high accuracy. Based on a diagrammatic formulation ofBIMSRG theory, the second objective of the present paper isthus to describe the third version (v3.0.0) of the ADG code thatautomatically (1) generates all valid BIMSRG(n) diagramsand (2) evaluates their algebraic expressions in a matter ofseconds. This is achieved in such a way that equations caneasily be retrieved for both the flow equation and the Magnusexpansion formulations of BIMSRG.Expanding on this work, the first future objective is to nu-merically implement BIMSRG(2) (eventually BIMSRG(3))equations and perform ab initio calculations of mid-massopen-shell nuclei. a e-mail: [email protected] b e-mail: [email protected] c e-mail: [email protected] d e-mail: [email protected] PROGRAM SUMMARY
Program title:
ADG
Licensing provisions: GNU General Public License Version3 or laterProgramming language: Python 3Repository and DOI: github . com/adgproject/adg DOI: . . Nature of problem: As formal and numerical developmentsin many-body-perturbation-theory-based ab initio meth-ods make higher orders reachable, manually produc-ing and evaluating all the diagrams becomes rapidlyuntractable as both their number and complexity growquickly. Thus, derivations are prone to mistakes andoversights.Solution method: BIMSRG diagrams are encoded as squarematrices known as oriented adjacency matrices ingraph theory, and then turned into graph objects us-ing the
NetworkX package. These objects are usedto eventually evaluate many-body expressions on apurely diagrammatic basis. The new capabilities addto those available in version (v2.0.0) of the code, i.e.the production and evaluation of HF-MBPT, diagonaland o ff -diagonal BMBPT diagrams. The intrinsic cost to solve the many-body Schrödinger equa-tion scales exponentially with the particle number A. Thisposes a great challenge to push ab initio calculations based onessentially exact methods such as configuration interaction [1–3] (CI) or quantum Monte-Carlo (QMC) approaches [4–6]beyond the lightest nuclei. During the last two decades theuse of methods that systematically expand the exact solu-tion with respect to a simple yet su ffi ciently rich A-body a r X i v : . [ nu c l - t h ] F e b reference state has emerged as an invaluable paradigm toextend nuclear many-body theory to heavier masses. Indeed,truncating the expansion provides a solution whose numer-ical cost scales polynomially with A. The most prominentexamples of such methods are many-body perturbation the-ory (MBPT) [7–12], self-consistent Green’s function (SCGF)theory [13–17], coupled-cluster (CC) theory [18–21] andthe in-medium similarity renormalization group (IMSRG)approach [22–28].Expansion methods, whenever based on a symmetry-conserving Slater-determinant reference state, e ffi ciently de-scribe weakly correlated closed-shell nuclei. Contrarily, theyfail to describe strongly-correlated open-shell nuclei due tothe inherent degeneracies of elementary excitations out ofsuch a Slater-determinant reference state. Within the last tenyears, extending these methods to symmetry-breaking refer-ence states has shown to be a powerful idea to systematicallyovercome the limitation to closed-shell systems. In singlyopen-shell nuclei, U(1) global-gauge symmetry associatedwith particle-number conservation is relaxed through the useof a Bogoliubov quasiparticle vacuum as the reference state.In doubly open-shell nuclei, SU(2) rotational symmetry as-sociated with angular-momentum conservation is (further)allowed to break via the use of deformed Slater determinantsor Bogoliubov vacua as reference states. The breaking ofU(1) symmetry, specifically, entails an explicit extension ofthe formalism based on the more general Bogoliubov algebra.This has been successfully achieved in Gorkov’s extension ofSCGF theory (GGF) [29–31], in Bogoliubov coupled-cluster(BCC) theory [32, 33] and in Bogoliubov many-body per-turbation theory (BMBPT) [12, 34–37]. The present worksimilarly generalizes the single-reference IMSRG methodto Bogoliubov in-medium similarity renormalization group (BIMSRG) theory.In recent years, various flavors of the IMSRG paradigmhave led to unprecedented ab initio calculations of atomicnuclei up to one hundred interacting particles [38], of globalbenchmarks [39] or of deformed nuclei [40]. Calculationsof open-shell nuclei have been performed via the multi-reference IMSRG (MR-IMSRG) formulated with respect to acorrelated vacuum, thus requiring a more elaborate formalismthan for the standard, i.e. single-reference, flavors [41–43].The presently developed BIMSRG method provides a simplealternative whose advantages and drawbacks compared toMR-IMSRG can be gauged in future calculations.The development of novel many-body theories applicableto a large set of nuclei and the need for accurate implementa-tions, i.e., for advanced truncation schemes, takes the deriva-tion of the working equations to the edge of what is humanlydoable. The design of BIMSRG at the BIMSRG(3) trunca-tion level, for instance, is a perfect example with a total of82 new diagrams to derive even when exploiting symmetriesto reduce their number. Consequently, the general discussion of the BIMSRG formalism is accompanied here by the useof a new, extended version of the
ADG code [36, 44] that iscapable of performing the sophisticated formal derivations inan automated fashion.The present work is structured as follows. The generalBIMSRG formalism is introduced in Sec. 2 before being spe-cialized to the use of normal-ordered operators with respectto a Bogoliubov vacuum in Sec. 3. Section 4 provides thediagrammatic formulation of the BIMSRG that is used todesign the automated derivation of the working equationsat arbitrary BIMSRG(n) truncation levels through the
ADG code described in Sec. 5. The conclusions and outlook aregiven in Sec. 7, and an appendix provides the details for theBIMSRG(2) approximation. { c † p , c p } of the one-body Hilbert space H is given by H ≡ (cid:88) pq t pq c † p c q + (cid:88) pqrs v pqrs c † p c † q c s c r + (cid:88) pqrstu w pqrstu c † p c † q c † r c u c t c s , (1)where t pq denotes matrix elements of the kinetic energywhereas v pqrs and w pqrstu denote anti-symmetric matrix ele-ments of two- and three-nucleon interactions, respectively.The generalization to Hamiltonians with even higher many-body, e.g. four-body, interactions is straightforward (but te-dious).In the following a constraint term involving the particle-number operator A ≡ (cid:88) p c † p c p , (2)is required to control the average particle-number in themany-body states of interest, hence we work with the grandpotential , Ω ≡ H − µ A , (3)instead of the Hamiltonian. In practice two separate Lagrange multipliers µ N and µ Z are introducedto account for neutron and proton chemical potentials, such that neutronand proton number are conserved individually. β † p ≡ (cid:88) k U pk c † k + V pk c k , (4a) β p ≡ (cid:88) k U ∗ pk c k + V ∗ pk c † k . (4b)They obey standard Fermionic anti-commutation relations { β † p , β † q } = , (5a) { β p , β q } = , (5b) { β p , β † q } = δ pq . (5c)The set of transformation coe ffi cients { U , V } in Eq. (4) definesa unitary Bogoliubov transformation satisfying U † U + V † V = , (6a) UU † + V (cid:63) V T = , (6b) U T V + V T U = , (6c) UV † + V (cid:63) U T = . (6d)The set of quasi-particle operators introduced in Eq. (4) alsodefine a many-body Bogoliubov vacuum | Φ (cid:105) through β p | Φ (cid:105) = , ∀ p . (7)Due to the mixing of single-particle creation and annihilationoperators in Eq. (4), | Φ (cid:105) breaks U(1) global-gauge symmetryassociated with the conservation of particle number, i.e., it isnot an eigenstate of the particle-number operator. In practice,transformation coe ffi cients { U , V } are typically obtained bysolving Hartree-Fock-Bogoliubov (HFB) mean-field equa-tions.Using | Φ (cid:105) as the reference state, excitations are syste-matically obtained via the action of an arbitrary number ofquasi-particle creation operators | Φ pq (cid:105) ≡ β † p β † q | Φ (cid:105) , (8a) | Φ pqrs (cid:105) ≡ β † p β † q β † r β † s | Φ (cid:105) , (8b) ... Combining the reference state with the complete set of itselementary excitations provides a basis of Fock space F . The number of quasi-particle excitations is restricted to be even due tothe conservation of number parity in the present work. Such a constraintshould be relaxed when, e.g., building states associated with a systemcontaining an odd number of particle via the action of an excitationoperator acting on the state describing an even neighbor. Whenever the reference state reduces to a Slater determinant, the many-body states introduced in Eq. (8) themselves reduce to particle-numberconserving Slater determinants obtained via arbitrary numbers of n-particle / p-hole excitations of the reference state. The associated basisof Fock space becomes a direct sum of bases of Hilbert spaces H N associated with definite particle numbers. U ( s )of the grand potential ΩΩ ( s ) ≡ U ( s ) Ω U † ( s ) , (9)that is parameterized by the continuous variable s ∈ R . Thesimilarity transformation (Eq. (9)) can be re-cast into a first-order ordinary di ff erential equation (ODE)dd s Ω ( s ) = [ η ( s ) , Ω ( s )] , (10)involving an anti-Hermitian generator η ( s ) that indirectlyparametrizes the transformation. The ODE on the similaritytransformed grand potential is to be solved with the initialcondition Ω (0) = Ω .The generator is chosen such that the reference state isdecoupled from its quasi-particle excitations at the end of theflowlim s →∞ (cid:104) Φ pq | Ω ( s ) | Φ (cid:105) = , (11a)lim s →∞ (cid:104) Φ pqrs | Ω ( s ) | Φ (cid:105) = . (11b) ... The exact A-body ground-state energy E can then be ob-tained from the vacuum expectation value of the flowinggrand potentiallim s →∞ (cid:104) Φ | Ω ( s ) | Φ (cid:105) = E − µ A . (12)2.4 Magnus formulationA di ff erent way to implement IMSRG methods is to solvefor the unitary transformation itself. In the so-called Mag-nus formulation of the approach [46], the transformation isparameterized as U ( s ) ≡ e M ( s ) , (13)where the Magnus operator M ( s ) satisfies M † ( s ) = − M ( s )and M (0) =
0. One can now solve for M ( s ) by integratingthe ODE [47, 48]dd s M ( s ) = ∞ (cid:88) l = B l l ! ad ( l ) M ( η ) . (14) The specific generators allowing for such a decoupling are not de-tailed in the body of the text. However, such a discussion is providedin Appendix A for the BIMSRG(2) approximation. The present article departs from the conventional notation for theMagnus operator [46–48] to avoid confusion with the grand-canonicalpotential introduced above.
Here, the B l are Bernoulli numbers and ad M ( η ) characterizesthe adjoint action defined recursively throughad (0) M ( η ) ≡ η ( s ) , (15a)ad ( l ) M ( η ) ≡ [ M ( s ) , ad ( l − M ( η )] . (15b)Once the Magnus operator is obtained, any arbitrary op-erator O , including Ω , can be transformed consistently byusing the Baker-Campbell-Hausdor ff (BCH) formula O ( s ) ≡ e M ( s ) Oe − M ( s ) = ∞ (cid:88) k = k ! ad ( k ) M ( O ) . (16)2.5 Particle-number adjustmentBecause U(1) is a symmetry of the nuclear many-body Ha-miltonian, i.e.,[ Ω, A ] = , (17)the identity[ Ω ( s ) , A ( s )] = up to truncations,such thatlim s →∞ (cid:104) Φ | A ( s ) | Φ (cid:105) = A . (19)Two di ffi culties arise to ensure that Eq. (19) is indeedfulfilled. First, for the targeted particle number A to be ob-tained in an exact calculation, the chemical potential µ mustbe fixed such that E − µ A is the lowest of all eigenvalues of Ω over Fock space. This is achievable only if the ground-stateenergy is strictly convex as a function of the particle numberin the neighborhood of A, which is generally but not alwaystrue for atomic nuclei. Second, approximations made duringthe flow are such that the symmetry properties embodied byEqs. (18) and (19) are explicitly violated in practice, i.e., theycould only be ensured by an exact BIMSRG calculation.These two features require a control of Eq. (19) through-out the flow. Following the strategies used in other expansionmethods based on Bogoliubov reference states [32, 35, 37],it can typically be achieved by an adjustment of the averageparticle number of the reference state along the flow. Provid-ing the reference state with an explicit s dependence, oneallows its average particle number to be adapted along theflow (cid:104) Φ ( s ) | A | Φ ( s ) (cid:105) = A aux ( s ) , (20)such that Eq. (19) is eventually fulfilled. Numerically, thiscan be achieved by iteratively solving constrained HFB calcu-lations and BIMSRG flow equations, where the constrainedvalue for the reference-state particle number is estimatedfrom the BIMSRG correction until convergence is achieved. The unitary transformation itself does not fulfill the symmetry suchthat the commutator with the unevolved particle-number operator doesnot vanish for s >
0, i.e., [ Ω ( s ) , A (0)] (cid:44) Because BIMSRG relies on the use of Bogoliubov productstates, algebraic equations are conveniently obtained by for-mulating the approach in terms of normal-ordered operatorsobtained via the use of Wick’s theorem [49] with respect tothe Bogoliubov vacuum | Φ (cid:105) .3.1 Normal-ordered operatorsLet us consider a generic particle-number-conserving opera-tor containing up to N -body terms O ≡ N (cid:88) n = o [2 n ] , (21)where the class o [2 n ] contains a single n -body operator o [2 n ] ≡ o nn = n ! n ! (cid:88) l ... l n o nnl ... l n l n + ... l n c † l . . . c † l n c l n . . . c l n + . (22)In Eq. (22), n -body matrix elements o nnl ... l n l n + ... l n are mode-2 n tensors, i.e., they are data arrays carrying 2 n indices associ-ated with the n particle creation and n annihilation operatorsthey multiply. They are fully anti-symmetric with respect tothe permutation of the n first, resp. n last, indices o nnl ... l n l n + ... l n = (cid:15) ( σ ) o nn σ ( l ... l n | l n + ... l n ) , (23)where (cid:15) ( σ ) refers to the signature of the permutation σ . Thenotation σ ( . . . | . . . ) denotes a separation between the n firstand the n last indices such that permutations are only consid-ered between members of the same group. If O is Hermitian,n-body matrix elements o nnl ... l n l n + ... l n satisfy o nn ∗ l ... l n l n + ... l n = o nnl n + ... l n l ... l n . (24)After normal ordering with respect to | Φ (cid:105) , the operator O can be expressed as O ≡ N (cid:88) n = O [2 n ] , (25)where the class O [2 n ] groups all terms containing a normal-ordered product of 2 n quasiparticle operators. We furtherdi ff erentiate these terms according to the number of quasi-particle creation and annihilation operators they contained,i.e., O [2 n ] ≡ N (cid:88) i , j = i + j = n O i j , (26) where O i j gathers terms with i , resp. j , quasiparticle creation,resp. annihilation, operators and reads O i j ≡ i ! j ! (cid:88) k ... k i + j O i jk ... k i k i + ... k i + j β † k . . . β † k i β k i + j . . . β k i + . (27)Matrix elements O i jk ... k i + j are mode-( i + j ) tensors, i.e., theyare data arrays carrying i + j indices associated with the i ( j )quasi-particle creation (annihilation) operators they multiply.They are fully anti-symmetric with respect to the permutationof the i first, resp. j last, indices O i jk ... k i k i + ... k i + j = (cid:15) ( σ ) O i j σ ( k ... k i | k i + ... k i + j ) , (28)and satisfy O i j ∗ k ... k i k i + ... k i + j = O jik i + ... k i + j k ... k i , (29)if O is Hermitian. Explicit expressions of such anti-symmetrizedmatrix elements up to the class O [6] can be found in Refs. [32,50].All operators at play in the present context, e.g., Ω ( s ), η ( s )and M ( s ), can now be written in such a normal-ordered form.To give one explicit example, the original grand potentialoriginating from the Hamiltonian introduced in Eq. (1) reads Ω ≡ Ω ( s = = Ω [0] + Ω [2] + Ω [4] + Ω [6] ≡ Ω + Ω + Ω + Ω + Ω + Ω + Ω + Ω + Ω + Ω + Ω + Ω + Ω + Ω + Ω + Ω , (30)i.e., it contains up to Ω [6] terms. By virtue of the normalordering, one trivially finds that (cid:104) Φ | Ω | Φ (cid:105)(cid:104) Φ | Φ (cid:105) = Ω . (31)Whenever | Φ (cid:105) corresponds to the HFB vacuum, it is notcoupled via Ω to two-quasi-particle excitations by virtue ofBrillouin’s theorem, i.e., Ω k k = Ω k k = , ∀ ( k , k ) . (32)3.2 Flow equationThe BIMSRG flow equation can now be written as a systemof coupled ODEs for the i j -components of the derivative ofthe grand potential, i.e., (cid:32) d Ω ds (cid:33) i j ( s ) ≡ [ η ( s ) , Ω ( s )] i j , (33) where [ , ] i j denotes the i j -component of the normal-orderedoperator resulting from the commutator. Since the Hermitianadjoints of Ω ( s ) and η ( s ) satisfy( Ω i j ( s )) † = Ω ji ( s ) , (34a)( η i j ( s )) † = − η ji ( s ) , (34b)the derivative of Ω is itself Hermitian, as expected. Conse-quently, the system of ODEs contains redundant equationsand one only needs to explicitly solve for i ≥ j before deduc-ing the results for i < j .3.3 Magnus formulationSimilarly, the ODE for the Magnus operator is transformedinto a set of coupled ODEs (cid:32) dMds (cid:33) i j ( s ) ≡ ∞ (cid:88) l = B l l ! (ad ( l ) M ( η )) i j , (35)where (ad ( l ) M ( η )) i j denotes the i j -component of the normal-ordered operator resulting from the l -fold nested commu-tators. Eventually, the i j -component of the normal-orderedtransformed operator associated with any observable O ofinterest is obtained through O i j ( s ) = ∞ (cid:88) k = k ! (ad ( k ) M ( O )) i j . (36)As before, the (anti-)Hermiticity of O ( s ) ( M ( s )) is used toavoid redundancies.3.4 Approximation schemesThe BIMSRG formalism discussed so far defines a way toexactly solve the A-body Schrödinger equation to accessnuclear ground-states observables. As such, it is compu-tationally intractable. For actual applications, appropriatetruncation schemes have to be devised to make calculationsfeasible.The ODE in Eq. (10) is characterized by the commu-tator [ η ( s ) , Ω ( s )] appearing on the right-hand side. Giventwo operators A and B containing terms up to class A [2 N A ] and B [2 N B ] , respectively, the commutator C ≡ [ A , B ] containsterms up to class C [2 N C ] with N C = N A + N B −
1. Consequently,the maximum class of Ω ( s ), and thus of η ( s ), increases ateach propagation step of the ODE and the system of coupledODEs (33) grows rapidly. Carrying the associated operators(or, equivalently, the coe ffi cient tensors defining them) andsolving the growing numbers of coupled ODEs becomesquickly intractable.Adapting the truncation scheme of standard IMSRG cal-culations and denoting it as BIMSRG(n), only terms O [2 m ] with m ≤ n in all the operators at play will be retainedthroughout the flow. The rationale behind this truncation canbe justified via perturbation theory arguments and relies onthe fact that the omitted terms are (naively) of higher orderthan the included ones. Of course, whenever the accuracy athand is insu ffi cient, one must go to the next truncation level.At present, the typical working truncation scheme employedin standard IMSRG calculations is IMSRG(2) but there is aconcerted e ff ort to move towards IMSRG(3).A similar truncation scheme can be applied to the Magnusformulation, i.e., operators in Eqs. (14) and (16), includingthose that appear in the construction of (ad ( l ) M ( η )) i j , are re-stricted to contributions O [2 m ] with m ≤ n . The summationson the right-hand sides of Eqs. (14) and (16) are cut o ff oncethe size of the terms drops below a desired numerical thresh-old [46] . The manual derivation of the BIMSRG tensor networks basedon a straightforward application of Wick’s theorem is alreadyquite tedious at the BIMSRG(2) truncation level. Given thatthe goal is to employ more refined truncation schemes, even-tually, a systematic and less error-prone methodology forderiving the BIMSRG(n) working equations is highly desir-able. To set the stage for the code that will automate this thederivation (see Sec. 5), we now introduce a diagrammaticmethod for evaluating BIMSRG commutators.4.1 Fundamental commutatorThe IMSRG method, whether formulated through Eq. (33)or through Eqs. (35)-(36), entirely relies on the (repeated)computation of the elementary commutator C ≡ [ A , B ] = AB − BA , (37)where A and B denote operator whose normal-ordered contri-butions A [2 n ] ( B [2 n ] ) range from n = n = N A ( n = N B ).More specifically, the goal is to compute each normal-ordered component of the operator CC i j ≡ [ A , B ] i j , (38)whose maximum class C N C is such that N C = N A + N B − .Thus, the objective is to work out the general case charac-terized by the triplet ( N A , N B ; N C ). As for the application ofpresent interest, the derivation of the BIMSRG(n) equationsrelies on restricting the procedure to N A = N B = n and tofurther limiting the output terms to C [2 m ] with m ≤ n , i.e., it Whenever N A ( N B ) is equal to 1, one has N C = N B ( N C = N A ); i.e. itis the only situation where the maximum class of C is not greater thanthose of A and B . corresponds to working out the restricted case characterizedby ( N A , N B ; N C ) = ( n , n ; n ). With this at hand, and at theprice of (repeatedly) replacing A and B with the operators ofinterest, the BIMSRG(n) algebraic equations are obtained.4.2 Rationale of the approachIn principle, the computation of Eqs. (37)-(38) is straightfor-ward: Indeed, it boils down to the application of the standardtime-independent Wick’s theorem to a product of two normal-ordered operators . Furthermore, the fact that one is actuallyinterested in the commutator of the two operators allowsone to reduce the computation to so-called connected terms ,i.e., to terms where at least one elementary contraction oc-curs between the two operators in the application of Wick’stheorem.A straightforward application of Wick’s theorem quicklybecomes cumbersome as N A and N B grow. Furthermore, itwould involve repetition since many algebraic contributionsare in fact identical. By identifying the corresponding pat-terns, one can design a diagrammatic representation of thevarious normal-ordered contributions to the commutator andevaluate their algebraic expressions such that a single dia-gram captures all these identical contributions at once. Thisapproach eventually allows us to design an optimal code thatautomatically generates all needed diagrams and algebraicexpressions in a matter of seconds. This will constitute thepresent extension of the ADG code [36, 44]. The diagram-matic framework introduced below shares similarities withthe e ff ective Hamiltonian diagrams occuring in BCC theorysince, beyond closed vacuum-to-vacuum diagrams, linkeddiagrams with external legs need to be considered [51].4.3 Diagrammatic representationThe diagrammatic representation of an arbitrary operator O is given by a sum of Hugenholtz vertices denoting its normal-ordered contributions O i j . In Fig. 1, such a representation isdisplayed for an operator O containing classes of terms upto O [6] . In the following, d O ij ≡ i + j describes the degree ofa given vertex. All terms O i j belonging to the same class O [2 n ] share the same degree d O ij = n .In the diagrammatic method presented below, a graphicalrepresentation of each of the three operators A , B and C occuring in Eq. (37) is needed. The corresponding verticesare introduced in Fig. 2. Rigorously speaking, one is dealing with the product of two operatorsthat are each made out of a sum of normal-ordered contributions. In the present paper, we restrict ourselves to operator of even vertexdegree, as is the case when working with operators that conserve thenumber of particles. The formalism and code described in the followingcan be straightforwardly extended to incorporate operators with odd-degree vertices. O [0] = O O [2] = O + O + O O [4] = O + O + O + O + O O [6] = O + O + O + O + O + O + O Fig. 1
Diagrammatic representation of all normal-ordered contribu-tions to an operator O containing classes of terms up to O [6] . Each termdenotes a fully anti-symmetric Hugenholtz vertex. The convention isthat reading an operator from left to right corresponds to an up-downreading of the associated diagram. A B C
Fig. 2
Vertices associated with operators A , B and C . The diagrammatic representation of O i j involves the ver-tex itself along with i ( j ) lines traveling out of (into) the vertexand representing quasi-particle creation (annihilation) opera-tors. The diagram is meant to be labelled in order to eventu-ally translate it into an algebraic expression. The canonicalrule to do so consists of first assigning indices k . . . k i con-secutively from the leftmost to the rightmost line above thevertex, while k i + . . . k i + j must be similarly assigned consec-utively to lines below the vertex. Next, the anti-symmetricmatrix element O i jk ... k i k i + ... k i + j is to be assigned to the vertex.This canonical labelling is illustrated in Fig. 3 for the diagramrepresenting the operator O . k k k k + O k k k k = k k k k + O k k k k Fig. 3
Left: canonical labelling of the legs and associated vertex inthe diagrammatic representation of the operator O . Right: sign rule toapply when departing from the canonical labelling through the crossingof the legs associated with the two creation operators. R − + k k = k k = k Fig. 4
Diagrammatic representation of the single non-zero elementarycontraction. The convention is that the left-to-right reading of a matrixelement corresponds to the up-down reading of the diagram.
In case the left-right labelling of the incoming and / or out-going legs does not match the canonical ordering employedin the amplitude O i jk ... k i k i + ... k i + j , the latter must be multipliedwith the sign ( − l p , where l p denotes the signature of thepermutation of the incoming and / or outgoing legs necessaryto bring them into the canonical ordering. This additionalrule is illustrated in Fig. 3 for O .Beyond operators, the elementary contractions occurringin the standard Wick’s theorem need to be represented di-agrammatically. Given that the operators are convenientlyexpressed in the quasi-particle basis associated with the Bo-goliubov vacuum, the four elementary contractions take thesimplest possible form R k k = (cid:104) Φ | β † k β k | Φ (cid:105)(cid:104) Φ | Φ (cid:105) (cid:104) Φ | β k β k | Φ (cid:105)(cid:104) Φ | Φ (cid:105)(cid:104) Φ | β † k β † k | Φ (cid:105)(cid:104) Φ | Φ (cid:105) (cid:104) Φ | β k β † k | Φ (cid:105)(cid:104) Φ | Φ (cid:105) ≡ (cid:32) R + − k k R −− k k R ++ k k R − + k k (cid:33) = (cid:32) δ k k (cid:33) , (39)such that the sole non-zero contraction R − + k k = δ k k needs tobe considered. The diagrammatic representation of this con-traction is provided in Fig. 4 and connects two up-going linesassociated with one annihilation and one creation operator,both carrying the same quasi-particle index. For simplicity,one can eventually represent the contraction as a line carryinga single up-going arrow along with one quasi-particle index.The complete set of diagrams that contribute to C = [ A , B ] can now be systematically generated and evaluatedby assembling the building blocks according to a specificset of rules. In principle, both contributions + AB and − BA need to be processed explicitly. However, given that the firstterm is considered for arbitrary values of N A and N B , thediagrams making up the second term can easily be deduced from those contributing to the first one by (i) adding a globalminus sign and (ii) performing a systematic exchange A ↔ B .Thus, the rules to produce and evaluate relevant diagrams areonly specified for + AB below while it is understood that thisterm is to be combined with the diagrams for − BA eventually.As explained below, summing the two terms making up thecommutator leads to the cancellation of a specific categoryof contributions.4.4 Diagram generationThe complete set of relevant diagrams are generated fromthe basic building blocks by applying the following set oftopological rules:1. A diagram contributing to C combines a top vertex A and a bottom vertex B connected by a number n l of ele-mentary contractions R − + . Because A and B are normalordered, contractions starting and ending at the samevertex cannot be of the R − + type and are therefore nec-essarily zero according to Eq. (39). Consequently, eachcontraction necessarily connects A and B .2. The fact that C actually represents the commutator of A and B immediately forces all relevant diagrams to be connected , i.e., in each diagrams A and B are connectedby at least one contraction and n l >
0. Indeed, terms with n l = AB and BA and cancelout. Consequently, contributions involving A and B along with the class of terms C [2( N A + N B )] can be omittedfrom the outset. As a result,(a) Normal-ordered contributions A kl and B mn can belimited to classes d A kl ≡ k + l = , . . . , N A and d B mn ≡ m + n = , . . . , N B , respectively.(b) C receives non-zero contributions C i j for d C ij ≡ i + j = , . . . , N A + N B − C i j , normal-ordered contributions to A and B must be contracted inall possible ways. Not all combinations of A kl and B mn can actually contribute to a given term C i j : The condition i − j = ( k − l ) + ( m − n ) , (40)must be fulfilled, which allows one to discard a largeset of combinations from the outset. Whenever Eq. (40)is satisfied, the number of contractions involved in thediagram denoted by C i j ( kl , mn ) is n l = k + m − i = l + n − j . (41) The fact that BIMSRG(n) equations arise from a commutator ensuresthat the associated results are size-extensive , i.e., that truncation errorsfrom omitting diagrams scale linearly with the system size.
4. Specific properties of A and B typically translate intoproperties of C such that only a subset of the normal-ordered contributions C i j has to be computed explicitly.For example, if A is anti-Hermitian and B is Hermitian, C is Hermitian and it is su ffi cient to only evaluate thesubset of terms { C i j , i ≥ j } .5. A convenient way to classify diagrams is to introduce themaximal vertex degree of a diagram, defined as d max ≡ max( d A kl , d B mn , d C ij ). In particular, the diagrams arisingfor the first time at order BIMSRG(n) correspond to allthe topologies characterized by d max = n .4.5 Diagram evaluationOnce the complete set of diagrams contributing to C i j isgenerated, the expressions of the associated amplitudes , i.e.,of the anti-symmetrized matrix elements C i jk ... k i k i + ... k i + j , areobtained by applying a set of algebraic rules to each con-tributing diagram1. Label outgoing (incoming) external lines with quasi-particle indices k . . . k i ( k i + . . . k i + j ).2. Label the n l internal lines with di ff erent quasi-particleindices, e.g., p , q , r , ... .3. Attribute to vertices A kl and B mn their associated canoni-cal amplitudes.4. Sum over the n l internal line labels.5. Include a factor ( n l !) − for the equivalent internal lines.Equivalent internal lines begin and end at the same twovertices. By virtue of the fact that each diagram containsonly two vertices and that there is no self contraction, the n l internal lines are necessarily equivalent in the presentcontext.6. Multiply by a phase factor ( − n c , where n c is the numberof line crossings in the diagram (vertices are not consid-ered line crossings).7. Sum over all distinct permutations P of labels of inequiva-lent incoming (outgoing) external lines, including a parityfactor ( − σ ( P ) associated with the signature σ ( P ) of thepermutation. Incoming (outgoing) external lines are in-equivalent if they do not originate from the same vertex.These algebraic rules are fully consistent with the ones forBCC [32] since they originate from the application of thestandard Wick’s theorem under the same conditions .4.6 Elementary exampleLet us demonstrate the use of both the topological and alge-braic rules for a simple example. For this purpose, we con-sider the diagram built out of A (22) and B (31) and contributing In BCC, there exists one more rule associated with so-called equivalentcluster amplitudes. p qk k k k A (22) k k pq B (31) pqk k Fig. 5
Fully-labeled form of diagram C (22 , to C , i.e., diagram C (22 , A (22) and B (31) consistently with Eq. (41), n l = p , q ) of equivalent in-ternal lines and two pairs ( k , k ) and ( k , k ) of inequivalentexternal lines. As drawn, the amplitudes are given in canon-ical form and the diagram displays no line crossing. As aresult, the algebraic expression of the anti-symmetrized ma-trix element is given by C k k k k (22 , = P ( k k / k ) 12 (cid:88) pq A k k pq B pqk k . (42)In Eq. (42), the permutation operator P ( s / s ) = P ( s / s )generically performs an appropriate anti-symmetrization ofthe matrix element it acts on. The operator does so by per-muting all indices in set s with all indices in set s whileincluding a sign associated with the signature of the permuta-tion. In the above example, the operator is defined by P ( k k / k ) ≡ − P k k − P k k , (43)where, e.g., P k k commutes indices k and k . Because theelementary commutator involves two operators, the requiredpermutation operators necessarily involve only two sets ofindices; i.e. s and s . This is at variance with, e.g., BCCwhere a larger number of operators may appear such thatpermutation operators involving more than two sets of indicesare called for [32]. Still, the content of the two sets ofindices presently depend on the classes of terms entering A , B and C . As an example, the complete list of permutationoperators involved at the BIMSRG(2) level is provided inEq. (A.2).Because C k k k k (22 ,
31) is anti-symmetric with respectto the exchange of any pair of labels within the triplet ( k , k , k ),it is su ffi cient to compute it explicitly for, e.g., the subset k < k < k . In, e.g., BCC, contributions may involve more than two operatorswith external legs, e.g., double-excitation amplitudes constructed fromone T , one T operator and a two-body vertex of the grand-canonicalpotential, as the diagram labelled D6c in Ref. [32]. As such, permuta-tion operators between more than two sets of indices appear in suchformalisms. N A , N B ; N C ) = (2 ,
2; 2)approximation to the elementary commutator under the hy-pothesis that C is Hermitian / anti-Hermitian.The diagrams contributing to { C i j , i ≥ j } in this approx-imation are displayed in Fig. 6. They correspond to all dia-grams with topologies characterized by d max ≤
2. Allowingthe exchange [ A ↔ B ] and exploiting the (anti-)Hermiticityof C , these 28 diagrams actually account for a total of 82diagrams.The algebraic expressions of the BIMSRG(2) flow equa-tions are obtained by applying the algebraic rules to thediagrams displayed in Fig. 6 and by identifying A → η ( s ) , B → Ω ( s ) , (44) C → d Ω d s ( s ) , for the case of direct integration or A → M ( s ) , B → ad ( l − M ( η ) , (45) C → ad ( l ) M ( η ) , in the case of the Magnus expansion. The evaluated expres-sions are collected in Appendix A.4.8 BIMSRG(n)While the BIMSRG(2) truncation is an e ffi cient workhorse,one eventually wishes to push to BIMSRG(3) in order toobtain high-accuracy results. The numerical code describedbelow allows one to obtain the corresponding equations in amatter of seconds. The additional 82 diagrams arising in the( N A , N B ; N C ) = (3 ,
3; 3) approximation are too numerous andthe equations too lengthy to be displayed in this article; theinterested reader is invited to run the
ADG code accordingly.For orientation, Tab. 1 lists the number of diagrams aris-ing from the commutator in the ( N A , N B ; N C ) = ( n , n ; n ) ap-proximation up to n =
9, i.e. all diagrams with topologiescharacterized by d max = n up to n = C can go Performing the exchange [ A ↔ B ] and exploiting the Hermitian / anti-Hermitian character of C , these extra 82 diagrams actually account fora total of 264 diagrams.0 C = C (02 , + C (04 , − (cid:104) ↔ (cid:105) C = C (11 , + C (02 , + C (22 , + C (13 , − (cid:104) ↔ (cid:105) C = C (11 , + C (02 , + C (13 , + C (02 , + C (04 , + C (13 , − (cid:104) ↔ (cid:105) C = C (11 , + C (31 , + C (22 , − (cid:104) ↔ (cid:105) C = C (31 , + C (11 , + C (22 , + C (13 , + C (22 , + C (02 , − (cid:104) ↔ (cid:105) C = C (22 , + C (04 , + C (13 , + C (11 , + C (22 , + C (02 , + C (13 , − (cid:104) ↔ (cid:105) Fig. 6
Diagrams contributing to the (2 ,
2; 2) commutator under the hypothesis that C is Hermitian or anti-Hermitian, i.e., only the subset { C ij , i ≥ j } is explicitly computed. up to a factor of 2 in the ideal case where the number of C kk terms is very small with respect to other diagrams. Togetherwith using the symmetry in the truncation, this yields an ideal reduction by a factor of 4 in the number of distinct diagramsthat have to be considered explicitly. N A = N B N A = N B and Hermiticity 4 24 82 208 452 830 1436 2320 3558 Table 1
Number of new diagrams appearing at each order of the BIMSRG(n) truncation. The naive counting refers to computing explicitly alldiagrams contributing to both + AB and − BA . The growing number of diagrams can be reduced by using the symmetry in the truncation and theHermiticity of C . ADG code
Given the set of diagrammatic rules, all diagrams contributingto C = [ A , B ] in the ( N A , N B ; N C ) scheme are automaticallyproduced and evaluated via the ADG code. The values of N A and N B do not have to be the same whereas N C can bechosen to take any value between 0 and N A + N B −
1; i.e. thecode operates is not limited to the specific constraints of theBIMSRG(n) truncation scheme.The code generates expressions according to the follow-ing algorithms.5.1 Diagram generation1. Select either A or B to be the top vertex and associate B (resp. A ) with the bottom vertex.(a) Select a valid vertex degree d C in { , , . . . , N C } stip-ulating the total number of external lines of the outputoperator.i. Partition d C between incoming and outgoing ex-ternal lines, i.e. pick integers i and j to select C i j such that d C ij ≡ i + j = d C .A. Connect the external legs to A and B ver-tices in all possible ways. This must be donewhile keeping the number of these connec-tions strictly smaller than 2 N A and 2 N B , re-spectively, in order to allow for at least onecontraction between A and B . Additionally,the numbers of external legs connected tovertex A and to vertex B have to be eitherboth even or both odd, restricting the num-ber of configurations.B. Add n l internal lines connecting vertex A tovertex B in all possible ways. The possiblevalues of n l must be such that Eq. (41) isfulfilled with the degrees d A kl and d B mn ofthe involved vertices A kl and B mn being in { , , . . . , N A } and { , , . . . , N B } , respec-tively.ii. Go back to i. and exhaust all valid partitions.(b) Go back to (a) and exhaust all valid vertex degrees d C . 2. Go back to 1. and repeat the operation with permutedvertices.5.2 Diagram evaluation1. If the top vertex corresponds to the operator B , insert aminus sign as the diagram is associated to the − BA termof the commutator.2. If incoming (resp. outgoing) external legs are tied toboth A and B vertices, write a permutation operator P associated to both groups of lines. Permutations are to betaken between the group of lines entering (resp. leaving)vertex A and the group entering (resp. leaving) vertex B .3. Write a symmetry factor 1 / ( n l !) where n l is the numberof internal lines connecting the two vertices.4. For each vertex, write the matrix element of the corre-sponding operator with the quasi-particle labels of theattached out-going and incoming lines in canonical order.5. Add the factor ( − n c where n c denotes the number oflines crossing.6. Write a summation over the quasi-particle labels associ-ated to the internal lines.5.3 Output of the programThe typical output associated to a BIMSRG diagram in thecode is Diagram 56 ( − BA):C (40 , = − (cid:88) p p B k k p p A p p k k All BIMSRG diagrams generated by
ADG are made in such a waythat both their graphical representation and algebraic expression are incanonical form. Additionally, there is a direct correspondence betweena diagram expression and its drawing as obtained by the code.2
ADG program
ADG has been designed to work on any computer with aPython3 distribution, and successfully tested on commonLinux distributions as well as MacOS. In addition to thePython base, the setuptools and distutils packages must al-ready be installed, which is the case in most current distri-butions. Having pip installed eases the process but is nottechnically required. The
NumPy , NetworkX and
SciPy li-braries are automatically downloaded during the installationprocess. Additionally, one needs a L A TEX distribution installedwith the PDFL A TEX compiler for
ADG to produce a pdf filefrom the output, if so desired.6.1 Installation
The easiest way to install
ADG is to obtain it from the PythonPackage Index by entering the following command pip3 install adg Provided setuptools is already installed, pip takes care ofdownloading and installing
ADG as well as
NumPy , SciPy and
NetworkX . Once a new version of
ADG is released, one caninstall it by entering the command pip3 install --upgrade adg
Once the
ADG source files are downloaded from the GitHubrepository , one must enter the project folder and either run pip3 install . or python3 setup.py install With this method, pip also takes care of downloading andinstalling NumPy , NetworkX and
SciPy .6.2 Run the program
The most convenient way to use
ADG is to run it in batch modewith the appropriate flags. For example, to run the programand generate BIMSRG diagrams at order 3, i.e. using N A = N B = N C =
3, one can use https://pypi . org/project/adg/ https://github . com/adgproject/adg Depending on the system, it might be necessary either to use the"–user" flag to install it only for a specific user or to run the previouscommand with "sudo -H" to install it system-wide. adg -o 3 -t BIMSRG -d -c where the -o flag is for the order, -t for the type of theory, -d indicates that the diagrams must be drawn and -c that ADG must compile the L A TEX output. The program supportsless conventional truncation schemes as well. Two integersafter the -o flag will be interpreted as settings for N A = N B and N C respectively, and three integers as N A , N B and N C respectively. A complete list of the command-line optionscan be found in the program’s documentation (see Sec. 6.4)or by typing adg -h Currently,
ADG supports HF-MBPT by using -t MBPT ,straight BMBPT by using -t BMBPT , o ff -diagonal BMBPTby using -t PBMBPT and BIMSRG by using -t BIMSRG .Though the algorithms described in the previous sections canbe used regardless of the truncation order, ADG has been arbi-trarily restricted to order 10 or less to avoid major overloadsof the system due to rapidly growing number of expressionsor diagrams. Users are advised to first launch calculations atlow orders (2, 3 or 4 typically) to develop an intuition for thetime and memory needs for specific applications.
As an alternative to the batch mode,
ADG can be run on aterminal by entering the command adg -i
A set of questions must be answered using the keyboard toconfigure and launch the calculation. The interactive modethen proceeds identically to the batch mode.6.3 Steps of a program runHere, we briefly describe the main steps of a typical
ADG run: – Select options by using the command-line flags or key-board input by the user in an interactive session. – ADG creates a list of adjacency matrices for the appro-priate theory and order using
NumPy , and feeds them to
NetworkX , which creates
MultiDiGraph objects. – If needed, checks are performed on the list of graphs toremove topologically equivalent or ill-defined graphs. – The list of topologically unique graphs is used to produce
Diagram objects that store the graph as well as some ofits associated properties depending on the theory (HFstatus, excitation level, etc.). The expressions associatedto the graphs are eventually extracted. – The program prints on the terminal the number of dia-grams per category and writes the L A TEX output file, thedetails of which depend on the options selected by theuser, as well as a list of adjacency matrices associated to the diagrams. Other output files may be produced,depending on the theory and the user’s input. – If asked by the user, the program performs the PDFL A TEXcompilation. – Unnecessary temporary files are removed and the pro-grams exits.6.4 Documentation
Once the source files have been downloaded, a quick startguide is available in the
README.md file. Once
ADG is in-stalled, it is possible to read its manpages through man adg or a brief description of the program and its options through adg -h
A more detailed HTML documentation can be generateddirectly from the source files by going into the docs directoryand running make html
The documentation is then stored in docs/build/html ,with the main file being index.html . A list of other sup-ported documentation formats is available by running make help
The full HTML documentation is available online under https://github . com/adgproject/adg and support forthe program can be obtained by reporting bugs or other issueson the GitHub repository’s tracker at https://github . com/adgproject/adg . A novel single-reference expansion many-body method ap-propriate to the ab initio description of superfluid nuclei isformulated in terms of a particle-number-breaking Bogoli-ubov vacuum: The Bogoliubov in-medium similarity renor-malization group (BIMSRG) formalism. Such an extensionof the standard, i.e., symmetry-conserving, single-referenceIMSRG approach parallels those developed recently withinthe framework of coupled cluster, self-consistent Green’sfunction and many-body perturbation theories. Furthermore,BIMSRG complements the already existing multi-referencein-medium similarity renormalization group (MR-IMSRG)formalism that is also applicable to open-shell systems. The derivation of BIMSRG equations relies on the com-putation of the commutator between two operators of poten-tially high ranks. To perform such a computation e ffi ciently,a diagrammatic method is designed on the basis of a normal-ordered form of the operators with respect to the Bogoliubovvacuum. While low-rank truncations of the commutator caneasily be worked out manually on the basis of the diagram-matic rules, it becomes impractical and error-prone to workout higher orders, e.g. BIMSRG(3) equations, in this way.To overcome this limitation, we introduce the third version(v3.0.0) of the ADG code, which can automatically (1) gen-erate all valid BIMSRG(n) diagrams and (2) evaluate theiralgebraic expressions in a matter of seconds. This is achievedin such a way that equations can easily be retrieved for boththe flow equation and the Magnus expansion variants of BIM-SRG.Based on this work, the first objective is to perform full-fledged ab initio calculations of singly open-shell nuclei atthe BIMSRG(2) level by employing a spherical implemen-tation that also makes use of angular-momentum couplingtools [52]. Later on, an extension to BIMSRG(3) can beenvisioned. Another future direction is to extend the
ADG code [36, 44] to automatically derive equations at play withinMR-IMSRG, in particular to systematically and safely accessthose defining the MR-IMSRG(3) approximation.
Acknowledgements
The authors thank S. Bogner and A. Schwenk for usefuldiscussions. This publication is based on work supported inpart by the Deutsche Forschungsgemeinschaft (DFG, Ger-man Research Foundation) – Projektnummer 279384907 –SFB 1245, the Max Planck Society, by the BMBF ContractNo. 05P18RDFN1, and the National Science Foundationunder Grants No. PHY-1614130.
Data Availability Statement
The datasets generated during and / or analysed during thecurrent study are available in the GitHub repository, https://github . com/adgproject/adg . Appendix A: BIMSRG(2)
Appendix A.1: BIMSRG(2) flow equationFollowing the strategy laid out in Sec. 4.7, the BIMSRG(2)flow equations for the normal-ordered pieces of the grandpotential readd Ω d s = (cid:88) pq η pq Ω pq + (cid:88) pqrs η pqrs Ω pqrs − [ η ↔ Ω ] , (A.1a)d Ω k k d s = P ( k / k ) (cid:88) p η k p Ω k p + (cid:88) pq η k k pq Ω pq + (cid:88) pq η pq Ω k k pq + P ( k / k ) 13! (cid:88) pqr η k pqr Ω k pqr − [ η ↔ Ω ] , (A.1b)d Ω k k d s = (cid:88) p η k p Ω pk + (cid:88) p η k p Ω pk + (cid:88) pq η k k pq Ω pq + (cid:88) pq η pq Ω pqk k + (cid:88) pqr η k pqr Ω pqrk + (cid:88) pqr η k pqr Ω pqrk − [ η ↔ Ω ] , (A.1c)d Ω k k k k d s = P ( k k k / k ) (cid:88) p η k p Ω k k k p + P ( k / k k k ) (cid:88) p η k k k p Ω k p + P ( k k / k k ) 12 (cid:88) pq η k k pq Ω k k pq − [ η ↔ Ω ] , (A.1d)d Ω k k k k d s = P ( k k / k ) (cid:88) p η k p Ω k k pk + P ( k k / k ) (cid:88) p η k k k p Ω pk + (cid:88) p η pk Ω k k k p + (cid:88) p η k k k p Ω pk + P ( k k / k ) 12 (cid:88) pq η k k pq Ω pqk k + P ( k k / k ) 12 (cid:88) pq η k pqk Ω k k pq − [ η ↔ Ω ] , (A.1e)d Ω k k k k d s = (cid:88) pq η pqk k Ω k k pq + P ( k / k ) P ( k / k ) 12 (cid:88) pq η k k pq Ω pqk k + (cid:88) pq η k k pq Ω pqk k + P ( k / k ) (cid:88) p η k p Ω pk k k + P ( k / k ) (cid:88) p η k k pk Ω pk + P ( k / k ) (cid:88) p η k k k p Ω pk + P ( k / k ) (cid:88) p η pk Ω k k pk − [ η ↔ Ω ] , (A.1f)where the required permutation operators are defined through P ( k / k ) ≡ − P k k , (A.2a) P ( k / k k ) ≡ − P k k − P k k , (A.2b) P ( k / k k k ) ≡ − P k k − P k k − P k k , (A.2c) P ( k k / k k ) ≡ − P k k − P k k − P k k − P k k + P k k P k k . (A.2d)The shorthand-notation [ η ↔ Ω ] indicates that additionalcontributions are to be obtained by exchanging η and Ω inall previous terms. The specific (anti-)symmetry propertiesof the matrix elements of (d Ω i j / d s ) i j can be used to onlycompute them for a subset of quasi-particle indices.Note that even when a canonical HFB vacuum is em-ployed at the beginning of the flow, i.e., Ω k k (0) = Ω k k (0) =
0, the corresponding evolved matrix elements are non-zerodue to RG-corrections from contracted normal-ordered two-body components.Appendix A.2: BIMSRG(2) Magnus expansionThe Magnus expansion relies on the computation of the l -foldnested commutator ad ( l ) M ( η ) defined recursively in Eq. (15).At each step of the nesting, the BIMSRG(2) corresponds tosetting A → M ( s ) / B → ad ( l − M ( η ) / C → ad ( l ) M ( η ) . (A.3)in the ( N A , N B ; N C ) = (2 ,
2; 2) approximation to the elemen-tary commutator. The corresponding equations could be eas-ily obtained by adapting Eq. (A.1). This procedure is to berepeated iteratively until the maximum value of l is reached. Appendix A.3: Computational scalingThe computational scaling of the BIMSRG tensor network isdriven by the commutator involving the highest particle rank,e.g., in the case of BIMSRG(2)[ η [4] , Ω [4] ] [4] . (A.4)Equation (A.4) requires N evaluations, where N charac-terizes the size of the quasi-particle basis or, equivalently,the dimension of the one-body Hilbert space. In standardIMSRG(2), the scaling is N as well, but certain types of gen-erators allow further reductions. For the frequently employedWhite and imaginary-time generators, the computational ef-fort is driven by the evaluation of particle-particle ladderdiagrams [26]. It scales as n p n h where n p and n h denote thenumber of particle and hole states, respectively, similar toCoupled Cluster with singles and doubles (CCSD). Since n p (cid:29) n h , the quasi-particle formulation involves a signifi-cantly higher cost even though the scaling is bound to bepolynomial in N at any finite truncation order. As such, theevaluation of BIMSRG(2) is comparable to the solution ofthe BCC equations at the singles and doubles level (BCCSD)once appropriate intermediates are defined [32].Appendix A.4: Perturbation theory diagnosticsThe decoupling of the Bogoliubov reference state from itsquasi-particle excitations may be quantified during the flowin terms of various tensor norms of operator matrix elements.The repeated evaluation of low-order MBPT energy correc-tions was also successfully used as a diagnostic in previousIMSRG studies. Since BIMSRG(2) resums all diagrams upto third order, second- and third-order grand-potential correc-tions must eventually vanish at convergence. This propertycan thus be used as a diagnostic tool.Limiting the flowing grand potential to Ω [4] ( s ) terms, theRayleigh-Schrodinger partition is introduced as Ω ( s ) = Ω ( s ) + Ω ( s ) , (A.5)with Ω ( s ) ≡ Ω ( s ) + ¯ Ω ( s ) , (A.6a) Ω ( s ) ≡ Ω ( s ) + Ω ( s ) + ˘ Ω ( s ) + Ω ( s ) + Ω ( s ) + Ω ( s ) + Ω ( s ) + Ω ( s ) , (A.6b)and where the splitting Ω ( s ) ≡ ¯ Ω ( s ) + ˘ Ω ( s ) (A.7)is introduced such that the unperturbed part is defined as adiagonal one-body operator. Low-order corrections to the ground-state grand potential are evaluated via BMBPT [12,35]. For example, the second-order reads asE (2)0 ( s ) − λ A (2)0 ( s ) ≡ (cid:88) pq Ω pq ( s ) Ω pq ( s )¯ Ω pp ( s ) + ¯ Ω qq ( s ) (A.8) − (cid:88) pqrs Ω pqrs ( s ) Ω pqrs ( s )¯ Ω pp ( s ) + ¯ Ω qq ( s ) + ¯ Ω rr ( s ) + ¯ Ω ss ( s ) . Even starting from a canonical HFB vacuum, non-canonicaldiagrams need to be evaluated consistently during the flowsince Ω , pq ( s ) (cid:44) s >
0. Canonical and non-canonicalthird-order diagrams are additionally needed and are conve-niently obtained using
ADG [36].Appendix A.5: GeneratorsSo far, no specific form of the generator has been assumed. Inthe following, generators appropriate to BIMSRG(2) calcu-lations are derived. Given that the flowing grand potential islimited to Ω [4] ( s ) terms, it is su ffi cient to drive the couplingsto two and four quasi-particle excitations (cid:104) Φ | Ω ( s ) | Φ pq (cid:105) = Ω pq ( s ) (A.9a) (cid:104) Φ | Ω ( s ) | Φ pqrs (cid:105) = Ω pqrs ( s ) , (A.9b)to zero as s → ∞ . To do so, one must employ generators ofthe form η ( s ) ≡ η ( s ) + η ( s ) + η ( s ) + η ( s ) , (A.10)with matrix elements defined by η pq ( s ) ≡ G pq Ω pq ( s ) , (A.11a) η pq ( s ) ≡ G pq Ω pq ( s ) , (A.11b) η pqrs ( s ) ≡ G pqrs Ω pqrs ( s ) , (A.11c) η pqrs ( s ) ≡ G pqrs Ω pqrs ( s ) . (A.11d)In Eq. (A.11), the tensor G has been included to enforceanti-Hermiticity of η , i.e., G pq = − G qp , (A.12a) G pqrs = − G rspq . (A.12b)Specific choices for η relate to fixing G .In order to analyse the performance of the above ansatzfor the generator, a perturbative analysis [26] of the flow equa-tion can be performed on the basis of the so-called Epstein-Nesbet (EN) partitioning of the grand potential Ω ( s ) = ˜ Ω ( s ) + ˜ Ω ( s ) (A.13) The numbering of the orders used here corresponds to the usual con-vention among MBPT practitioners, at variance with what is used in the
ADG code [36] mentioned below.6 where˜ Ω ( s ) ≡ Ω ( s ) + (cid:88) p Ω pp ( s ) β † p β p + (cid:88) pq Ω pqpq ( s ) β † p β † q β q β p , (A.14a)˜ Ω ( s ) ≡ (cid:88) p (cid:44) q Ω pq ( s ) β † p β q + (cid:88) pqrs ( pq ) (cid:44) ( rs ) Ω pqrs ( s ) β † p β † q β s β r + Ω ( s ) + Ω ( s ) + Ω ( s ) + Ω ( s ) , (A.14b)such that diagonal matrix elements of Ω ( s ) have been in-cluded into the unperturbed grand potential O ( λ ). Countinginsertions of ˜ Ω ( s ) ( ˜ Ω ( s )) as order O ( λ ) ( O ( λ )), the pertur-bative expansion of Eqs. (A.1b) and (A.1d) yieldsdd s Ω pq = − ( Ω pp + Ω qq ) η pq − Ω pqpq η pq + O ( λ ) ≡ − ∆ pq η pq + O ( λ ) , (A.15a)dd s Ω pqrs = − ( Ω pp + Ω qq + Ω rr + Ω ss ) η pqrs − ( Ω pqpq + Ω prpr + Ω psps + Ω qrqr + Ω qsqs + Ω rsrs ) η pqrs + O ( λ ) ≡ − ∆ pqrs η pqrs + O ( λ ) , (A.15b)where ∆ pq ( s ) ≡ (cid:104) Φ pq | Ω ( s ) | Φ pq (cid:105) − (cid:104) Φ | Ω ( s ) | Φ (cid:105) = Ω pp ( s ) + Ω qq ( s ) + Ω pqpq ( s ) , (A.16a) ∆ pqrs ( s ) ≡ (cid:104) Φ pqrs | Ω ( s ) | Φ pqrs (cid:105) − (cid:104) Φ | Ω ( s ) | Φ (cid:105) = Ω pp ( s ) + Ω qq ( s ) + Ω rr ( s ) + Ω ss ( s ) + Ω pqpq ( s ) + Ω prpr ( s ) + Ω psps ( s ) + Ω qrqr ( s ) + Ω qsqs ( s ) + Ω rsrs ( s ) , (A.16b)define EN energy denominators that are fully symmetricunder the permutation of any pair of quasi-particle indices.Inserting Eq. (A.11) into Eq. (A.15), one obtainsdd s Ω pq = − ∆ pq G pq Ω pq + O ( λ ) , (A.17a)dd s Ω pqrs = − ∆ pqrs G pqrs Ω pqrs + O ( λ ) , (A.17b)which can be formally integrated while ignoring the s depen-dence of the EN energy denominators as Ω pq ( s ) = Ω pq (0) e − s ∆ pq G pq , (A.18a) Ω pqrs ( s ) = Ω pqrs (0) e − s ∆ pqrs G pqrs . (A.18b)Under the assumption that ∆ pq and ∆ pqrs are positive , Ω pq ( s )and Ω pqrs ( s ) are exponentially suppressed such that the de-coupling invoked in Eq. (11) is indeed realized for the BIM-SRG(2) approximation when s → ∞ . It is to be noted that To ensure the suppression in Eq. (A.18), one can slightly modifyEq. (A.11) according to η pq ≡ sgn( ∆ pq ) G pq Ω pq , (A.19a) the rate at which the matrix elements Ω pq ( s ) and Ω pqrs ( s )are suppressed increases with the size of the associated ENdenominator | ∆ pq | and | ∆ pqrs | , as is customary for an RGevolution. References
1. P. Navratil, S. Quaglioni, I. Stetcu, and B. R. Barrett, J.Phys. G , 083101 (2009).2. R. Roth, J. Langhammer, A. Calci, S. Binder, andP. Navrátil, Phys. Rev. Lett. , 072501 (2011).3. B. R. Barrett, P. Navrátil, and J. P. Vary, Prog. Part. Nucl.Phys. , 131 (2013).4. A. Gezerlis, I. Tews, E. Epelbaum, S. Gandolfi,K. Hebeler, A. Nogga, and A. Schwenk, Phys. Rev. Lett. , 032501 (2013).5. J. Carlson, S. Gandolfi, F. Pederiva, S. C. Pieper, R. Schi-avilla, K. E. Schmidt, and R. B. Wiringa, Rev. Mod. Phys. , 1067 (2015).6. J. E. Lynn, I. Tews, J. Carlson, S. Gandolfi, A. Gezerlis,K. E. Schmidt, and A. Schwenk, Phys. Rev. C , 054007(2017).7. R. Roth and J. Langhammer, Phys. Lett. B , 272(2010).8. J. Langhammer, R. Roth, and C. Stumpf, Phys. Rev. C , 054315 (2012).9. J. D. Holt, J. Menéndez, J. Simonis, and A. Schwenk,Phys. Rev. C , 024312 (2014).10. A. Tichai, J. Langhammer, S. Binder, and R. Roth, Phys.Lett. B , 283 (2016).11. A. Tichai, E. Gebrerufael, K. Vobig, and R. Roth, Phys.Lett. B , 448 (2018).12. A. Tichai, R. Roth, and T. Duguet, Front. Phys. , 164(2020).13. W. H. Dickho ff and C. Barbieri, Prog. Part. Nucl. Phys. , 377 (2004).14. C. Barbieri and M. Hjorth-Jensen, Phys. Rev. C ,064313 (2009).15. A. Cipollone, C. Barbieri, and P. Navrátil, Phys. Rev.Lett. , 062501 (2013).16. F. Raimondi and C. Barbieri, Phys. Rev. C , 054327(2019).17. V. Somà, Front. Phys. , 340 (2020).18. G. Hagen, T. Papenbrock, D. J. Dean, and M. Hjorth-Jensen, Phys. Rev. C , 34330 (2010).19. G. Hagen, T. Papenbrock, M. Hjorth-Jensen, and D. J.Dean, Rep. Prog. Phys. , 096302 (2014). η pq ≡ sgn( ∆ pq ) G pq Ω pq , (A.19b) η pqrs ≡ sgn( ∆ pqrs ) G pqrs Ω pqrs , (A.19c) η pqrs ≡ sgn( ∆ pqrs ) G pqrs Ω pqrs . (A.19d)7
20. S. Binder, J. Langhammer, A. Calci, and R. Roth, Phys.Lett. B , 119 (2014).21. S. J. Novario, G. Hagen, G. R. Jansen, and T. Papenbrock,Phys. Rev. C , 051303 (2020).22. K. Tsukiyama, S. K. Bogner, and A. Schwenk, Phys.Rev. Lett. , 222502 (2011).23. K. Tsukiyama, S. K. Bogner, and A. Schwenk, Phys.Rev. C , 061304(R) (2012).24. H. Hergert, S. K. Bogner, S. Binder, A. Calci, J. Lang-hammer, R. Roth, and A. Schwenk, Phys. Rev. C ,034307 (2013).25. S. K. Bogner, H. Hergert, J. D. Holt, A. Schwenk,S. Binder, A. Calci, J. Langhammer, and R. Roth, Phys.Rev. Lett. , 142501 (2014).26. H. Hergert, S. K. Bogner, T. D. Morris, A. Schwenk, andK. Tsukiyama, Phys. Rept. , 165 (2016).27. S. R. Stroberg, A. Calci, H. Hergert, J. D. Holt, S. K.Bogner, R. Roth, and A. Schwenk, Phys. Rev. Lett. ,032502 (2017).28. S. R. Stroberg, H. Hergert, S. K. Bogner, and J. D. Holt,Annu. Rev. Nucl. Part. Sci. , 307 (2019).29. V. Somà, T. Duguet, and C. Barbieri, Phys. Rev. C ,064317 (2011).30. V. Somà, C. Barbieri, and T. Duguet, Phys. Rev. C ,011303 (2013).31. V. Somà, A. Cipollone, C. Barbieri, P. Navrátil, andT. Duguet, Phys. Rev. C , 061301(R) (2014).32. A. Signoracci, T. Duguet, G. Hagen, and G. Jansen, Phys.Rev. C , 064320 (2015).33. T. M. Henderson, J. Dukelsky, G. E. Scuseria, A. Signo-racci, and T. Duguet, Phys. Rev. C , 054305 (2014).34. T. Duguet and A. Signoracci, J. Phys. G , 015103(2017), [Erratum: J.Phys.G 44, 049601 (2017)].35. A. Tichai, P. Arthuis, T. Duguet, H. Hergert, V. Somá,and R. Roth, Phys. Lett. B , 195 (2018).36. P. Arthuis, T. Duguet, A. Tichai, R. D. Lasseri, and J. P.Ebran, Comput. Phys. Commun. , 202 (2019).37. P. Demol, M. Frosini, A. Tichai, V. Somà, and T. Duguet,Annals Phys. , 168358 (2021).38. T. D. Morris, J. Simonis, S. R. Stroberg, C. Stumpf,G. Hagen, J. D. Holt, G. R. Jansen, T. Papenbrock,R. Roth, and A. Schwenk, Phys. Rev. Lett. , 152503(2018).39. S. R. Stroberg, J. D. Holt, A. Schwenk, and J. Simonis,Phys. Rev. Lett. , 022501 (2021).40. J. M. Yao, B. Bally, J. Engel, R. Wirth, T. R. Rodríguez,and H. Hergert, Phys. Rev. Lett. , 232501 (2020).41. H. Hergert, S. Binder, A. Calci, J. Langhammer, andR. Roth, Phys. Rev. Lett. , 242501 (2013).42. H. Hergert, S. K. Bogner, T. D. Morris, S. Binder,A. Calci, J. Langhammer, and R. Roth, Phys. Rev. C , 041302(R) (2014). 43. E. Gebrerufael, K. Vobig, H. Hergert, and R. Roth, Phys.Rev. Lett. , 152503 (2017).44. P. Arthuis, A. Tichai, J. Ripoche, and T. Duguet, Comput.Phys. Commun. , 107677 (2021).45. P. Ring and P. Schuck, The Nuclear Many-Body Problem (Springer Verlag, New York, 1980).46. T. D. Morris, N. M. Parzuchowski, and S. K. Bogner,Phys. Rev. C , 034331 (2015).47. W. Magnus, Commun. Pure Appl. Math. , 649 (1954).48. S. Blanes, F. Casas, J. Oteo, and J. Ros, Phys. Rep. ,151 (2009).49. G. C. Wick, Phys. Rev. , 268 (1950).50. J. Ripoche, A. Tichai, and T. Duguet, Eur. Phys. J. A ,02 (2020).51. I. Shavitt and R. J. Bartlett, Many-Body Methods inChemistry and Physics: MBPT and Coupled-ClusterTheory , Cambridge Molecular Science (Cambridge Uni-versity Press, 2009).52. A. Tichai, R. Wirth, J. Ripoche, and T. Duguet, Eur. Phys.J. A56