THE SIMILARITY RENORMALIZATION GROUP AND FACTORIZATION TECHNIQUES By Boyao Zhu A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Physics – Doctor of Philosophy Computational Mathematics, Science & Engineering – Dual Major 2023 ABSTRACT In this work, we explore the use of factorization techniques as a means to identify low- rank structures in nuclear interactions and many-body methods, which we aim to leverage in order to reduce computational cost and memory requirements of modern nuclear many-body calculations. Using the Singular Value Decomposition (SVD), we show that we can construct accurate low-rank models of modern nucleon-nucleon interactions, and we develop a new variant of the Similarity Renormalization Group (SRG), that allows us to evolve nuclear interactions directly in the SVD-factorized form. Next, We extend these developments to three-nucleon interactions, using randomized SVD (R-SVD) algorithms to mitigate the impact of the much greater basis dimensions. While we find evidence of low-rank structures in the three-nucleon interaction, we also identify obstacles to a persistent rank reduction that we trace back to SRG-induced interactions. Finally, we explore applications of the SVD for factorization and rank reduction within the context of the In-Medium SRG (IMSRG). We develop rank-reduced IMSRG flow equations, which, while promising, still pose computational challenges stemming from the treatment of the so-called particle-hole terms. As a first step towards tackling this problem, we perform a detailed breakdown of the scaling and importance of all contributions in the IMSRG flow in applications for a schematic model as well as infinite neutron matter based on realistic interactions from chiral Effective Field Theory. ACKNOWLEDGMENTS I extend my deepest gratitude to all those who have provided unwavering support and encouragement throughout the completion of my dissertation in the field of computational physics. Their invaluable assistance has been instrumental in turning this academic endeavor into a tangible reality. First and foremost, I express my profound appreciation to my esteemed supervisor, Prof. Heiko Hergert, whose expert guidance and profound insights have been invaluable throughout this research journey. His unwavering encouragement and belief in my abilities have been a constant source of motivation and inspiration. I am also immensely thankful to the esteemed committee members, Prof. Scott K Bogner, Prof. Sean N. Liddick, Prof. Mark Iwen, Prof. Matthew Hirn, and Prof. Longxiu Huang from the Department of Physics and Center of Computational Mathematics, Science, and Engineering at Michigan State University. To Prof. Huang, I warmly welcome your in- volvement and eagerly anticipate benefiting from your expertise in the final stages of my work. Gratitude is extended to my colleagues Jacob Davison, Roland Wirth, and Jiangming Yao for their engaging discussions and insightful ideas. Their diverse perspectives have enriched my understanding of the subject matter and made this academic journey all the more rewarding. I would like to express my deep appreciation and send my best wishes to my cohort for their constant spiritual support and companionship during this journey. To each member, including Hao Lin, Xingze Mao, Tong Li, Mengzhi Chen, Avik Sakar, Zhite Yu, Omokuyani Udiani, Xiaoyi Sun, Zhouyou Fan, Rui Zhang, Kang Yu, Shuyue Xue, Ruixin Zhang, Lijiang Xu, Siyuan Ma, Min Feng, and anyone I met, your unwavering presence and encouragement iii have meant the world to me, and I am truly grateful for your invaluable support. I would be remiss not to express my heartfelt thanks to my family for their unwavering support and belief in my abilities. Their love and encouragement have been a constant source of strength throughout the ups and downs of this scholarly pursuit. Finally, my appreciation goes to all the participants who contributed to this research, making it possible to achieve profound and meaningful outcomes. With profound gratitude to all who have played an integral role in my academic journey, I am sincerely thankful for their invaluable contributions and encouragement. iv TABLE OF CONTENTS Part I Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 1 Introduction to ab initio Nuclear Many-Body Theory . . . . . . . . . . . 1 2 Chapter 2 The Similarity Renormalization Group . . . . . . . . . . . . . . . . . . . . 35 Chapter 3 The In-Medium SRG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 Chapter 4 Factorization Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Part II Factorization in Few-Body Systems . . . . . . . . . . . . . . . . . . . . . . . . . 74 Chapter 5 SVD-SRG evolution in Two-Nucleon Interactions . . . . . . . . . . . . . . 75 Chapter 6 SVD-SRG for Three-Nucleon Interactions . . . . . . . . . . . . . . . . . . 112 Part III Factorization in Many-Body Systems . . . . . . . . . . . . . . . . . . . . . . . . 129 Chapter 7 SVD and IMSRG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 Chapter 8 SVD-IMSRG for Schematic Models . . . . . . . . . . . . . . . . . . . . . . 141 Chapter 9 SVD-IMSRG for Infinite Matter . . . . . . . . . . . . . . . . . . . . . . . . 156 Chapter 10 Further Developments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 Part IV Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175 Chapter 11 Conclusions and Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178 v Part I Preliminaries 1 Chapter 1 Introduction to ab initio Nuclear Many-Body Theory The aim of theoretical nuclear physics is to accurately describe the structure and dy- namics of atomic nuclei. Over the past decade, new developments in many–body theory, in particular new insights into the nuclear interactions and rapid growth of many–body tech- nology, has greatly expanded the capabilities of ab initio1 [21, 40] nuclear structure and reaction theory in low–energy nuclear physics. 1.1 The nuclear many-body problem Nuclear structure and dynamics are governed by the time-independent and time-dependent many-body Schr¨odinger equations, and H ∣ΨA k ⟩ = Ek ∣ΨA k ⟩ i̵h d dt ∣Ψ(t)⟩ = H ∣Ψ(t)⟩ , (1.1) (1.2) 1This means ”from first principles”, i.e., starting from interactions between constituent nucleons and deriving all other observables from there. 2 respectively. The main task of nuclear many-body theory is to provide precise and accurate solutions to these equations. Here, ∣ΨA k ⟩ denotes the k-th eigenstate of the system with eigen- value Ek, and ∣Ψ(t)⟩ is a general time-dependent state describing the target and projectile(s) in a nuclear reaction. The nuclear Hamiltonian H arises from the modeling of the strong interaction at low energy using nucleonic degrees of freedom H = T + V [2] + V [3] + ⋯ (1.3) T denotes the (intrinsic) kinetic energy, and V [2] and V [3] are two-nucleon (NN) and three-nucleon (3N) interactions, systematically constructed from chiral effective field the- ory (χEFT), see Section 1.2, where χEFT is a systematic framework for describing the strong interaction at low momentum, using pions and nucleons as degrees of freedom, and incorporating the chiral symmetry of Quantum Chromodynamics . Before then, many phe- nomenological potentials, such as Argonne and CD-Bonn 2, describe nucleus from the bottom up by solving the Schr¨odinger equation. Nuclear many-body physics, when trying to solve Schr¨odinger equation, faces two major challenges. The first is that the nuclear forces between nuclei cannot be directly derived from QCD, the underlying theory of strong interaction. The second major challenge of nuclear theory lies in the difficulty of obtaining accurate results in an efficient way, due to the fact that the many-body Hilbert space grows in a combinatoric way as the number of particles increases. Thus, an exact diagonalization of the nuclear Hamiltonian, referred to as a no- 2they are rooted in pion or boson-exchange pictures to some degree, and give a very precise and accurate fit to NN scattering data, but lack the more systematic underpinning of chiral EFT, especially when it comes to consistent 3N (and higher) forces, as well as other operators like transitions and currents. 3 core Full Configuration Interaction (FCI) calculation, scales exponentially with the number of particles A, and quickly becomes unfeasible as A increases. 1.2 Interactions from Chiral Effective Field Theory Efforts to construct a fundamental theory of the interaction between nucleons date back as far as eight decades. With the discovery of the neutron by Chadwick [17] in 1932, it became evident that the atomic nuclei are built up from protons and neutrons. The concept of strong nuclear interactions was first introduced by Heisenberg [36], in order to explain why electromagnetic forces cannot bind the nucleus together. In 1935, the meson exchange theory by Hideki Yukawa [112] suggested that the nucleons create the force by exchanging particles between each other, in analogy to the theory of the electromagnetic interaction where the exchange of the particle, known as photon, is the cause of the force. By the 1960s, the systematic measurements of nucleon-nucleon (N N ) observables and phase shift [91] analysis had revealed the existence of a repulsive core, an attractive intermediate- ranged force, as well as strong tensor and spin-orbit forces [67], see Figure 1.1. Yukawa’s meson exchange, which included these features was an excellent phenomenology [66] for describing nuclear interactions. However, with the rise of Quantum Chromodynamics (QCD), we desire to understand nuclear interactions on a more fundamental level. Within the Standard Model of particle physics, QCD provides the theoretical for describ- ing the strong interaction, while nuclear interactions are considered as residuals of strong forces. According to the theory, the fundamental degrees of freedom are quarks and gluons, which carry an abstract color charge. Quarks interact via gluon exchange, but unlike other exchange particles, gluons also can also interact with themselves, because QCD is a is a 4 3S1 3S1 − 3D1 Figure 1.1: The Argonne V18 [110] nuclear potential in coordinate (left) and momentum (right) representations [41]. In the left panel, the radial dependencies of the central (VC), tensor (VT ) and spin-orbit interactions (VLS) in the S = 1, T = 0 spin- and isospin-channel are shown, while the right panel shows its momentum space matrix elements in the deuteron partial waves. non-Abelian gauge field theory. A consequence of this is that the QCD coupling constant becomes weak at high energies, a property that is referred to asymptotic freedom. Thus, QCD processes at high energy and momentum transfer can be described with perturbative methods. At low energies, where nuclear physics resides, the QCD coupling is large and the theory is inherently non-perturbative. The main approach for performing non-perturbative calculations in this regime is Lattice QCD (LQCD), which simulates QCD processes on a discrete space-time lattice. However, the range of applications of LQCD is limited due to the enormous computational cost, as well as algorithmic challenges, in particular related to extracting meaningful results from a noisy signal. Therefore, QCD calculations of finite nuclei are not feasible, and even a direct derivation 5 ���������������-�����������������[��]�(�)[���]������� of nuclear forces from QCD is extremely challenging. In the 1980s, Steven Weinberg applied the concepts of effective field theory (EFT) to de- velop a systematic and model-independent treatment of low-energy QCD that also accounts for the spontaneously broken chiral symmetry (see [102–104] for more detailed discussion), with pions and nucleons as effective degrees of freedom. This chiral EFT (χEFT) is estab- lished through a low-momentum expansion and is valid only at low energies, suppressing high-momentum contributions by introducing regularizations. According to Weinberg, one can construct effective Lagrangians that consist of all interactions that are consistent with the symmetries of low-energy QCD and organized by an expansion in powers of (Q/Λ). Here, Q is a typical momentum of the interacting system, and Λ is the so-called breakdown scale of the theory. Applying this Lagrangian to a certain process, an infinite number of Feynman diagrams appear, but they can be organized according to the power counting in (Q/Λ), as shown in Figure 1.2. These interactions consist of (multi-)pion exchanges between nucleons, denoted by dashed lines, as well as nuclear contact interactions. While the chiral interactions [15] constructed in this framework are nowadays the de facto standard input to ab initio many-body calculations, open challenges remain. Because the typical cutoffs Λ for χEFT interactions are lower than those of phenomeno- logical interactions like AV18 and CD-Bonn [65], they tend to be better behaved in ab initio calculations and results converge more rapidly. However, such calculations are still compu- tationally heavy. Merely using an even lower cutoff to obtain even softer χEFT potentials is problematic, because it will eventually affect their accuracy [51]. To solve this problem, we perform renormalization group (RG) evolutions of the nuclear interaction, which “soften” them while preserving relevant physical information. RG methods are natural colleagues to any EFT, since they allow us to dial the EFTs resolution scales, and are (in principle) able 6 Figure 1.2: Hirarchy of nuclear forces at increasing orders up to N4LO in chiral expansion arranged by the particle rank of the interaction in the Weinberg scheme [23]. Solid and dashed lines refer to nucleons and pions, respectively. Solid dots, filled circles, filled squares, filled diamonds, and open squares refer to vertices from the Lagrangian of dimension ∆ = 0, 1, 2, 3 and 4, respectively. 7 to event connect theories with different degrees of freedom [60, 61]. 1.3 Elements of Many-Body Theory Let us now review some of the essential ideas of many-body theory that will become relevant in later chapters. We will keep the discussion light on mathematical details and refer to Refs. [25, 82, 89] for details. 1.3.1 Many-Particle Spaces and Basic Concepts In nuclear physics, we are solely dealing with systems of fermions, protons and neutrons, which are usually regarded as identical particles with a different isospin quantum number when dealing with the strong interaction. In coordinate representation, an abstract one- particle state ∣i⟩ has the following wave function, φi(x) ≡ ⟨x ∣ i⟩ (1.4) where x ≡ (r, s, t) denotes the spatial (r) and internal coordinates (spin s, isospin t). The overlap of two single-particle states is easily obtained using the identity resolution 1 = ∫ dx ∣x⟩ ⟨x∣ (1.5) where 1 designates the identity operator in the product of the coordinate and internal space(s) and the integral symbol indicates the integration over the 3 dimensional Euclidean 8 space and the summation over discrete internal variable(s). Thus ⟨i ∣ j⟩ = ∫ dx ⟨i ∣ x⟩ ⟨x ∣ j⟩ = ∫ dxφ∗ i (x)φj(x) where ⟨i ∣ x⟩ and ⟨x ∣ i⟩ are complex conjugates. We typically work in single-particle bases that are orthonormalized, so that ⟨i ∣ j⟩ = δij . (1.6) Any single particle state ψ(x) ≡ ⟨x ∣ ψ⟩ can be expressed as a linear combination of the vectors in a complete basis {φi}i∈N or ψ(x) = ∑ i ciφi(x) ∣ψ⟩ = ∑ i ci ∣i⟩ (1.7a) (1.7b) where the summation generally extends over the full basis. The expansion coefficients ci are given as ci = ⟨i ∣ ψ⟩ = ∫ dx φ∗ i (x)ψ(x) (1.8) where the resolution of identity has been used again, Eq. (1.5). A complete orthonormal basis {∣i⟩} or {φi(x)} could thus be regarded as a basis for a one-particle Hilbert space V1.3 Consider next a system of N identical fermions, e.g., a system of N nucleon. Its state 3Technically, the states ∣i⟩ and {φi(x)} are different vector spaces, with the former being abstract while the latter being a product of square-integrable functions with a spin-isospin spinor space. However, these spaces are merely related by a natural isomorphism, as indicated in Eq. (1.4) and we thus use the identical symbol V1 for both. The same convention will be used for corresponding many-particle spaces. 9 space FN can be interpreted as an antisymmetric outer product of N single-particle Hilbert spaces FN = AVN (1.9) where A denotes the so-called antisymmetrizer, i.e., the projector onto the totally antisym- metric subspace of VN , and VN = V1 ⊗ V1 ⊗ ⋅ ⋅ ⋅ ⊗ V1 ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ N-times = V ⊗N 1 The product space VN consists of many-body states where none of the particles are inter- acting with each other, i.e., simple product states of the form φi1(x1) ⊗ φi2(x2) ⊗ ⋯ ⊗ φiN (xN ) (1.10a) or ∣i1⟩ ⊗ ∣i2⟩ ⊗ ⋯ ⊗ ∣iN ⟩ . (1.10b) For compactness, this is usually written as φi1(x1)φi2(x2)⋯φiN (xN ) (1.11a) or ∣i1⟩ ∣i2⟩ ⋯ ∣iN ⟩ ≡ ∣i1i2⋯iN ⟩ (1.11b) 10 A general normalized basis vector in FN can be constructed from the single-particle basis as ∣ΦK⟩ = √ N !A ∣i1i2⋯iN ⟩ (1.12) where K designates the index set K ≡ {i1, i2, ⋯, iN } . (1.13) The antisymmetrizer A has the explicit form A = (N !)−1 ∑ P ∈SN (−1)π(P )P (1.14) where the sum extends over N ! permutations P that constitute the symmetric group SN , π(P ) is the signature of P — i.e., the number of pairwise permutations it consists of — and (−1)π(P ) is referred to as the parity of P . Antisymmetrized product states of the form (1.12) are also referred to as Slater determinants as well as basis configurations in the following. The configurations ∣ΦK⟩ form a basis of the N -particle space FN so that a general state of any many-fermion system can be expanded as follows: where ∣Ψ⟩ = ∑ K cK ∣ΦK⟩ cK = ⟨ΦK ∣ Ψ⟩ . (1.15) (1.16) The single-particle states appearing in K are said to be occupied in this configuration while those not present are said to be unoccupied. 11 1.3.2 Occupation Number Representation To uniquely specify any given basis configuration, we need to specify a fixed ordering of the single-particle basis states, e.g., such that indices are in ascending order, and we need to know which single-particle states are occupied and which are not in ∣ΦK⟩. We can then introduce the occupation number representation, where a given ∣ΦK⟩, can be uniquely defined as ∣ΦK⟩ = ∣{i1, i2, ⋯, iN }⟩ ≡ ∣n1n2⋯nN ⋯⟩ (1.17) where for fermionic systems the occupation number nk for any configuration can only take on two values, 0 and 1. Clearly, nk = 1 if ∣ik⟩ is occupied in ∣ΦK⟩ for ik ∈ K. For example, a five-particle state ∣{i1, i2, ⋯, i5}⟩ with i1 = 0, i2 = 1, i3 = 4, i4 = 7 and i5 = 9, has the occupation number representation ∣Φ01479⟩ = ∣{0, 1, 4, 7, 9}⟩ = ∣110010010100⋯0⋯⟩ . (1.18) As we can see, the index-based notation for fermionic systems is usually more concise than the occupation number notation. The latter, however, reveals that basis configurations can be represented as bit strings. This is very convenient for numerical implementation, in particular for implementing the action of particle creation and annihilation operators, which will be discussed in the next subsection. 1.3.3 Creation and Annihilation Operators We are now in the proper position to introduce the time-independent operators ˆai and ˆa† i , which annihilate and create a particle in the single-particle state ∣i⟩, respectively. We 12 define the action of a fermionic annihilation operator as follows: ˆai ∣{ij2j3⋯jN }⟩ = ∣{j2, j3⋯jN }⟩ if i ∉ {j2, j3, ⋯, jN } , ˆai ∣{j1j2j3⋯qN }⟩ = 0 if i ∉ {j1, j2, j3, ⋯, jN } . (1.19) The annihilation operator annihilates a particle in the state i assuming there is a particle present, i.e., the state i is occupied (and distinct from all other occupied states). If there is no particle in that state, the annihilation operator makes the whole N -particle state vanish. In the first equation of Eq. (1.19) we have made an assumption that the state i occurs at the first position in the configuration for simplicity. If this is not the case, we can use the antisymmetrization of the fermionic state to obtain ∣{i1i2i3i4⋯ iN }⟩ = (−1)π(P ) ∣{j1j2j3j4⋯ jN }⟩ (1.20) where (−1)π(P ) is the parity of the permutation P that reorders the set {i1, i2, . . . , iN } into the set {j1, j2, . . . , jN } . For example ˆai ∣{ijkl⋯}⟩ = ∣{jkl⋯}⟩ ˆai ∣{jikl⋯}⟩ = ˆai(−1) ∣{ijkl⋯}⟩ = − ∣{jkl⋯}⟩ ˆai ∣{jkil⋯}⟩ = ˆai(−1)2 ∣{ijkl⋯}⟩ = ∣{jkl⋯}⟩ . With a fermionic state in the occupation number representation one can write Eq. (1.20) in the following compact form ˆai ∣n1n2⋯ni⋯⟩ = (−1)s(i)ni ∣n1n2⋯ ni − 1 ⋯⟩ (1.21) 13 where s(i) is the number of occupied states that are preceding the state i, which is also equivalent to the number of pairwise permutations that are needed to move the state i from its position to the first position: s(i) = ∑ k εF , and occupied or hole states with ε ≤ εF . With these conventions, we can write the Fermi vacuum as ∣Φ0⟩ = ( ∏ εi≤εF ˆa† i ) ∣∅⟩ (1.49) where the index i designates different single-particle states up to the Fermi level. Let us now modify our notation a bit and reserve ranges of indices for particular kinds of single-particle states4. We use letters ijkl⋯ for states below the Fermi level (holes) and 4The convention used in the following originates in quantum chemistry, see [89], for example. 19 i, j, k, . . . a, b, c, . . . p, q, r, . . . hole state particle state generic state Table 1.1: Single-Particle Labeling Convention. abcd⋯ for states above the Fermi level (particles), while general single-particle state are indicated by letters pqrs⋯ (cf. Table 1.1). As suggested above, we can now write a given one–particle–one–hole (1p1h) state as ∣Φa i ⟩ = ˆa† aˆai ∣Φ0⟩ , a (2p2h) state as ∣Φab ij ⟩ = ˆa† aˆa† bˆajˆai ∣Φ0⟩ and a general (ApAh) excited state as ∣Φabc⋯ ijk⋯⟩ = ˆa† aˆa† bˆa† c⋯ˆakˆajˆai ∣Φ0⟩ . We also see that ˆa† i ∣Φ0⟩ = 0 , ˆaa ∣Φ0⟩ = 0 , (1.50) (1.51) (1.52) because a creation operator can not create a particle under the Fermi level since the state was occupied in ∣Φ0⟩, and an annihilator cannot destroy a particle above the Fermi level since there are no particles present. 20 1.3.7 Normal Ordering and Wick’s Theorem Next, we introduce the idea of normal ordering the Hamiltonian with respect to the vacuum state ∣∅⟩. For two arbitrary creation or annihilation operators, we define XY ≡ {XY } + XY , (1.53) where the first term on the right {⋯}5 designates the normal-ordered product operator, and the brace indicates the contraction of X and Y . The normal-ordered operators are defined as the product of the same operators whose ordering is rearranged in such a way that all the creators are to the left of all the annihilators, keeping only the signs from anticommuting the operators, but not the Kronecker delta from Eq. (1.38), which is taken care of by the contraction. For example. we have {a† paqa† r} = {−a† pa† raq} = a† ra† paq . It is obvious that the expectation values of arbitrary normal-ordered operators with respect to the true vacuum ∣∅⟩ must vanish ⟨∅∣ {a† p⋯aq} ∣∅⟩ = 0 , (1.54) because the rightmost operator will always be an annihilator that kills the vacuum. For the contractions, we have four elementary combinations: 1. X = ap, Y = aq: 5This is to be distinguished from the anti-commutator {X, Y }. apaq = apaq − {apaq} = apaq − apaq = 0 , (1.55) 21 2. X = a† p, Y = a† q: 3. X = a† p, Y = aq: 4. X = ap, Y = a† q: pa† a† q = a† pa† q − {a† pa† q} = a† pa† q − a† pa† q = 0 , paq = a† a† paq − {a† paq} = a† paq − a† paq = 0 , apa† q = apa† q − {apa† q} = apa† q + a† qap = {ap, a† q} = δpq . (1.56) (1.57) (1.58) In the first three cases, the operator string is already in normal order, so it makes sense that the contractions vanish. In the fourth case, ⟨∅∣ {apa† q} ∣∅⟩ = 0 implies that the contraction is equal to the vacuum expectation value of the operator, ⟨∅∣ apa† q ∣∅⟩ = apa† q ⟨∅ ∣ ∅⟩ , (1.59) which vanishes unless (reading from right to left) we create and then annihilate a fermion in the same state p. The generalization of these considerations to longer strings of creators and annihilators is provided by Wick’s theorem [89, 108], which provides a convenient algorithm for expanding an operator string in products of normal-ordered operators and contractions as follows: X1⋯XN ≡ {X1⋯XN } + X 1X 2 {X3⋯XN } − X 1X 3 {X2X4⋯XN } + all other single contractions 22 + (X 1X 2X 3X 4 − X 1X 3X 2Y 4) {X5⋯XN } + all other double contractions + all triple contractions ⋯ + all full contractions . (1.60) Here we have assumed that the string consists of an even number of operators, so that it can be fully contracted. For example, pa† a† qasar = {a† pa† qasar} (1.61) because the operator is already in normal ordered form and all contractions vanish in this case. Another example is apaqa† sa† r = {apaqa† sa† r} − apa† r {a† saq} + apa† s {a† raq} + aqa† r {a† sap} − aqa† s {a† rap} + apa† raqa† s − apa† saqa† r = − δpr {a† saq} + δps {a† raq} + δqr {a† sap} − δqs {a† rap} + δprδqs − δpsδqr . (1.62) A very useful corollary of Wick’s theorem exists for products of normal-ordered operators (see [41, 108]): {X1⋯XN } {Y1⋯YM } ≡ {X1⋯XN Y1⋯YM } 23 + X N Y 1 {X1⋯XN −1Y2⋯YM } − X N Y 2 {X1⋯XN −1Y1Y3⋯YM } + all other single contractions + (X N Y 1X N −1Y 2 − X N Y 2X N −1Y 1) {X1⋯XN −2Y3⋯YM } + all other double contractions + all triple contractions ⋯ . (1.63) We can expand the product by pairwise contracting one operator from each of the strings. It is easy to see that if there are N operators in the first and M operators in the second string, the expansion will contain products of length ∣N − M ∣ through N + M . Products of three or more normal-ordered operators can be evaluated by using Eq. (1.63) recursively. As examples, we evaluate {a† paq}{a† ras} ={a† pa† rasaq} + aqa† r{a† pas} ={a† pa† rasaq} + δqr{a† pas} (1.64) and {apaq}{a† ra† s} = − {a† ra† saqap} − aqa† r{a† sap} + aqa† s{a† rap} + apa† r{a† sap} − apa† s{a† rap} + aqa† rapa† s − apa† raqa† s = − {a† ra† saqap} − δqr{a† sap} + δqs{a† rap} + δpr{a† sap} − δps{a† rap} 24 + δqrδps − δprδqs . (1.65) We now turn to define the normal product of operators with respect to the Fermi vacuum ∣Φ0⟩ in an analogous way as for the vacuum ∣∅⟩. As discussed in Section 1.3.6, the killing condition for ∣Φ0⟩ turns into two separate conditions for hole and particle states, Eqs. (1.52) and (1.51). To make the analogy with the true vacuum, we can introduce new creation and annihilation operators for the Fermi vacuum by defining b† i ≡ ai bi ≡ a† i a ≡ a† b† a ba ≡ aa (1.66a) (1.66b) (1.66c) (1.66d) where i and a denote hole and particle indices, respectively, as introduced in Sec. 1.3.6. Equation (1.66a) states that annihilation of a fermion in the state i below the Fermi level is equivalent to creating a hole in the Fermi vacuum. Similarly, creating a particle in the state i below the Fermi level is analogous to destroying a hole in Eq. (1.66b). For the unoccupied states above the Fermi level the creation and annihilation operators remain unchanged. Since the killing condition in the new creators and annihilators is once again simply given by bp ∣Φ0⟩ = 0 (1.67) for all single-particle basis states, the discussion of the elementary contractions is exactly 25 the same as above, and the only non-vanishing contraction is bpb† q = δpq . Reverting to the original operators, this means for particle indices that bab† b = aaa† b = δab , and for hole indices bib† j = a† i aj = δij . (1.68) (1.69) (1.70) We can generalize these contractions to arbitrary single particle states with the help of occupation numbers: Then apa† q = (1 − np)δpq , a† paq = npδpq . (1.71) (1.72) 1.3.8 Hamiltonians and General Observables in Second-Quantized Form Let H be a general Hamiltonian for an A-body system, consisting of up to three-body interaction terms: H = H [1] + H [2] + H [3] . (1.73) In this work, we will consider both schematic Hamiltonians (usually without three-body parts) and realistic nuclear Hamiltonians consisting of kinetic terms as well as two- and 26 three-body interactions, for which H [1] = T [1] , H [2] = T [2] + V [2] , H [3] = V [3] . (1.74) (1.75) (1.76) There are one- and two-body kinetic energy terms because we will be using the intrinsic kinetic energy Tint = T − Tcm = (1 − 1 A ) ∑ i p2 i 2m − 1 A A ∑ i≠j pi ⋅ pj m . (1.77) Note that A is generally to be understood as a particle-number operator [38], although this subtlety doesn’t matter when we work with Hilbert spaces with fixed particle numbers as in this work. In second-quantized form, we can write the individual contributions in normal ordering with respect to the true vacuum as H [1] = ∑ pq ⟨p∣ h[1] ∣q⟩ a† paq , H [2] = 1 4 ∑ pqrs ⟨pq∣ h[2] ∣rs⟩ a† pa† qasar , H [3] = 1 36 ∑ pqrstu ⟨pqr∣ h[3] ∣stu⟩ a† pa† qa† rauatas (1.78) (1.79) (1.80) (note the usual convention for the index ordering in the matrix elements and the annihilation operators). Here, we have assumed that the matrix elements defining the two- and three- 27 body terms are antisymmetrized and obey ⟨pq∣ h[2] ∣rs⟩ = − ⟨qp∣ h[2] ∣rs⟩ = − ⟨pq∣ h[2] ∣sr⟩ = ⟨qp∣ h[2] ∣sr⟩ , (1.81) ⟨pqr∣ h[3] ∣stu⟩ = ⟨qrp∣ h[3] ∣stu⟩ = ⟨rpq∣ h[3] ∣stu⟩ = − ⟨qpr∣ h[3] ∣stu⟩ = − ⟨prq∣ h[3] ∣stu⟩ = − ⟨rqp∣ h[3] ∣stu⟩ = . . . (additional permutations) (1.82) Using a single Slater determinant ∣Φ⟩ as the reference state, we can use Wick’s theorem (cf. Section 1.3.7) to express the Hamiltonian in terms of normal–ordered operators: H = E + ∑ pq fpq{a† paq} + 1 4 ∑ pqrs Γpqrs{a† pa† qasar} + 1 36 ∑ pqrstu Wpqrstu{a† pa† qa† rauatas} . (1.83) To obtain the individual normal–ordered contributions in Eq. (1.83), let us first look at the one–body term (1.78): H [1] = ∑ pq ⟨p∣ h[1] ∣q⟩ a† paq = ∑ pq ⟨p∣ h[1] ∣q⟩ ({a† paq} + npδpq) = ∑ pq = ∑ pq ⟨p∣ h[1] ∣q⟩ {a† paq} + ∑ p ⟨p∣ h[1] ∣p⟩ np ⟨p∣ h[1] ∣q⟩ {a† paq} + ∑ i ⟨i∣ h[1] ∣i⟩ . (1.84) Thus, we obtain contributions to the normal-ordered zero- and one-body parts. For the two-body operator (1.79), Wick’s theorem (1.60) yields pa† a† qasar = {a† pa† qasar} + {a† pa† qasar} + {a† pa† qasar} + {a† pa† qasar} + {a† pa† qasar} 28 + {a† pa† qasar} + {a† pa† qasar} = {a† pa† qasar} + nqδqs{a† par} − nqδqr{a† pas} − npδps{a† qar} + npδpr{a† qas} + npnq (δprδqs − δpsδqr) , (1.85) (1.86) so the two–body operator can be decomposed into a sum of a normal-ordered two–body operator, a collection of one–body operators and two zero–body terms (scalars). Inserting this expansion into Eq. (1.79) H [2] = = 1 4 1 4 ∑ pqrs ∑ pqrs ⟨pq∣ h[2] ∣rs⟩ a† pa† qasar ⟨pq∣ h[2] ∣rs⟩ {a† pa† qasar} + 1 4 ( ∑ pqrs nqδqs ⟨pq∣ h[2] ∣rs⟩ {a† par} − nqδqr ⟨pq∣ h[2] ∣rs⟩ {a† pas} − npδps ⟨pq∣ h[2] ∣rs⟩ {a† qar} + npδpr ⟨pq∣ h[2] ∣rs⟩ {a† qas} + npnqδprδqs ⟨pq∣ h[2] ∣rs⟩ − npnqδpsδqr ⟨pq∣ h[2] ∣rs⟩) = 1 4 ∑ pqrs ⟨pq∣ h[2] ∣rs⟩ {a† pa† qasar} + + 1 4 1 4 ∑ pqi ∑ ij (⟨pi∣ h[2] ∣qi⟩ − ⟨pi∣ h[2] ∣iq⟩ − ⟨ip∣ h[2] ∣qi⟩ + ⟨ip∣ h[2] ∣iq⟩) {a† paq} (⟨ij∣ h[2] ∣ij⟩ − ⟨ij∣ h[2] ∣ji⟩) = 1 4 ∑ pqrs ⟨pq∣ h[2] ∣rs⟩ {a† pa† qasar} + ∑ pqi ⟨pi∣ h[2] ∣qi⟩ {a† paq} + 1 2 ∑ ij ⟨ij∣ h[2] ∣ij⟩ (1.87) where i and j refer to holes states (see Table 1.1). The calculation for the three-body part is analogous but much more tedious due to the much greater number of possible contractions 29 that need to be evaluated — without going into details, we see that it contributes to the zero-, one- and two-body parts of the normal-ordered H. Ultimately, we obtain E = ∑ i ⟨i∣ t ∣i⟩ + 1 2 ∑ ij ⟨ij∣ v ∣ij⟩ + 1 6 ∑ ijk ⟨ijk∣ w ∣ijk⟩ (1.88a) fpq = ⟨p∣ t ∣q⟩ + ∑ i ⟨pi∣ v ∣qi⟩ + 1 2 ∑ ij ⟨pij∣ w ∣qij⟩ Γpqrs = ⟨pq∣ v ∣rs⟩ + ∑ i ⟨pqi∣ w ∣rsi⟩ Wpqrstu = ⟨pqr∣ w ∣uts⟩ . (1.88b) (1.88c) (1.88d) All summations are limited to states that are occupied in the reference state ∣Φ⟩. We see that, as the name suggest, E is just the expectation value of H in the reference state, and f is the usual mean-field Hamiltonian, as used in the Hartree-Fock method, for instance. The zero-, one-, and two-body components of the Hamiltonian all encompass in-medium contributions from the 3N interaction. Other observables can be handled exactly like the Hamiltonian. The second-quantized form of a general k-body operator, can be written as O[k] = ( 2 ) 1 k! ∑ p1⋯pk q1⋯qk ⟨p1⋯pk∣ o[k] ∣q1⋯qk⟩ a† p1 ⋯a† pk aqk⋯aq1 (1.89) where we once again assume that the matrix elements have been antisymmetrized. Using Wick’s theorem (1.60), we can normal order such operators with respect to arbitrary Slater determinant reference states. 30 1.3.9 Full Configuration Interaction Theory In principle, we now have all ingredients for setting up the Hamiltonian matrix and trying to solve the (stationary) many-body Schr¨odinger equation through an exact diagonalization — this is the so-called Full Configuration Interaction (FCI) method [89], which is also known as No-Core Full Configuration (NCFC) or the No-Core Shell Model (NCSM) in nuclear physics [5, 74]. To perform an FCI calculation, one first needs to select an appropriate single-particle basis for describing bound systems, usually the spherical harmonic oscillator, and construct a complete orthonormal basis of Slater determinants by distributing the nucleons of a specific nucleus over the available single-particle levels in all possible configurations. As explained in Section 1.3.6, it is useful to organize the basis in terms of particle-hole excitations of a reference Slater determinant ∣Φ0⟩ to aid the convergence of the expansion for the exact eigenstates of H. Schematically, we have ∣Ψ⟩ = C0 ∣Φ0⟩ + ∑ ai C a i ∣Φa i ⟩ + ∑ abij C ab ij ∣Φab ij ⟩ + ⋯ ≡ (C0 + ˆC) ∣Φ0⟩ (1.90) where the ellipsis three-particle-three-hole (3p3h) and higher excitations, and where we have introduced the so-called correlation operator ˆC = ∑ ai i a† C a aai + 1 4 ∑ abij Our requirement of orthonormality gives ij a† C ab aa† bajai + ⋯ . (1.91) ⟨Ψ0 ∣ Ψ0⟩ = ∑ ai ∣C a i ∣2 + ∑ abij ∣C ab ij ∣2 + ⋯ ≡ ∑ I ∣CI∣2 = 1 , (1.92) 31 where the collective index I runs over all sectors of the many-body basis. Expanding the eigenstates in the many-body Schr¨odinger Equation (1.1) in terms of the Slater determinant basis configurations, we obtain H ∑ I CI ∣ΦI⟩ = E ∑ I CI ∣ΦI⟩ Projecting Eq. (1.93) onto ⟨ΦI ′∣, we see that ⟨ΦI ′∣ H ∑ I CI ∣ΦI⟩ = E ⟨ΦI ′∣ ∑ I CI ∣ΦI⟩ By rearranging the formula, we obtain ∑ I CI ⟨ΦI ′∣ H ∣ΦI⟩ = E ∑ I CI ⟨ΦI ′∣ ∣ΦI⟩ = E ∑ I CIδII ′ = ECI ′ (1.93) (1.94) (1.95) (1.96) (1.97) which precisely corresponds to the eigenvalue problem in linear algebra, commonly expressed as HC = EC. The elements of the Hamiltonian are given by HII ′ = ⟨ΦI∣ H ∣ΦI ′⟩ , and they can be evaluate with the help of Wick’s theorem: Note that aa† a† b⋯ajai = {a† aa† b⋯ajai} (1.98) 32 because we cannot contract particle (a, b, . . .) and hole (i, j, . . .) indices, so we can write the matrix elements as HII ′ = ⟨Φ0∣ {a† i1 ⋯a† ik aak⋯aa1}H{a† b1 ⋯a† bk′ ajk′ ⋯aj1} ∣Φ0⟩ , (1.99) where I and I ′ are understood as kpkh and k′pk′h basis configurations, respectively. In a practical many-body calculation, we first need to truncate the single-particle basis, which leads to a finite number of Slater determinants and thus a Hamiltonian matrix with finite size. In principle, eigenstates can then be obtained by directly diagonalizing the trun- cated matrix, usually with Lanczos [59] or Arnoldi methods [3], since we are mostly interested in a few low-lying states. However, the basis dimension for the matrices and eigenvectors scales as (N A) = N ! A!(N −A)! with the single-particle basis size N and the nucleon number A, hence the diagonalization quickly becomes impractical even for modern supercomputers as A grows. Therefore, we need to truncate the basis of Slater determinants. The No-Core Shell Model (NCSM) in nuclear physics uses an (oscillator) energy based truncation instead that can guarantee an exact factorization of the center of mass and in- trinsic wave functions, but does not change the exponential scaling of the method [5]. Another option is to restrict the type of excitations that are included in the basis. If only singles and doubles excitations are considered, for example, then we get the configuration singles doubles (CISD) approximation ∣ΨCISD 0 ⟩ = C0 ∣Φ0⟩ + ∑ ia C a i ∣Φa i ⟩ + ∑ ijab C ab ij ∣Φab ij ⟩ , (1.100) for which the basis dimension scales as O(N 2 p N 2 h) with the number of hole (Nh = A) and the 33 number of particle states (Np = N − Nh = N − A). Lanczos or Arnoldi iterations can then be performed at a cost of O(N 4 p N 4 h), but there are alternative methods with similar scaling that are more efficient for incorporating correlations into the target many-body state, like Coupled Cluster [89] or the In-Medium Similarity Renormalization Group [39, 41], which will be discussed in later chapters of this work. Finding optimal ways to capture many-body correlations is a widely studied topic (see, e.g., [40] and references therein), and there is presently no one-size-fits-all solution. For specific applications, the choice of approximation depends on the strengths and weaknesses of each method. 34 Chapter 2 The Similarity Renormalization Group Nuclear many-body calculations are challenging because of the influence of correlation effects that arise from the strong repulsion of the nuclear force at small distances as well as the tensor force, as mentioned in chapter 1. This leads to a coupling of many-body states with high and low momentum. To address this issue, we can use renormalization group (RG) techniques to “tame” these correlation effects. RG methods are naturally complementary to EFT approaches, and in our case, to the sequence of EFTs for the strong interaction. This is because they enable the systematic adjustment of resolution scales and cutoffs in these theories, allowing for the connection of various levels within the EFT hierarchy. Additionally, RG methods also offer a range of diagnostic tools for checking the fundamental consistency of power counting schemes in EFTs. In the field of nuclear many-body theory, the most commonly used method is the Sim- ilarity Renormalization Group (SRG), which will be discussed in more detail in this chap- ter. This technique differs from the Wilsonian [109] RG approach, which eliminates high- momentum degrees of freedom through decimation. The SRG method separates the physics of low and high momenta through continuous unitary transformations. It is important to note that this concept is not limited to RG applications, as it can be used to modify a many- 35 body Hamiltonian or other relevant observables according to specific requirements, such as extracting eigenvalues [37, 39] or imposing certain operator structures [10, 13, 49, 93, 94], which will be discussed in the next chapter. 2.1 Flow Equation We utilize the method of SRG transformations to “soften” the interactions between nu- cleons. The transformation process involves a series of unitary transformations, which allows for the evaluation of the Hamiltonian’s evolution or flow as a function of the parameter s, Hs = UsH0U † s ≡ Trel + Vs (2.1) where Trel is the relative kinetic energy1, , H0 is the initial Hamiltonian, and Us is the unitary transformation. Note that we make the choise to absorb induced operators from the evolution of Trel into Vs. Instead of making an ansatz for Us, we implement SRG transformations through the use of the flow equation formalism, as introduced by Wegner [100, 101]. We find the flow equation by taking the derivative of Eq. (2.1) with respect to s, dHs ds = = = dUs ds dUs ds dUs ds H0U † s + UsH0 s UsH0U † U † U † s Hs + HsUs dU † s ds s + UsH0U † dU † s ds s Us dU † s ds (2.2) 1Since all fundamental interactions are translationally invariant, the center-of-mass kinetic energy (Tcm) commutes with the other operators in Eq.(2.1). 36 where we have plugged in 1 = UsU † s = U † s Us in the second line, and we also have used which implies d ds (1) = d ds (UsU † s ) = 0 dUs ds U † s = −Us dU † s ds Defining the anti-Hermitian operator Eq. (2.2) can be rewritten as ηs ≡ dUs ds U † s = −ηs , dHs ds = ηsHs + Hsη† s = ηsHs − Hsηs = [ηs, Hs] (2.3) (2.4) (2.5) (2.6) The SRG flow equation for the Hamiltonian describes the way a system changes over time when it is subject to the dynamical generator ηs. Note that because of the unitary nature of the transformation, the Hamiltonian’s spectrum is conserved unless we need to introduce approximations. By integrating the flow equation (2.6) for s → ∞, we can achieve the desired transforma- tion of the Hamiltonian by selecting an appropriate generator ηs. Wegner originally proposed a type of generators known (henceforth referred to as Wegner generators), given by ηs = [H d s , H od s ] (2.7) 37 which is constructed by dividing the Hamiltonian into a diagonal part (H d s ) we want to keep as s → ∞, and an off-diagonal part (H od s ) we want to suppress. Several options for splitting the Hamiltonian are available for consideration [1, 9, 51, 53], and one can even generalize this idea and consider the commutator of an arbitrary Hermitian operator Gs with Hs, e.g., for unmixing symmetries [49]. In the following, we use H d s = Trel, a constant, for which the flow equation becomes dHs ds = [[Trel, Hs] , Hs] (2.8) The Eq. (2.8) represent operator equations that are independent of any specific basis. The equation can be evaluated numerically by projecting onto a suitable basis and converting it into a matrix differential equation. In the case of two-nucleon system, a partial-wave momentum basis is a convenient, accurate, and straightforward choice, since our generator is based on Trel. 2.2 SRG Decoupling One of the major consequence of the SRG is that the Hamiltonian will be driven to a diagonal or band-diagonal form in a chosen basis. To see this, we consider the Wegner generator and split the Hamiltonian matrix into diagonal and off-diagonal pieces (2.9) H d s = diag(Hs) H od s = Hs − H d s 38 Suppressing the s dependence of the quantities, Eq. (2.6) can then be written as d ds ⟨i∣ H ∣j⟩ = ⟨i∣ [η, H] ∣j⟩ (⟨i∣ η ∣k⟩ ⟨k∣ H ∣j⟩ − ⟨i∣ H ∣k⟩ ⟨k∣ η ∣j⟩) = ∑ k (⟨i∣ η ∣k⟩ ⟨k∣ H d + H od ∣j⟩ − ⟨i∣ H d + H od ∣k⟩ ⟨k∣ η ∣j⟩) = ∑ k = −(Ei − Ej) ⟨i∣ η ∣j⟩ + ∑ k (⟨i∣ η ∣k⟩ ⟨k∣ H od ∣j⟩ − ⟨i∣ H od ∣k⟩ ⟨k∣ η ∣j⟩) (2.10) where we have used the completeness of the basis, 1 = ∑k ∣k⟩ ⟨k∣, in the second equality above and, in the third equality, we have used that H d ∣k⟩ = Ek ∣k⟩ by construction. Then we also have ⟨i∣ η ∣j⟩ = ⟨i∣ [H d, H od] ∣j⟩ = (Ei − Ej) ⟨i∣ H od ∣j⟩ . (2.11) Inserting this into the Eq. (2.10), we obtain d ds ⟨i∣ H ∣j⟩ = −(Ei − Ej)2 ⟨i∣ H od ∣j⟩ + ∑ k (Ei + Ej − 2Ek) ⟨i∣ H od ∣k⟩ ⟨k∣ H od ∣j⟩ (2.12) For the diagonal elements (i = j), Eq. (2.12) becomes d ds ⟨i∣ H ∣i⟩ = 2 ∑ k (Ei − Ek)∣⟨i∣ H od ∣k⟩∣2 (2.13) 39 where ⟨i∣ H od ∣i⟩ = 0, and d ds ∑ i ∣⟨i∣ H ∣i⟩∣2 = 2 ∑ i ⟨i∣ H ∣i⟩ d ds ⟨i∣ H ∣i⟩ Ei(Ei − Ek)∣⟨i∣ H od ∣k⟩∣2 Ei(Ei − Ek)2∣⟨i∣ H od ∣k⟩∣2 (2.14) = 2 ∑ i,k = 4 ∑ i≠k (Ei − Ek)2∣⟨i∣ H od ∣k⟩∣2 = 2 ∑ i≠k ≥ 0 . Due to unitary, the trace of a matrix is fixed under the evolution Tr H 2 = const. = ∑ i,j ∣⟨i∣ H ∣j⟩∣2 = ∑ i ∣⟨i∣ H ∣i⟩∣2 + ∑ i≠j ∣⟨i∣ H ∣j⟩∣2 Therefore, we obtain d ds ∣⟨i∣ H ∣j⟩∣2 = − ∑ i≠j d ds ∣⟨i∣ H ∣i⟩∣2 ∑ i ≤ 0 . (2.15) (2.16) This means the off-diagonal matrix elements will decrease and, the transformation caused by η effectively suppresses H od. Thinking further, this means that the norm of the absolute value of H od will become small for some s0. For even larger values s > s0, we can then disregard the second term of the flow equation (2.12) compared to the first term, and this implies dEi ds = d ds ⟨i∣ H d ∣i⟩ = 2 ∑ k (Ei − Ek)∣⟨i∣ H od ∣k⟩∣2 ≈ 0 (2.17) 40 and d ds ⟨i∣ H od ∣j⟩ ≈ (Ei − Ej)2 ⟨i∣ H od ∣j⟩ (2.18) respectively. This means that the diagonal elements of the matrix will remain nearly constant in the asymptotic region, Es ≈ Es0, s > s0 , (2.19) which allows us to integrate Eq. (2.18): ⟨i∣ H od s ∣j⟩ ≈ ⟨i∣ H od s0 ∣j⟩ e−(Ei−Ej )2(s−s0), s > s0 . (2.20) It is revealed that the characteristic decay scale of each element is dependent on the square of the energy difference between the states it connects, e.g. (Ei−Ej)2: States with higher energy differences will be decoupled before those with lower energy difference, thereby ensuring that the Wegner generator produces a proper RG transformation. 2.3 Evolution of Nuclear Interactions For the sake of simplicity, let us consider a nuclear Hamiltonian that consists of the intrinsic kinetic energy and a two-nucleon interaction (cf. Section 1.3.8): Hint = Tint + V [2] . (2.21) Tint (Eq. (1.77)) can be written as a pure two-body operator Tint = 2 A ∑ i Λ1 > Λ2 and integrating out high-momentum modes, effectively renormalizing it into smaller and smaller blocks in momentum space. This method is also known as the Vlow-k approach and was first introduced in nuclear physics during the early 2000s [10, 11]. On the other hand, Figure 2.1b displays the SRG evolution of the N N interaction to a band-diagonal form using the flow equation (2.8), where the generator is constructed from the relative kinetic energy Trel in the two-nucleon system [9, 10]. Here, 43 Figure 2.1: The diagram depicts two different RG evolutions of N N potentials in momentum space. In the Vlow-k approach (a), the interaction is renormalized into blocks of decreasing cutoff Λ. In contrast, the SRG evolution (b) transforms the interaction to a band-diagonal form with decreasing width λ. In the diagram, k and k′ designate the relative momenta of the initial and final states, respectively. At each Λi and λi, the matrix elements outside of the corresponding blocks or bands are negligible. Therefore, it can be inferred that high and low momentum states are effectively decoupled [10]. the SRG evolution of a system has been parameterized by a momentum with dimensions of λ = s−1/4, rather than the usual parameter s, which is an indicator of the width of the band in momentum space: ∣q′ − q∣ ≲ λ . (2.26) This suggests low- and high-momenta are decoupled as λ is decreased through the SRG. Physically, this implies that long-ranged components of the N N interactions are explicitly resolved, while the short-range components can be replaced by contact interactions [11, 33, 47, 61]. Thus, we can think of λ as a momentum resolution scale. When an N N interaction is evoved through an RG process, it leads to a universal long-range interaction, namely one- pion exchange (OPE), which is valid in the range of 1.5 fm−1 to 2.5 fm−1. Further evolution to lower values of λ removes pieces of OPE and eventually results in a pion-less theory that 44 Λ0Λ1Λ2k’kλ0λ1λ2k’k(a)(b) is purely parameterized in terms of contact interactions. Using the orthogonality and completeness relations for the states (2.25) ⟨qLSJM T MT ∣ ∣q′L′S′J ′M ′T ′M ′ T ⟩ = δ(q − q′) qq′ δLL′δSS′δJJ ′δM M ′δT T ′δMT M ′ T , (2.27a) 1 = ∑ LSJM T MT ∞ ∫ 0 dqq2 ∣qLSJM T MT ⟩ ⟨qLSJM T MT ∣ , (2.27b) we can represent the flow equation (2.24) in partial-wave form as (− λ5 4 ) d ⟨qL∣ v ∣q′L′⟩ dλ = − (q2 − q′2)2 ⟨qL∣ v ∣q′L′⟩ ∞ + ∑ ¯L ∫ 0 dpp2(q2 + q′2 − 2p2) ⟨qL∣ v ∣p ¯L⟩ ⟨p ¯L∣ v ∣q′L′⟩ (2.28) where we used ̵h = m = 1 and a prefactor −λ5/4 appears due to the variables change from the flow parameter s to the resolutions scale λ. By discretizing the relative momentum variable using uniform or Gaussian quadrature meshes [81] , we can convert this integro-differential equation into a matrix flow equation. For example, the matrix elements of the relative kinetic energy operator can be represented as ⟨qiL∣ t ∣qjL′⟩ = q2 i δqiqj δLL′ and the discretization transforms the continuous integration into the summation ∞ ∫ 0 dqq2 ⇒ ∑ i wiq2 i where the weights wi depend on the choice of mesh. To simplify the calculation, we incor- porate the weights and q2 from the integral measure into the interaction matrix element, 45 defining ⟨qiL∣ ¯v ∣qjL′⟩ ≡ √ wiwjqiqj ⟨qiL∣ v ∣qjL′⟩ . (2.29) As a result, the discretized flow equation can be expressed as d dλ ⟨qiL∣ ¯v ∣qjL′⟩ = − 4 λ5 ⟨qiL∣ [[t, ¯v], t + ¯v] ∣qjL′⟩ . (2.30) As an illustration of a practical application, snapshots of the SRG transformation of the chiral N3LO N N interactions by Entem and Machleidt with a cutoff Λ = 500 MeV [22, 68] are shown in Figure 2.2. The left column displays the matrix elements of the initial interaction s = 0, or λ = ∞, in the 3S1 partial wave, and the corresponding deuteron wave function. The middle column and the right column are the interaction evolved to 3 fm−1 and then to 2 fm−1, where the off-diagonal matrix elements are suppressed and the interaction is primarily confined to a block of states with ∣q∣ ∼ 2 fm−1. Furthermore, it is noted that the chiral interaction initially has much weaker off-diagonal matrix elements, while the AV18 matrix elements extend up to ∣q∣ ∼ 20 fm−1. In nuclear physics terminology, AV18 is considered as a harder interaction than the chiral interaction due to its repulsive core. The positive impact of low-momentum potentials in nuclear many-body calculations is illustrated in Figures 2.3a and 2.3b. The figures show ground-state energies of 4He and 6Li from exact diagonalization of the Hamiltonian matrix in the No-Core Shell Model (NCSM), one of the nuclear FCI approaches discussed in Section 1.3.9. The original potential is once again the chiral N N interaction by Entem and Machleidt that we discussed above, and it is now supplemented with a chiral 3N interaction — for details, we refer to [27, 50]. Although this interaction is rather soft in nuclear physics terminology, the NCSM convergence with respect to the harmonic oscillator (HO) many-body basis is slow, all while matrix dimensions 46 (a) λ = ∞ (b) λ = 3 fm−1 (b) (c) λ = 2 fm−1 Figure 2.2: The SRG evolution of the chiral N3LO nucleon-nucleon interaction by Entem and Machleidt [22, 68], with initial cutoff Λ = 500 MeV. In the top row, we display the momentum-space matrix elements of the interaction in the 3S1 partial wave at different values of the SRG resolution scale. The initial interaction is shown in the left column at s = 0, or λ = ∞. In the bottom panel, we display S− and D−wave components of the deuteron wave function that results from solving the Schr¨odinger equation with the corresponding SRG- evolved interaction [41]. are growing rapidly with both Nmax and the number of nucleons, as shown in Fig. 2.3c. However, SRG evolution significantly improves convergence and allows reliable extrapolation of partly converged result to the Nmax → ∞ limit, reducing the computational cost by several orders of magnitude. Because of the SRG softening of the interaction and the use of techniques like importance sampling [84, 87], exact diagonalization approaches can today reach nuclei with mass A ∼ 20 [40], but the fast growth of the basis with the number of nucleons remains an obstacle in calculating even larger nuclei. 2.4 Induced Interactions While SRG decoupling of low- and high-momentum components of nuclear interactions improves the convergence of many-body calculations, there is a price to pay for in practical 47 02468100.00.10.20.30.4r(cid:64)fm(cid:68)ΦL(cid:72)r(cid:76)(cid:64)fm(cid:45)3(cid:144)2(cid:68)L(cid:61)0L(cid:61)202468100.00.10.20.30.4r(cid:64)fm(cid:68)ΦL(cid:72)r(cid:76)(cid:64)fm(cid:45)3(cid:144)2(cid:68)L(cid:61)0L(cid:61)202468100.00.10.20.30.4r(cid:64)fm(cid:68)ΦL(cid:72)r(cid:76)(cid:64)fm(cid:45)3(cid:144)2(cid:68)L(cid:61)0L(cid:61)2 (a) (b) (c) Figure 2.3: The convergence of the No-Core Shell Model (NCSM) calculations of the ground- state energies of 4He (a) and 6Li (b) is analyzed as a function of the basis size parameter Nmax at different SRG resolution scales λ. Panel (c) illustrates the growth of the NCSM basis dimension with Nmax for various nuclei. See Refs. [27, 50] for details. applications: The evolution induces three-body and higher many-body forces. To see this, we again consider the Hamiltonian in second-quantized form (see Section 1.3.8): H = ∑ pq T [1] pq a† paq + 1 4 ∑ pqrs (T [2] + V [2]) pqrs pa† a† qasar + 1 36 ∑ pqrstu V [3] pqrstua† pa† qa† rauatas . (2.31) To understand how the SRG generates many-body forces in an A-body space, we do not need to solve the flow equation (2.24) using the Hamiltonian given by Eq. (2.31) — we can demonstrate the problem by considering only the two-body term: dVs ds = ⎡ ⎢ ⎢ [∑ ⎢ ⎣ pq T [1] pq a† paq ∑ pqrs T [2] pqrsa† pa† qasar, ∑ tuvw tuvwa† V [2] t a† uawav] , ∑ ijkl ijkla† V [2] i a† jalak ⎤ ⎥ ⎥ ⎥ ⎦ + . . . . (2.32) Without explicitly calculating the above, let’s consider the intermediate steps, the calculation of commutators: For the commutator of a one- and a two-body operator, application of 48 2468101214161820Matrix Size [Nmax]−29−28−27−26−25−24−23−22−21−20Ground State Energy [MeV]Originalexpt.VNN = N3LO (500 MeV)Helium-4Softenedwith SRGVNNN = N2LOground-state energy24681012141618Matrix Size [Nmax] −36−32−28−24−20−16−12−8−40481216Ground-State Energy [MeV]Lithium-6expt.Originalh- Ω = 20 MeVSoftened with SRGground-state energyVNN = N3LO (500 MeV)VNNN = N2LOλ = 1.5 fm−1λ = 2.0 fm−102468101214Nmax101102103104105106107108109Matrix dimension4He6Li12C16O Wick’s theorem (1.63) yields [a† paq, a† ra† sauat] =a† paqa† ra† sauat − a† ra† sauata† paq =a† pa† ra† sauataq + δqra† pa† sauat − δqsa† pa† rauat − a† pa† ra† sauataq − δpta† ra† sauaq + δpua† ra† sataq =δqra† pa† sauat − δqsa† pa† rauat − δpta† ra† sauaq + δpua† ra† sataq , (2.33) so we see that while each individual product in the commutator creates a three-body term (cf. Section 1.3.8), those terms cancel and we once again obtain a two-body operator, albeit one that appears more complicated. For the commutator of two two-body operators, however, [a† pa† qasar, a† t a† uawav] =a† pa† qasara† t a† uawav − a† t a† uawava† pa† qasar =a† pa† qa† t a† uawavasar − δrta† pa† qa† uawavas + other single contractions + (δrtδsu − δruδst) a† pa† qawavas − a† pa† qa† t a† uawavasar + δpva† t a† ua† qawasar + other single contractions − (δpvδqw − δpwδqv) a† t a† uasar , (2.34) i.e., the four-body operators cancel, but a new three-body operator remains. This would be true even if we started with a pure N N interaction, and as soon as after a single integration step. The induced 3N interaction will then induce up to 5N interactions in the next inte- gration step, and so on — in general, the SRG flow induces up to A-nucleon interactions. However, these interactions are only potentially significant when analyzing A-nucleon sys- tems [34, 35, 50, 51], and the expectation is that their importance diminishes with increasing A, so that the expected hierarchy of interactions is preserved, namely that V [2] > V [3] > V [4]. 49 Figure 2.4: The ground state energy of 3H as a function of the SRG resolution scale λ for chiral N N and N N + 3N interaction [34]. Three different typles of interactions are considered: N N − only, N N + 3N − induced, and N N + 3N − full. The experimental binding energy [99] is represented by a black dotted line. Figure adapted from Ref. [41], with original data courtesy of K. Hebeler. For a practical calculation of 3H, for example, truncating the SRG flow equations the three-body level, should preserve all observables. The 3H ground-state energies for different chiral N N and N N + 3N interactions, are presented in Figure 2.4. We see that the energy changes significantly when we truncate the generator and evolving Hamiltonian at the two- body level (N N − only), while it remains constant when the evolution is performed at the three-body level and induced 3N interactions are included (N N + 3N − induced). The same remains true when we include an initial 3N force whose parameters have been fit to reproduce the experimental 3H ground-state energy (N N + 3N − full). 50 1.5234571015λ [fm−1]−8.5−8.4−8.3−8.2−8.1E [MeV]NN-onlyNN + 3N-inducedNN + 3N-full0.00010.0010.010.1s [fm4]−8.5−8.4−8.3−8.2−8.1Expt.550/600 MeV Chapter 3 The In-Medium SRG In our discussion of the SRG in chapter 2, we mentioned that we can contemplate utilizing the SRG methodology not solely for preprocessing nuclear interactions but also for imposing other desired structures on the Hamiltonian and other observables — we could even use it to calculate eigenvalues and eigenstates if we were to implement the SRG for the Hamiltonian matrix and chose H d to be the diagonal of that matrix (see Section 2.2). This is not practical, of course, because of the exponential growth of the many-body basis (cf. Section 1.3.9) with the single-particle basis size and particle number. However, we can actually implement the SRG evolution directly at the level of the operator, and this leads us to the so-called In-Medium SRG (IMSRG), which is the subject of this chapter. 3.1 The Role of Normal Ordering In order to implement SRG flow equations on the operator level, we first need to choose a set of “basis” operators. Ideally, we would want this set to be complete under the application of commutators, i.e., the operators would form an algebra, but unfortunately this is not possible in practice. Because our goal is to compute energies and other observables of many-body systems, a 51 basis of second-quantized operator strings is a natural choice, and we could write H = ∑ pq Hpqa† paq + 1 4 ∑ pqrs Hpqrsa† pa† qasar + ⋯ η = ∑ pq ηpqa† paq + 1 4 ∑ pqrs ηpqrsa† pa† qasar + ⋯, (3.1) (3.2) dropping the particle-rank superscripts because it is clear how many particles the operators act on from the number of indices. In Section 2.4, we saw that each evaluation of the commutator on the RHS of the operator flow equation (2.6) leads to an increase in the particle rank of H(s), ηpqrsHtuvw[a† pa† qasar, a† t a† uawav] ∑ pqrstuvw = ∑ pqrstuv ηpqrvHvstua† pa† qa† rauatas + 3N terms + N N terms . (3.3) (3.4) It should be noted that there are no induced 4N interaction, and that commutators involving at least one one–body operator do not result in any change in the particle rank. We cautioned that the impact of truncating these terms to “close” the flow equations can lead to significant truncation errors, and that these errors are expected to become more substantial in in systems comprising A ≥ 3 particles. Luckily, we also have a potential solution to this issue at our disposal: normal ordering. In Section 1.3.7, we saw that the process of normal ordering an operator with respect to a reference state, e.g., a Slater determinant, shifts parts of any second-quantized operators into terms of lower particle rank. Likewise, Wick’s theorem 1.63 tells us that a large number of terms with contractions appear in the expansion of products — and therefore commutators — of normal ordered operators. For this reason, we choose to work with a basis of normal- 52 ordered operators, B ≡ {{a† paq}, {a† pa† qasar}, {a† pa† qasar}, ⋯} . (3.5) The discussion in Section 1.3.8 then tells us how to represent the relevant operators for the IMSRG in the basis B. We still must introduce truncations to the operator expressions, since each integration step in the IMSRG flow equations will technically induce new normal ordered operators of higher-body rank. We usually close the flow equations by truncating all operators consis- tently at a given operator level k, including the initial normal ordered Hamiltonian and observables of interest. This defines the so-called NOkB approximations. Empirical evidence suggests that in the case of low–momentum interactions, neglecting the normal–ordered three-body component of the Hamiltonian results in a deviation of only 1-2% in the energies of the ground and absolute excited states of light and medium–mass nuclei [7, 29, 31, 86, 87]. This NO2B approximation will serve as our primary workhorse in the following. It offers an efficient approach to incorporate the effects of 3N forces in nuclear many–body calculations without incurring the computational expense of explicitly considering three-body operators. 3.2 In-Medium SRG Flow Equations 3.2.1 The IMSRG(2) Scheme As discussed in Section 3.1, we close the IMSRG flow equations by truncating all induced operators above a given rank k, consistent with the NOkB approximation for representing the initial Hamiltonian. In the IMSRG context, this defines the IMSRG(k) truncation scheme. 53 Because of the empirical success of NO2B approximation for low-momentum N N +3N Hamil- tonians, where the omission of normal–ordered 3N term in exact calculations causes only small deviations of about 1–2% in ground-state energies of in oxygen, calcium, and nickel isotopes [6, 7, 86], we adopt k = 2. Let’s assume, then, that for each flow parameter s η(s) ≈ η[1](s) + η[2](s) H(s) ≈ E(s) + f (s) + Γ(s) d ds H(s) ≈ d ds E(s) + d ds f (s) + d ds Γ(s) (3.6) (3.7) (3.8) This is the so–called IMSRG(2) truncation, the foundation of all the results presented in the remaining part of this dissertation. It is similar to Coupled Cluster with Singles and Doubles (CCSD) [30, 64, 89] and the ADC(3) methods used in Self–Consistent Green’s Function Theory [4, 90] to describe the same level of many–body correlations. Therefore, we anticipate similar convergence properties and results for observables obtained with the three techniques. By inserting Eq. (3.6) − (3.8) into the operator flow equation (2.6) and calculating com- mutators from Section 1.3.7, we obtain dE ds dfpq ds (np − nq)ηpqfqp + = ∑ ab 1 2 ∑ pqrs ηpqrsΓrspqnpnq ¯nr ¯ns (3.9) = ∑ r (1 + Ppq)ηprfrq + ∑ rs (nr − ns)(ηrsΓsprq − frsηsprq) + 1 2 (nrns¯nt + ¯nr ¯nsnt)(1 + Ppq)ηtprsΓrstq ∑ rst (3.10) 54 dΓpqrs ds {(1 − Ppq)(ηptΓtqrs − fptηtqrs) − (1 − Prs)(ηtrΓpqts − ftrηpqts)} = ∑ t + 1 2 (1 − nt − nu)(ηpqtuΓturs − Γpqtuηturs) ∑ tu (nt − nu)(1 − Ppq)(1 − Prs)ηtpurΓuqts + ∑ tu (3.11) where we have suppressed the s–dependence has been suppressed for brevity and intro- duced the permutation symbol Pij to interchange the indices in any expression Pijg(⋯, i, ⋯, j) ≡ g(⋯, j, ⋯, i) . (3.12) As in chapter 1, na ∈ {0, 1} designates the particle occupation number and we introduce ¯na = 1 − na the hole occupation number for convenience. Since the reference state is fixed, they do not change during the evolution. In Eqs. (3.9)–(3.11), we note that the occupation number factors effectively restrict the summations to particular combinations of particle and/or hole indices. The terms on the RHS of the energy flow equation (3.9), for example, are restricted to one-body matrix elements between particle and hole (ph) states, and two- body matrix elements with a pphh structure. To calculate nuclear ground–state energies, we integrate equations (3.9-3.11) from s = 0 to s → ∞. We begin with the initial components of a normal–ordered Hamiltonian in NO2B, which are given by equations (1.88). By integrating the flow equations, we incorporate many–body correlations directly into the normal–ordered Hamiltonian. Looking at Eq. (3.9) again and noting that η is constructed from H, we notice a similarity to second-order en- ergy corrections, except that here they would be evaluate for the flowing Hamiltonian H(s). Similarly, the RHSs of Eqs. (3.10) and (3.11) have the structure of many-body corrections 55 to f and Γ, but since H(s) rather than H(0) are inserted in these expressions, the IMSRG evolution constructs highly nested versions of these corrections — from a Many-Body per- turbation Theory perspective, certain classes of terms are included to all orders, which is why the IMSRG is a non-perturbative method [39]. 3.2.2 Computational Scaling The computational scaling of the IMSRG(2) scheme is dominated by the two–body flow equation, which naively requires O(N 6) operations, putting it in the same category as CCSD and ADC(3). However, matrix products can be used to express large portions of the flow equations, allowing for the use of optimized linear algebra libraries. To further reduce com- putational cost, particle and hole states can be distinguished, as Nh, the number of hole states, is typically much smaller than Np, the number of particle states. The best scaling achievable in IMSRG(2) depends on the choice of generator η: If the one– and two–body parts only consist of ph and pphh type matrix elements and their Hermitian conjugates, the scaling is reduced to O(N 2 hN 4 p ), matching the cost of solving CCSD amplitude equations. We will come back to this discussion in later parts of this dissertation. 3.3 Decoupling After establishing the IMSRG flow equations, it is necessary to determine the decoupling approach, which involves separating the Hamiltonian into diagonal components to be pre- served and off–diagonal components to be suppressed, as described in Sections 2.2 and 2.3. To accomplish this, we utilize the matrix representation of the Hamiltonian in an A-body Slater determinant basis (see Section 1.3.9). It should be noted that the Hamiltonian ma- 56 Figure 3.1: Schematic view of IMSRG decoupling in the many-body Hilbert space spanned by a Slater determinant reference ∣Φ⟩ and its particle-hole excitations ∣Φa⋯ i⋯ ⟩. trix is not actually constructed in this representation, because its size is prohibitive. Our Slater determinant basis comprises the reference determinant ∣Φ⟩ and all of its particle-hole excitations, as discussed in Sections 1.3.6 and 1.3.9: ∣Φ⟩ , {a† aai} ∣Φ⟩ , {a† aa† bajai} ∣Φ⟩ , ⋯ (3.13) By utilizing Wick’s theorem, it can be demonstrated that the particle–hole excited Slater determinants are orthogonal to the reference state and each other (see Sections 1.3.6–1.3.9). In the Hilbert space spanned by this basis, the initial Hamiltonian (in NO2B approximation) has the structure illustrated in the left panel of Figure 3.1. Its matrix representation has a band–diagonal form and can potentially couple excitations involving npnh and (n ± 2)p(n ± 2)h. In order to proceed, we must decompose the Hamiltonian into suitable diagonal and off– diagonal parts on the operator level, as outlined in references [56–58]. A broad interpretation of diagonality is not recommended, as it may result in strong in–medium 3N, ⋯ interactions, 57 (cid:104)i|H(0)|j(cid:105)(cid:104)i|H(∞)|j(cid:105)s→∞−→|Φabcijk(cid:105)|Φabij(cid:105)|Φai(cid:105)|Φ(cid:105)|Φ(cid:105)|Φai(cid:105)|Φabij(cid:105)|Φabcijk(cid:105)|Φabcijk(cid:105)|Φabij(cid:105)|Φai(cid:105)|Φ(cid:105)|Φ(cid:105)|Φai(cid:105)|Φabij(cid:105)|Φabcijk(cid:105) which can violate the truncation of the IMSRG(2). To avoid this, we employ a minimal decoupling scheme that exclusively aims to decouple the one–dimensional block consisting of the reference state from all particle–hole excitations. The right panel of Figure 3.1 provides an illustration of the descired outcome of the evolution. In the ideal case, if we could implement the minimal decoupling scheme without any approximations, we would extract a single eigenvalue and eigenstate of the many–body Hamiltonian for the nucleus of interest in the limit as s → ∞. The eigenvalue would be equivalent to the zero–body piece of H(∞), or simply E(∞), while the eigenstate could be obtained by applying the unitary IMSRG transformation to the reference state, U †(∞) ∣Φ⟩. However, in practice, truncations are inevitable, and we can only obtain an approximate eigenvalue and mapping. It should be noted that the choice of reference state plays a crucial role in determining which eigenvalue and eigenstate of the Hamiltonian can be extracted using our minimal decoupling scheme. An empirical guideline is that the IMSRG flow will connect the reference state to the eigenstate that has the largest overlap with it. In general, if we aim to obtain the exact ground state, a Hartree-Fock Slater determinant is typically the best choice, as it minimizes both the absolute energy and the correlation energy. By utilizing Wick’s theorem to analyze the matrix elements between the reference state and its excitations, we can observe that the Hamiltonian couples the 0p0h block to 1p1h excitations block through specific matrix elements and their Hermitian conjugates: ⟨Φ∣ H{a† aai} ∣Φ⟩ = E ⟨Φ∣ {a† aai} ∣Φ⟩ + ∑ pq fpq ⟨Φ∣ {a† paq}{a† aai} ∣Φ⟩ + 1 4 ∑ pqrs Γpqrs ⟨Φ∣ {a† pa† qasar}{a† aai} ∣Φ⟩ = ∑ pq fpqδpiδaqnp¯nq 58 = fia (3.14) where the contributions arising from the zero–body and two–body components of the Hamil- tonian vanish since they represent the expectation values of normal-ordered operators in the reference state — in the latter case, because we cannot fully contract a one- and a two-body operator (cf. Section 1.3.7). Similarly, the 0p0h and 2p2h blocks are coupled through certain matrix elements and their conjugates: ⟨Φ∣ H{a† aa† bajai} ∣Φ⟩ = Γijab (3.15) The matrix elements of the two–body Hamiltonian are responsible for coupling npnh and (n ± 2)p(n ± 2)h states, which includes not only the coupling between ∣Φ⟩ and excitations, but the entire outermost side diagonals of the Hamiltonian matrix. This observation implies that we can manipulate the Hamiltonian to take the form depicted in the top right panel of Figure 3.1 by defining its offdiagonal part as Hod ≡ ∑ ai fai{a† aai} + 1 4 ∑ abij Γabij{a† aa† bajai} + H.c. (3.16) In the following section, we will demonstrate that the IMSRG flow exponentially reduces the matrix elements of Hod and successfully accomplishes the desired decoupling in the limit s → ∞. 59 3.4 Generators and Decay scales Several commonly used generators for the single-reference scenario have the ansatz [37, 41] η = ∑ ai ηai{a† aai} + 1 4 ∑ abij ηabij{a† aa† bajai} − H.c. (3.17) where H.c. denotes the Hermitian conjugate. The matrix elements of the generator are often either directly proportional to those of the off-diagonal Hamiltonian we introduced above, or simple functions of those matrix elements. The most commonly used IMSRG generators are: 1. The White generator, named after S.R. White whose work [107] motivated this ansatz: ηai = fai ∆ai , ηabij = fabij ∆abij , where the quantities ∆ai = ⟨Φa i ∣ H ∣Φa i ⟩ − ⟨Φ∣ H ∣Φ⟩ = fa − fi − Γaiai , (3.18) (3.19) ∆abij = ⟨Φab i j∣ H ∣Φab i j⟩ − ⟨Φ∣ H ∣Φ⟩ = fa + fb − fi − fj + Γabab + Γijij − Γaiai − Γajaj − Γbibi − Γbjbj . (3.20) are the differences of the Hamiltonian’s diagonal matrix elements for 1p1h or 2p2h states and the diagonal element for the reference state ∣Φ⟩ . They play the role of the energy differences that govern the decay scales for the matrix elements of f and Γ, analogous to the free-space SRG discussion in Section 2.2. 2. The White-ArcTan generator, also motivated by [107]. Since the matrix elements of 60 the White generator diverge when the energy differences (3.19) and/or (3.20) vanish, causing the IMSRG flow to stall, it can be useful to work with a “regularized” form of the White generator instead, which is given by ηai = 1 2 ∑ ai arctan 2fai ∆ai , ηabij = 1 2 ∑ abij arctan 2Γabij ∆abij . (3.21) Note that for vanishing energy differenes, the matrix elements approach π 4 . 3. The imaginary time generator [37, 39, 71], which is inspired by imaginary-time evo- lution techniques that are frequently used in Quantum Monte Carlo methods. It is defined by ηai = sgn(∆ai) fai , ηabij = sgn(∆abij) fabij . (3.22) 4. Finally, the Wegner generator, which is given as usual by η = {Hd, Hod} (3.23) with Hod defined in Eq. (3.16), and Hd ≡ H − Hod. Its one- and two-body parts have a similar structure as the IMSRG(2) flow equations (3.10) and (3.11), so we do not include it here, but refer to [37, 39, 41] for the full details. We note that this generator does not reduce to the form (3.17) in general — the Wegner generator can have additional nonzero matrix elements. One can perform a perturbative analysis along the lines of our discussion in Section 2.2 to determine the asymptotic decay scales for the off-diagonal Hamiltonian’s matrix elements 61 that result from each of the generators [39, 41]. It turns out that Γabij(s) ≈ Γabij(0) exp (−∣∆abij∣n s) , (3.24) where n = 0 for the White generator, n = 1 for the imaginary-time generator, and n = 2 for the Wegner generator. Note that this implies that s has different units in each case. 3.5 Magnus Formulation of IMSRG While the computational scaling of the IMSRG(2) approach is modest and comparable to other state-of-the-art many-body methods, we still may need a large number of integration steps to achieve “convergence” to the asymptotic limit of the flow for s ≫ ∞. It is also necessary to use high-order ODE solvers like fourth-fifth order Runge-Kutta (RK45) when integrating the IMSRG flow equations (3.9)–(3.11) in order to avoid an accumulation of integration-step errors, and such solvers require multiple copies of the solution vector, which in our cases can consist of hundred of millions of evolving matrix elements fpq(s), Γpqrs(s). This is a substantial computational bottleneck of the technique (refer to, for instance, [14, 42, 43]). Moreover, to calculate the expectation values of observables O in addition to the Hamil- tonian, we would in principle need to consistently evolve general operators using the flow equation dO ds = [η, O] , (3.25) with the same generator η[H(s)] as in Eq. (2.6). Since storing the generators at every point in the flow trajectory is simply too expensive, we are forced to solve the equation (3.25) 62 simultaneously with the flow equation for the Hamiltonian. This not only increases the dimension of the ODE system, but may cause additional stiffness, especially if the operators evolve at different characteristic scales than the Hamiltonian. These limitations can be overcome by applying Magnus expansion methods [8, 69, 71] to the the differential equation for the unitary transformation itself, dU (s) ds = −η(s)U (s) . (3.26) Mathematically, the solution to this ODE is given by the s- or path-ordered exponential : s s U (s) = Sexp ∫ 0 1 n! ∫ 0 ≡ ∑ n ds′η(s′) ds1 ∫ 0 s ds2⋯ ∫ s 0 dsnS{η(s1)⋯η(sn)} (3.27) (3.28) Due to the dynamical changes that occur in the generator during the evolution, the expres- sion (3.27) should be merely understood as a short-hand for the Dyson series (3.28). The s–ordering operator S guarantees that the integrands in question always feature the flow parameters in a descending order, whereby s1 > s2 > ⋯. However, this form of the solution is of limited practical value in calculations due to several reasons. First, it lacks specific instructions on how the series should be truncated. Second, we would again need to store η for multiple s–values to evaluate it. Finally, it is unclear how the Hamiltontian and other observables can be consistently transformed in a fully linked and size-extensive manner if we inevitably need to introduce truncations [71]. In 1954, Magnus proposed that under unspecified conditions the Dyson series can be 63 resummed into a proper exponential of an operator, U (s) = eΩ(s) (3.29) where the Magnus operator Ω(s)† = −Ω(s) is anti-Hermitian and satisfies Ω(0) = 0, and he developed formal techniques for the [69]. The review paper [8] discusses some of the formal conditions for the existence of Ω(s) as well as the convergence of the series — how- ever, these are difficult to check mathematically in our context. Assuming that the IMSRG transformation can be expressed as in Eq. (3.29), we can obtain a flow equation for Ω(s): dΩ(s) ds = ∞ ∑ 0 Bk k! adk Ω(η) (3.30) where Bk are Bernoulli numbers, and adk Ω(η) = [Ω, adk−1 Ω (η)], ad0 Ω(η) = η . (3.31) Similary to the standard IMSRG(2), both η and Ω as well as any of the commutators in the nested series on the RHS of Eq. (3.30) are truncated at the two–body level. This leads to the development of a new computational approach, referred to as the Magnus(2) formulation of the IMSRG. The nested commutators generated by adk Ω are evaluated recursively until satisfactory convergence of RHS of Eq. (3.30) is achieved [8, 18, 24]. At every integration step, the Baker–Campbell–Hausdorff (BCH) formula is utilized to construct the Hamiltonian and other observables, H(s) = eΩH(0)e−Ω = ∞ ∑ k=0 1 k! adk Ω(s)(H(0)) (3.32) 64 O(s) = eΩO(0)e−Ω = ∞ ∑ k=0 1 k! adk Ω(s)(O(0)) (3.33) The Magnus formulation offers a significant advantage in that the flow equations for Ω(s) can be solved using a straightforward first–order Euler step method without compro- mising accuracy. This leads to considerable memory savings and moderate reduction in CPU time. Although integration–step errors accumulate in Ω(s) with a first–order method, the transformation remains explicitly unitary upon exponentiation, and the transformed H(s) = U (s)H(0)U †(s) is unitarily equivalent to the initial Hamiltonian, except for the truncations made during the evaluation of BCH formula [71]. Maybe even more importantly, we can store the final Ω(∞), and use it to compute any observable of interest by simplify applying Eq. (3.33) at will. 65 Chapter 4 Factorization Techniques In this chapter, we will provide an overview of various decomposition techniques. Decom- position methods are essential tools that allow us to represent large and complex matrices as a product of smaller, more manageable matrices, vectors, or tensors. These decompo- sition techniques provide powerful tools that we hope to use to reduce the computational complexity of the SRG and IMSRG, and can lead to more efficient and accurate simulations of quantum many-body systems. In the following sections, we will delve into each of these techniques in more detail. 4.1 Singular Value Decomposition Consider an m × n matrix, A, with real entries, and let r be the smaller of the two di- mensions. For any such matrix, a singular value decomposition (SVD) exists in the following form: A m×n = U m×r Σ r×r †1 V r×n (4.1) where U and V are orthonormal matrices and Σ is a diagonal matrix. The columns (uj)r j=1 and (vj)r j=1 of U and V are the left and right singular vectors of A, respectively. The singular values of A are the diagonal entries (σj)r j=1, which we usually assume to be arranged in 1In mathematical formalism, ∗ is sometimes used to represent the conjugate transpose of a matrix. 66 descending order so that σ1 ≥ σ2 ≥ ⋯ ≥ σr ≥ 0. To put it another way, U = [u1u2⋯ur], V = [v1v2⋯vr], Σ = diag(σ1, σ2, ⋯, σr) (4.2) The factorization (4.1) represents A as a sum of r rank-one matrices A = r ∑ j=1 σjujv† j , (4.3) and we are primarily interested in the cases where the singular values σj decay rapidly to zero, resulting in the swift convergence of the sum. In such instances, we can approximate Ak ≈ A, where Ar is defined by the truncated sum Ak = k ∑ j=1 σjujv† j = UkΣkV † k (4.4) with k ≤ r. It was proven by Eckart and Young [20] that the truncated SVD gives the optimal, i.e., most precise, rank-k approximation to a matrix A. 4.2 Randomized SVD The classical algorithms for computing the SVD (4.1) can be computationally expen- sive for large matrices like the ones we work with in (nuclear) many-body theory since its asymptotic cost is O(mnr), where again r = min(m, n). Naively, we would first need to cal- culate the full SVD before applying a truncation to get the optimal low-rank approximation. However, randomized sampling can significantly speed up the algorithms for obtaining the low-rank SVD factorization, with only minor losses in accuracy [32]. 67 In essence, randomized factorization algorithms construct low-rank approximation of a matrix by applying the desired factorization to a smaller matrix derived from the original matrix with random projection techniques. In this section, we focus on constructing an approximate low-rank singular value decomposition using randomization. To obtain a sample of the column space of A, we compute Y = AΩ, where Ω is an n × (k + p) Gaussian random matrix, k is the desired rank, and p a small oversampling parameter. We then obtain an orthonormal column matrix Q of size m × (k + p) by performing a QR factorization on Y . If Y captures a good portion of the range of A and assuming that k is large enough relative to the numerical rank of A, we can approximate A with QQ†A. We can then compute the factorization of the smaller (k + p) × n product matrix Q†A and multiply the result by Q to obtain an approximation of A. The matrix product P = QQ† is a projector onto the range of Q and hence onto that of Y , which is obtained via a compact QR factorization of Y . The equality of QQ†A and A is achieved when Y captures the entire range of A. Mathematically, if R(Q)2 ≈ R(A), then QQ†A ≈ A [32]. We may prove this informally by splitting the SVD of A, A = [ U1 U2 ] ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Σ1 0 0 Σ2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ V † 1 V † 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (4.5) where the left singular matrix U is split into two parts U1 of size m × k and U2 of size m × (r − k), the singular value matrix Σ into two diagonal matrices Σ1 and Σ2 with size k × k and (r − k) × (r − k), respectively, and the transpose of the right singular matrix V † into 1 of size k × n and V † V † 2R designates the range, or column space, of a matrix. 2 of size (r − k) × n. Now suppose we take a Gaussian independent 68 identically distributed matrix Ω of size n × (k + p) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ V † 1 V † 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ V † 1 Ω V † 2 Ω Ω ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (4.6) Y = AΩ = [ U1 U2 ] = [ U1 U2 ] = [ U1 U2 ] ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Σ1 0 0 Σ2 Σ1 0 Σ1Ω1 0 Σ2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ Σ2Ω1 where we define Ω1 ≡ V † 1 Ω and Ω2 ≡ V † 2 Ω of size k × (k + p) and (r − k) × (k + p) respectively. Since R(Q) = R(Y )3, the projector onto the range of Y is equivalent to that of Q. Thus, we have PY = QQ† and the error of A and Ar is A − QQ†A = A − PY A = (1 − PY )A (4.7) In [32], it is demonstrated that if the matrix Ω1 has a full row rank, then the error of this approximation can be bounded by 2 ∥(I − PY )A∥ ≤ ∥Σ2∥2 + ∥Σ2Ω2Ω† 1∥ (4.8) where ∥⋅∥ is either the spectral or Frobenius norm. Obviously, if A has exactly rank k, then Σ2 = 0, causing the right hand side of the equation to be 0. Consequently, QQ†A = A, and Y will cover the range of A. As k gets closer to the rank of A, Y will cover more of R(A), and the product QQ†A will converge to A. Therefore, if we can find an orthonormal matrix 3This is a result of the orthogonal matrix Q being obtained via a compact QR factorization of Y . 69 Q such that ∥(I − QQ†)A∥ ≤ ϵ , (4.9) we can factorize Q†A of size (k +p)×n and multiply the result by Q to obtain an approximate factorization of A with the same approximation error as stated in (4.8) or in (4.9). If k is much smaller than min(m, n), then Q†A is considerably smaller than A. We can then compute the SVD of full rank k + p of Q†A and multiply the result by Q to obtain an approximate SVD of A. A ≈ Q(Q†A) = Q( ˜U DV †) = (Q ˜U )DV † (4.10) where D is the truncated singular value matrix of Q†A, and thus of A, and Q ˜U is the left singular value matrix of A. 4.3 Power Iteration We now introduce a power iteration method that enhances the accuracy of the low-rank approximation for matrices with large and slowly decaying singular values. To implement this method, instead of forming the matrix AΩ, we construct a matrix Y = ((AA†)qA)Ω, where q is an integer power greater than or equal to 1. We then obtain the orthogonal matrix Q from Y . It is important to note that both A and (AA†)qA share the same eigenvectors and corresponding eigenvalues. By substituting the SVD A = U ΣV †, we see that (AA†)qA = (U Σ2U †)qA = U Σ2qU †U ΣV † = U Σ2q+1V † . (4.11) 70 For matrices A with slowly decaying singular values, (AA†)qA with q ≥ 1, exhibits a notice- ably accelerated decay rate for its own singular values. 4.4 Algorithms for Low-Rank Approximation via Ran- domized SVD To compute a truncated, low–rank singular value decomposition of a matrix A of size m × n, the basic randomized algorithm involves the following steps: 1. Generate a random matrix Ω of size n × l with l = k + p. 2. Form a sample matrix Y of size m × l by multiplying Y = AΩ. 3. Apply power iteration if necessary. 4. Orthogonalize the sample matrix to obtain the matrix Q = orth(Y ). 5. Project the original matrix A into a lower-dimensional matrix B = Q†A, where B is an l × n matrix significantly smaller than A. 6. Compute the SVD of the smaller matrix B as B = ˜U ΣV †. 7. Form the matrix U by multiplying U = Q ˜U . 8. Approximate the rank–k SVD of A by setting Ur = U (∶, 1 ∶ k) , Σr = Σ(1 ∶ r, 1 ∶ k) , Vr = V (∶, 1 ∶ k) , so that A ≈ UkΣkV † k . 71 Algorithm 1 RSVD Algorithm Version I Input : A ∈ Rm×n, integer rank parameter k < min(m, n), an integer oversampling parameter p > 0, an integer power iteration parameter q ≥ 0, and an integer re– orthogonalization amount parameter s ≥ 1. Output : Matrices Uk ∈ Rm×k, Σk ∈ Rk×k, and Vk ∈ Rk×n. [Y, ⋅] = qr(Y, 0); if mod((2j − 2), s) == 0 then end if Z = A†Y ; if mod((2j − 1), s) == 0 then 1: Set l = k + p and initialize a matrix R ∈ Rn×l with Gaussian random entries; 2: Form samples matrix Y = AR and perform optional power iteration; 3: for j = 1 to q do 4: 5: 6: 7: 8: 9: 10: 11: 12: end for 13: Orthonormalize the columns of Y in [Q, ⋅] = qr(Y, 0); 14: Construct the smaller matrix B = Q†A derived from A; 15: Form the even smaller l × l matrix T = BB†; 16: Perform eigendecomposition of l × l matrix [U, D] = eig(T ); 17: Form the approximate low-rank SVD components of A using the results of the eigende- end if Y = AZ [Z, ⋅] = qr(Z, 0); √ composition Σk = D, Uk = QU, Vk = B†U Σ−1 k ; 18: Extract components corresponding to the k largest by magnitude singular values Uk = Uk(∶, (p + 1) ∶ l); Σk = Σk((p + 1) ∶ l, (p + 1) ∶ l); Vk = Vk(∶, (p + 1) ∶ l); Next, we present alternative randomized algorithms for approximating low–rank SVD adapted from [98]. Algorithm 1 introduces a modification that obviates the necessity to perform the SVD computation on the substantial l × n matrix B. This size can persist when the number of columns within matrix A is extensive. Instead, the algorithm leverages the eigendecompo- sition of the more compact k × k matrix BB†. In contrast, Algorithm 2 takes a diffrent approach by employing the QR factorization of B† to generate a reduced matrix R, which is then used for the subsequent SVD calculation. It’s important to highlight that even the downscaled matrix B, which corresponds to the 72 Algorithm 2 RSVD Algorithm Version II Input : A ∈ Rm×n, integer rank parameter k < min(m, n), an integer oversampling parameter p > 0, an integer power sampling parameter q ≥ 0, and an integer re– orthogonalization amount parameter s ≥ 1. Output : Matrices Uk ∈ Rm×k, Σk ∈ Rk×k, and Vk ∈ Rk×n. [Y, ⋅] = qr(Y, 0); if mod((2j − 2), s) == 0 then end if Z = A†Y ; if mod((2j − 1), s) == 0 then 1: Set l = k + p and initialize a matrix R ∈ Rn×l with Gaussian random entries; 2: Form samples matrix Y = AR and utilize optional power scheme; 3: for j = 1 to q do 4: 5: 6: 7: 8: 9: 10: 11: 12: end for 13: Orthonormalize the columns of Y in [Q, ⋅] = qr(Y, 0); 14: Compute the smaller matrix Bt = A†Q; 15: Obtain the small l × l matrix R using a compact QR factorization [Q, R] = qr(Bt, 0); 16: Take the SVD of the l × l matrix R, [U, Σk, V ] = svd(R); 17: Form the approximate low rank SVD components of A using the results of the SVD of end if Y = AZ [Z, ⋅] = qr(Z, 0); R. Uk = QV , Vk = QU ; 18: Extract component corresponding to the k largest by magnitude singular values Uk = Uk(∶, 1 ∶ k) Σk = Σk(1 ∶ k, 1 ∶ k); Vk = Vk(∶, 1 ∶ k); larger matrix A, might still pose challenges for multiplication due to its size. In scenarios where this hold true, a customized adaptation of Algorithm 1 prove useful. This involves assessing the matrix-matrix product BB† column-wise, utilizing standard basis vectors like BB†ej, particularly when matrix B becomes excessively large. This methodology can also be extended to the computation of Vk through the application of the matrix-matrix product B†U Σ−1. 73 Part II Factorization in Few-Body Systems 74 Chapter 5 SVD-SRG evolution in Two-Nucleon Interactions Ab initio nuclear many–body theory faces challenges in the coming decade due to the increase in computational and storage costs required for calculations of heavy, neutron- rich, and structurally complex nuclei. This chapter1 explores the factorization of nuclear interactions as a first step towards solving this issue. SVDs of nucleon–nucleon interactions are performed in partial wave representation, and the dependence of the singular value spectrum on interaction characteristics is studied. Next, an SRG evolution is developed and implemented directly in terms of the relevant singular vectors. Results show that low– resolution interactions allow for the truncation of the SVD at low rank, and a small number of relevant components are sufficient to capture the nuclear interaction and perform an accurate SRG evolution. Importantly, the rank is uniform across all partial waves and almost independent of the basis choice. This chapter also discusses strategies for mitigating the growth of the rank that is caused by the subsequent transformation of the interaction between the center–of–mass and laboratory frames. The low rank approximation to the SRG–evolved interactions is tested in many–body calculations using the in–medium SRG, and the implementation of the SRG using the singular vectors of the interaction is shown 1This chapter is adapted from my publication [113]. 75 to not affect the evolution of other observables, which is verified by analyzing nuclear mean- square radii. 5.1 Free-Space SRG Throughout this chapter, we will be studying SRG evolutions in a momentum-space partial-wave representation, as introduced in chapter 2. We will use the standard generator for performing momentum decoupling (also see [10, 35, 96]) η(s) = [T, H(s)] (5.1) where T is the relative (or intrinsic) kinetic energy — we drop the label for the sake of brevity. We keep the kinetic energy constant during the flow, hence all s-dependent contributions from evolving T are absorbed into the interaction V (s): H(s) = T + δT (s) + ¯V (s) ≡ T + V (s) with δT (0) = 0, ¯V (0) = V (0) (5.2) (5.3) This partitioning will be important later. Substituting Eq. (5.1) and Eq. (5.2) into Eq. (2.6) and setting dT ds = 0 yields the following equation for V (s) (cf. Eq. (2.24)): d ds V (s) = [η(s), T ] + [η(s), V (s)] = [[T, V (s)] , T ] + [[T, V (s)] , V (s)] . (5.4) 76 This evolution can be conveniently implemented in momentum space, where T will be rep- resented by a diagonal matrix. As discussed in previous chapters, general observables can in principle be calculated by co-evolving them with the Hamiltonian according to d ds O(s) = [η(s), O(s)] , (5.5) but each additional observable increases the size of the ODE system by the full dimension of the momentum basis and may affect its overall stiffness. Because we are dealing with matrix-based SRG evolutions here, several methods exist to address this issue. One such method is to obtain the eigenvalue solutions of the initial and final Hamiltonian matrices and construct U (s) in Eq. (2.1) directly as [2, 26, 88]. U (s) = ∑ n ∣ψn(s)⟩ ⟨ψn(0)∣ . (5.6) Another way to obtain U (s) is by directly solving the ODE d ds U (s) = η(s)U (s), U (s = 0) = I , (5.7) and a third option is to employ the Magnus expansion (see Section 3.5), but this approach was found to be susceptible to convergence problems [96]. 77 5.2 SRG evolution of SVD factors Let us assume that we have successfully performed an SVD decomposition of the initial Hamiltonian H(0), H(0) = ∑ i ∣ui(0)⟩ σi ⟨vi(0)∣ , and consider its SRG evolution. The evolved Hamiltonian can be written as H(s) = ∑ i U (s) ∣ui(0)⟩ σi ⟨vi(0)∣ U †(s) ∣ui(s)⟩ σi ⟨vi(s)∣ , ≡ ∑ i (5.8) (5.9) where we have used the invariance of singular values under unitary transformations and defined evolved singular vectors as ∣ui(s)⟩ and ∣vi(s)⟩. Inserting Eq. (5.9) into the flow equation (2.6), d ds ∑ i (∣ui(s)⟩ σi ⟨vi(s)∣) = [η(s), H(s)] (η(s) ∣ui(s)⟩ σi ⟨vi(s)∣ − ∣ui(s)⟩ σi ⟨vi(s)∣ η(s)) . = ∑ i The LHS of this equation can be rewritten using (5.10) d ds (∣ui(s)⟩ σi ⟨vi(s)∣) = ( d ds ∣ui(s)⟩) σi ⟨vi(s)∣ + ∣ui(s)⟩ σi ( d ds ⟨vi(s)∣) . (5.11) By comparing Eq. (5.10) and Eq. (5.11), we can then easily derive flow equations for the singular vectors: d ds ∣ui(s)⟩ = η(s) ∣ui(s)⟩ , (5.12a) 78 Figure 5.1: Singular value spectra of H and V in the proton–proton (solid curves) and neutron–proton 1S0 partial waves (dashed curves), for a chiral N3LO two–nucleon interaction (EM, cutoff Λ = 500 MeV/c, see text and Ref. [22]). The units of the singular vectors result from adopting scattering units ̵hc/mc2 = 1 as well as momentum–space discretization discussed in Sec. 2.3. d ds ∣vi(s)⟩ = η(s) ∣vi(s)⟩ . (5.12b) When attempting to apply the SRG flow to the SVD of the Hamiltonian, a significant obstacle arises due to the unbounded nature of the kinetic energy term. This prevents us from identifying a natural low-rank truncation for the components of H(0), as illustrated in Figure 5.1. Fortunately, Eq. (5.4) suggests that we should consider the SVD of the evolving part of the Hamiltonian anyway, which is V (s) as defined above. At s = 0, the SVD of the interaction is given by V (0) = ∑ ij ∣ui(0)⟩ σij ⟨vj(0)∣ , σij = σiδij . (5.13) 79 020406080100120140i107105103101101i [fm2]1S0HppHnpVppVnp As an example, Figure 5.1 also displays the singular value spectrum of the proton–proton (pp) and neutron-proton (np) 1S0 partial waves of a realistic chiral N3LO interaction [22]. This spectrum shows that the interaction possesses a low–rank structure. We note that the inclusion of the Coulomb interaction can increase the rank of the interaction in the pp channels, but not to an extent where truncations would be impractical. We also note that the singular values of the interaction decay exponentially, while the kinetic term in the Hamiltonian grows only quadratically. As a result, the singular value spectrum of the generator (5.1) will also decay exponentially. Therefore, the SRG evolution via Eq. (5.4) will not alter the low–rank structure of the interaction, except for potential truncation artifacts. For s > 0, we have V (s) = ∑ ij ∣ui(s)⟩ (σij + δσij(s)) ⟨vj(s)∣ (5.14) where the set of singular values pertaining to the initial interaction is invariant during unitary evolution but a δσij(s) arises from the incorporation of induced contributions from the evolution of the kinetic energy. By examining the expressions on both sides of Eq. (5.4), we can formulate a flow equation for δσij(s). Substituting Eq. (5.14) into the equation, we obtain d ds V (s) = [η(s), V (s)] + ∑ ij ∣ui(s)⟩ ( d ds δσij(s)) ⟨vj(s)∣ (5.15) Using the orthogonality of the singular vectors, we obtain d ds δσij(s) = ⟨ui(s)∣ [η(s), T ] ∣vj(s)⟩ (5.16) In principle, it is necessary to solve this flow equation alongside Eq. (5.12), but it turns out that it can be solved analytically. By expanding the right–hand side and substituting 80 Eq. (5.12) to convert to derivatives of the singular vectors, we can obtain a closed solution: d ds δσij(s) = ⟨ui(s)∣ [η(s), T ] ∣vj(s)⟩ = ⟨ui(s)∣ (η(s)T − T η(s)) ∣vj(s)⟩ = ⟨ui(s)∣ η(s)T ∣vj(s)⟩ − ⟨ui(s)∣ T η(s) ∣vj(s)⟩ = − ⟨ui(s)∣ η†(s)T ∣vj(s)⟩ − ⟨ui(s)∣ T η(s) ∣vj(s)⟩ = − ( d ds ⟨ui(s)∣) T ∣vj(s)⟩ − ⟨ui(s)∣ T ( d ds ∣vj(s)⟩) (5.17) where we recognize that the last line is merely the derivative of ⟨ui(s)∣ T ∣vj(s)⟩. Integrating both sides, we obtain the closed solution δσij(s) = ⟨ui(0)∣ T ∣vj(0)⟩ − ⟨ui(s)∣ T ∣vj(s)⟩ , (5.18) i.e., it is necessary to calculate the matrix elements of the kinetic energy within the repre- sentation defined by the initial and final singular vectors. However, we can rewrite the first term in Eq. (5.18) as follows: ⟨ui(0)∣ T ∣vi(0)⟩ = ⟨ui(0)∣ U †(s)U (s)T U †(s)U (s) ∣vi(0)⟩ = ⟨ui(s)∣ U (s)T U †(s) ∣vi(s)⟩ = ⟨ui(s)∣ T + δT (s) ∣vj(s)⟩ . (5.19) Plugging this back into Eq. (5.18), the contributions from T cancel and we see that δσij(s) = ⟨ui(s)∣ δT (s) ∣vj(s)⟩ , (5.20) 81 i.e., we only need to evaluate the s–dependent portion of the kinetic energy within the (truncated) basis of evolved singular vectors. We remark that while the evolved SVD representation of V (s) will have a non-diagonal singular value matrix σ + δσ(s) in the end, we can readily diagonalize it to obtain a “final” set of singular vectors if we so desire, since its dimensions match the rank of the truncated SVD and are therefore small. The construction SRG transformation U (s) can be readily constructed from the initial and final singular vectors ∣ui(s)⟩ and ∣vi(s)⟩: U (L)(s) = ∑ i U (R)(s) = ∑ i ∣ui(s)⟩ ⟨ui(0)∣ , ∣vi(s)⟩ ⟨vi(0)∣ . (5.21a) (5.21b) In principle, we need to differentiate between left and right unitary transformations, depend- ing on whether we construct U (s) from the ∣ui(s)⟩ or ∣vi(s)⟩. Since we are usually interested in the evolution of Hermitian operators, the difference between the left and right singular vectors is merely an irrelevant phase factor. By means of the U (s) constructed from the singular vectors, it is feasible to transform any arbitrary observables O according to O(s) = U (s)O(0)U †(s) . (5.22) 5.3 Applications in the Two-Nucleon System We are now ready to study the interplay of the SVD and SRG in the two-nucleon system. We will first perform SVDs of interactions that have been evolved using the SRG before 82 applying the SVD-based SRG formalism developed in the previous section. 5.3.1 Momentum Discretization and Conventions Throughout this section, we will work with the discretized momentum-space partial-wave representation introduced in Section 2.3. In particular, we absorb the weights and q2 i factors from the integration measure into the interaction matrices (see Eq. (2.29)), so that they do not appear explicitly in Eq. (5.12). This has the added benefit that it simplifies the comparison of truncated SVDs among different mesh and basis choices, and it may facilitate the interpretation of the singular components as mere representations of particular basis operators for the interaction in the future. When comparing with the existing literature on the SRG transformation of matrix elements (e.g., [10, 35, 96]), we will simply revert back to the original basis. Figure 5.2 presents our findings regarding the independence of the singular values on the selected mesh for the proton–proton (pp) and neutron–proton (np) partial waves of the N N interaction, using the chiral N3LO interaction with Λ = 500 MeV/c by Entem and Machleidt as an example. We refer to this interaction as the EM interaction in the subsequent discussion. We found that the singular value spectra for two equidistant meshes with the same maximum momentum but different spacings are practically identical for σi ≥ 10−10 fm−2. The precision with which the input matrix elements are stored affects the singular value spectrum’s extension. Specifically, reducing the number of decimal digits from 10 to 6 leads to an elongated tail of unphysical singular values of the order of σi ≈ 10−6 fm−2. These contributions to the interaction would not be considered in practical applications anyway — based on the computed spectra, a natural threshold σt ≈ 10−3 fm−2 emerges that should ensure high accuracy in the pp and np 1S0 partial waves. This threshold is discussed further 83 Figure 5.2: Singular value spectra of V in the proton-proton and neutron-proton 1S0 partial waves, for the EM interaction. The solid and dashed curves are obtained on equidistant meshes with the same qmax = 7 fm−1, but different spacings ∆q. The dotted curve illustrates that the tail of the singular value spectrum is related to the precision of the input data, which is reduced from 10 digits for the other data sets to 6 digits after the decimal. All conventions are the same as in Fig. 5.1. in next subsection. 5.3.2 Singular Value Spectra of Nucleon-Nucleon Interaction With the technical discussion out of the way, we can proceed to analyze the SVDs of N N interactions, such as the EM [22] and AV18 [110] potentials. As mentioned before, these interactions are accurate representations of nucleon-nucleon scattering data, but they have different characteristics. The EM potential is derived using chiral Effective Field Theory (EFT) and has a moderate initial cutoff. The interaction is non–local due to the adopted regularization scheme and is best represented in momentum space. On the other hand, the AV18 potential is designed for use in coordinate–space Quantum Monte Carlo (QMC) 84 020406080100120140i109107105103101101i [fm2]1S0pp, q=0.10 fmpp, q=0.05 fmpp, q=0.05 fm, 6 digitnp, q=0.10 fmnp, q=0.05 fmnp, q=0.05 fm, 6 digit calculations and is constructed to be as local as possible. Consequently, it has a strong, repulsive core in the interaction, which makes its use problematic in many–body methods that rely on basis expansions. Although newer versions of these interactions have been developed, EM and AV18 are still relevant for analysis because they contain all the essential features found in more recently developed interactions, as indicated in recent reviews [23, 80, 83]. As shown in Figures 5.1 and 5.2, the 1S0 partial waves of the EM interaction exhibit a low rank structure, as evident from the singular value spectra. In the proton-proton (pp) channel, a ”kink” is observed, indicating a natural threshold value of σt = 10−3 fm−2 for truncating the SVD. This threshold value is approximately three to four orders of magnitude smaller than the largest singular values in the pp and neutron-neutron (np) channels for all partial waves up to j = 9, which are σmax(pp) ≈ σmax(np) ≈ 3 fm−2. The number of components of the interaction in the pp channel that have singular values above the threshold is approximately 35, while in the np channel, only about 10 components are above the threshold. The reason for this difference is due to the inclusion of the Coulomb interaction in the initial matrix elements, which is long–ranged and lacks an inherent scale. In the absence of a cutoff or regulator, the matrix representation of (VC)ij ∼ (qi − q′ j)−2 would diverge on the diagonal, making it impossible to truncate the SVD in the pp channels. In this study, a hard cutoff is imposed at RC = 15 fm, which serves two purposes: in phase–shift calculations, it provides the matching scale between the momentum–space wave function obtained from the matrix elements of the nuclear and cutoff Coulomb potentials VN + VC and the asymptotic Coulomb wave function [65, 97]. In nuclear structure calculations, the impact of imposing this cutoff on Hartree–Fock and IM–SRG ground–state energies is below 1 keV. 85 Figure 5.3: Truncated SVD rank r of partial waves up to j = 9 for the EM interac- tion, using a singular value threshold σt = 10−3 fm−1. Points for partial waves with dif- ferent l and s can overlap. The maximal singular values for the interaction channels are{σmax(nn), σmax(np), σmax(pp)} = {3.00 fm−1, 2.97 fm−1, 3.06 fm−1}, respectively. It is evident that changing the value of RC will affect the number of components that exceed a pre-selected singular value threshold. If we are willing to accept a small impact on ground–state energies, which remains negligible compared to other sources of uncertainty, we can reduce the cutoff to RC ≈ 10 fm. Doing so would lower the rank of the truncated interactions studied to approximately 20 in the pp partial waves. However, any further decrease in RC would cause a rapid deterioration in the energies and phase shifts. The comparatively simple nature of VC as compared to nuclear interactions could facil- itate its efficient inclusion through means other than a truncated SVD in a chosen configu- ration space. However, in the context of the present study, it is important to note that VC cannot be treated completely separately as it couples with the other terms in the Hamilto- nian during the SRG evolution, even though it evolves weakly. One could try to create an SRG evolution from T and the nuclear interaction VN and apply it to VC subsequently, as 86 02468j05101520253035rnnnppp (a) (b) Figure 5.4: Singular values of the neutron-proton 1SS partial wave for the SRG-evolved EM (panel a) and AV18 interactions (panel b) at different resolution scales λ. with other observables. However, when we attempted this, it resulted in effects on the order of several percent on the ground–state energies of medium–mass nuclei, which is comparable to or greater than other theoretical uncertainties for these quantities. This topic warrants further detailed exploration in the future. As we move to higher partial waves, we see a consistent pattern across the different channels as shown in Figure 5.3. The number of singular values above the threshold σt = 10−3 fm−2 is between 5 and 10 for np and nn partial waves up to j = 9. On the other hand, pp partial waves contain about 35 singular values above the threshold, which is due to the presence of the Coulomb interaction. However, for off–diagonal pp partial waves like 3P2 − 3F2, the number of components is comparable to that in the np and nn channels. This is because the Coulomb interaction cannot contribute to partial waves with l ≠ l′, which further supports that increased number of relevant components in the l = l′ partial waves is due to the presence of VC. The relative number of components in different partial waves is similar across different channels. Next, we consider the impact of a free–space SRG transformation on the singular value distribution. Figure 5.4 shows the singular values in the neutron-proton 1S0 partial wave for 87 0102030405060i108106104102100i [fm2]EM1S0(np)10.0 fm14.0 fm13.0 fm12.5 fm12.0 fm11.8 fm11.5 fm11.2 fm1050100150200i108106104102100102i [fm2]AV181S0(np)10.0 fm14.0 fm13.0 fm12.5 fm12.0 fm11.8 fm11.5 fm11.2 fm1 the EM and AV18 interactions after they have been evolved to different resolution scales λ. For the EM interaction, we observe that the number of singular values beyond the threshold σt does not change substantially as we evolve from the initial interaction to λ = 1.8 fm−1. However, if we continue to evolve the interaction to λ = 1.5 fm−1 and λ = 1.2 fm−1, the growth of the rank accelerates, and its value becomes twice the value of the original interaction. At low resolution scales, the evolution starts to affect the momentum transfers (cf. Eq. (2.26)) that arise from the dynamics of one–pion exchange. In a projective renormalization group (RG) approach, we would say that we are “integrating out” the pion [10, 11]. However, since the SRG is unitary, it cannot eliminate interaction strength, but only redistribute it. Within the two–body sector, the SRG increasingly diagonalizes the V (s) by redistributing interaction strength toward the diagonal. As the band’s width becomes narrow enough, the strength begins to move along the diagonal to higher momenta [10], and previously insignificant components of the flowing interaction may now become relevant. It is crucial to consider the observed growth of the rank when we perform the SRG evolution directly in terms of the SVD factors, as explained in Section 5.3.1. Typically, we evolve the interaction to resolution scales around λ ≈ 2.0 fm−1, which has demonstrated to be an advantageous value for nuclear many–body calculations, thereby avoiding a substantial increase in rank. However, based on our present findings, it is recommended to “over– sample” and include additional components in the process to capture the renormalization group (RG) flow accuracy and prevent any loss of unitary in the two–body system. In Figure 5.4b, we show the singular value spectrum of the AV18 interaction. The unevolved AV18 has a higher rank than the EM interaction because it has a greater extension in momentum space (also see Fig. 1.1). This is due to the AV18 interaction’s hard core, which allows it to easily connect incoming the outgoing momenta that couple by up to 88 Figure 5.5: Detail view of the dominant singular values in selected S waves of the EM and AV18 interactions at different SRG resolution scales λ. 20 fm−1. When the AV18 interaction undergoes an SRG evolution, there is a significant change: the dominant singular values become smaller, but the rank rapidly expands. As a result, the interaction’s spectrum becomes so flat that possible truncation points for the SVD are between 150 and 200 components, which are substantial portions of the full rank of the matrix. Figure 5.5 presents a detailed view of the dominant singular values in the neutron–proton S waves of our interactions. In particular, for the EM interaction, we observe that only 2–3 singular values are highly dominant before the spectrum drops off rapidly. As we decrease λ, the relative dominance of one of these singular values increases compared to the others, before the growth of the rank becomes problematic beyond λ = 1.8 fm−1. This finding is consistent with the analysis of Bogner et al. [12], who showed that low-rank separable approximations 89 102101100101i [fm2]1S0(np)EM1S0(np)EM1S0(np)EM1S0(np)EM1S0(np)EM1S0(np)EM1S0(np)EM1S0(np)EM1S0(np)EM3S13S13S13S13S13S13S13S13S1051015i102101100101102i [fm2]1S0(np)AV181S0(np)AV181S0(np)AV181S0(np)AV181S0(np)AV181S0(np)AV181S0(np)AV181S0(np)AV181S0(np)AV18051015i3S13S13S13S13S13S13S13S13S110.0 fm14.0 fm13.0 fm12.5 fm12.0 fm11.8 fm11.5 fm11.2 fm1 to the N N interaction become more accurate as the resolution of the interaction is reduced. This is because the SVD can be considered as a generalization of such techniques. Interestingly, the analysis presented in Ref. [12] arrived at a similar conclusion for a low-resolution Vlow−k interaction that was obtained from AV18 through a projective RG decimation [10, 11]. Although only a few singular values dominate the spectrum at the beginning of our SRG evolution, the rapid flattening of the spectrum and the concurrent expansion of the rank seem to contradict the results obtained for Vlow−k. To address this discrepancy, we emphasize that the magnitude of a singular value is not always enough to determine whether a component of the interaction is significant for the physics we aim to explain or not. For example, suppose we have a local interaction with matrix elements in momentum space, which we can write schematically as Vreg(q, q′) = V (q, q′)f ( q − q′ λ ) (5.23) where f is a local regulator. In momentum space, the interaction between particles can be represented by a matrix that has a band diagonal structure. The parameter λ, which is named suggestively, determines the width of the band. However, both the strength of f and λ do not impose any limit on the total momentum of the particles (represented by q + q′). The AV18 potential is an example of a local potential with a strong repulsive core, and it has significant positive matrix elements at high momenta (see Figure 1.1 or Fig. 2 in Ref. [10]). The eigenvalues of the AV18 potential are dominated by the high–momentum region and are larger in magnitude than the negative eigenvalues from the attractive region at low momentum. The singular value spectrum reflects only the absolute value of the two types of eigenvalues, which get mixed together. 90 Figure 5.6: Matrix representation of the unitary transformation U (λ1, λ2) = ⟨uj(λ1)∣ ∣ui(λ2)⟩ between the initial and final resolution scales λ1, λ2 ∈ {2 fm−1, 4 fm−1, ∞}. Here, we specifi- cally focus on the 1S0(np) partial wave. The basis vectors are ordered by decreasing singular value. Note the logarithmic color scale. Additional information that may aid in determining the relevance of an interaction com- ponent is present within the structure of the associated singular vectors. As our primary focus lies in the behavior of these vectors under component–wise SRG evolution, we exam- ine the expansion of one set of singular vectors, denoted as ∣ui(λ)⟩ (or ∣vi(λ)⟩, in terms of another set at different λ value. The coefficients of this expansion correspond to the entries of the unitary evolution matrix (5.12), but with the constraint that neither of the scales is set to λ = ∞ (or, equivalently, s = 0 fm4). Figure 5.6 displays matrices for different values of (λ1, λ2), namely (∞, 4.0 fm−1), (∞, 2.0 fm−1) and (4.0 fm−1, 2.0 fm−1). Considering the EM interaction first, we note that the matrices 91 04080120050100juj()|ui(4.0)EM04080120uj()|ui(2.0)04080120uj(4.0)|ui(2.0)050100150i050100150200jAV18050100150i050100150200i1001011021030103102101100 appear similar , containing a washed out diagonal band in the upper left corner, a large cen- tral block, and a sharp diagonal in the lower right corner. It is relatively simple to connect the matrices shown in Figure 5.6 to our findings for the truncated SVD of the interaction. As depicted in Figure 5.4, the upper left and central blocks of the matrices contain the ∣ui⟩ of the singular values σi ≳ 10−6 and σi ≲ 10−6, respectively. The separation point between the blocks is likely determined by the precision settings of the ODE solver, as revealed in our analysis. Additionally, the lower right block of the matrices comprises the singular vectors that remain constant, indicating that V = 0 and H = T . The block structure is the result of the nonlocal regularization of the interaction. This regularization suppresses the momentum space matrix elements independently for incoming and outgoing momenta, as indicated by the equation Vreg = V (q, q′) exp (− ( 2n ) ) exp (− ( q Λ 2n ) ) q′ Λ (5.24) where n = 2 or n = 3 and Λ = 2.5 fm−1 [22]. Due to the matrices’ structure, the first r singular vectors remain nearly unconnected from the rest of the spectrum as we evolve. However, it should be noted that the central block extends towards lower indices in the U (∞, 2.0) matrix, indicating a slight increase in rank, as also noted in our discussion of Fig. 5.4 and 5.5. This phenomenon is less apparent in U (4.0, 2.0) because the two evolved bases are more alike. In the case of AV18, the structure of the unitary transformations is considerably more intricate. The constituents that will have the greatest impact on the low–momentum region are distributed throughout the basis, making it challenging to identify ahead of the evolution. As we evolve, the U (∞, λ) matrices become less organized, which corresponds to the increase in rank observed in Figure 5.4. The matrices indicate that we can expect to use around 160 92 Figure 5.7: Momentum-space matrix elements of the EM interaction at λ = 2.0 fm−1 in selected uncoupled partial waves. The top panels are by SVD-SRG, the bottom panels by conventional matrix-based SRG evolution. The matrix elements are given in scattering units (̵h = c = ̵h2/m = 1). Not the logarithmic color scale. to 170 out of 200 interaction components to apply the SRG using the interaction’s factorized form. However, there is a simplification in the structure of the U (4.0, 2.0) matrix. If we perform the SVD at λ = 4.0fm−1, we may be able to enhance the low–rank structure by safely removing the high–momentum components that are already decoupled, similar to what is done in Vlow−k. In the next subsection, we will see further evidence of this when we discuss the deuteron ground-state energy for AV18 (cf. Figure 5.12). The main conclusion from our analysis is that while the magnitude of the singular values is an important factor to consider when determining the necessary singular vectors for an effective low–rank representation of interaction in the low–momentum regime, it alone is not sufficient. In addition, other criteria that account for the momentum structure of the singular vectors may need to be taken into consideration. Specifically, when dealing with chiral interactions that employ nonlocal regulators, the selection process based solely on the singular values can be effective if the initial cutoff is not set too high, as there are generally not many strongly positive eigenvalues present in this scenario. The situation is more complex when dealing with local interactions, as the potential’s 93 02.04.06.02.04.06.02.04.06.02.04.06.02.04.06.02.04.06.0q[fm1]02.04.06.02.04.06.0q0[fm1]3S13S13S13D13D13D13S13D13S13D13S13D1r=15r=30exact1001011021031040.0104103102101100V(q,q0)[fm] Figure 5.8: Momentum-space matrix elements of the EM interaction at λ = 2.0 fm−1 in selected uncoupled partial waves. The top panels are obtained with SVD-SRG, the bottom panels with conventional matrix-based SRG evolution. The matrix elements are given in scattering units (̵h = c = ̵h2/m = 1). Note the logarithmic color scale. behavior and regularization at both short and long distances impact the interaction’s rank. Although the Coulomb interaction is singular in S waves as the distance between particles approaches zero, its rank is predominantly influenced by the regularization of its long–range behavior, thus allowing for an accurate low–rank expansion. In contrast, the AV18 inter- action cannot undergo significant rank reduction because it strong short-range coordinate- space potentials that parameterize (multi)meson exchange interactions, and we are sensitive to their effects because of AV18’s high resolution scale. Comparable issues are likely to arise with local or semi-local chiral interactions, as noted in previous studies [23, 80], although their lower resolution scale could help mitigate these concerns. 5.3.3 SVD-Based SRG Evolution In this section, we implement the SRG evolution of an SVD–factorized interaction, ab- breviated as SVD–SRG, using the framework developed in Section 5.2. As an example, we 94 02.04.06.02.04.06.02.04.06.02.04.06.0q[fm1]02.04.06.02.04.06.0q0[fm1]1S0(nn)1S0(np)1S0(pp)3P1(np)1S0(nn)1S0(np)1S0(pp)3P1(np)r=30r=30r=40r=30exactexactexactexact1.01011021031040.01041031021011.0V(q,q0)[fm] perform the SVD–SRG evolution of the EM interaction at λ = 2.0 fm−1. Figure 5.7 depicts the momentum space matrix elements in the deuteron channel for various SVD ranks. Using only 15 components per partial wave, some distortion is still observed. However, with 30 components, the SVD-SRG evolution reproduces the results from evolving with the complete matrix accurately, with absolute deviations of roughly 10−4 that are only visible upon close examination. It is worth noting that we did not attempt to fine–tune the rank r using the information from the previous section (e.g., Figure 5.3). Instead, the decision to employ the same rank for all partial waves was solely made for the sake of convenience. In Figure 5.8, we investigate the efficacy of the SVD–SRG approach in other selected partial waves. The matrix elements for the neutron–neutron and neutron–proton partial waves were obtained using 30 components, while 40 components were necessary for the 1S0 proton–proton partial wave due to the presence of the Coulomb interaction, as evidenced by the Coulomb tail along the diagonal, which is absent in other isospin channels. Weak oscillations and “fraying” can be observed around the edges of the primary structures, but the deviations’ absolute values are again at or below 10−4. Given the efficacy of the SVD–SRG on the matrix element level, we have utilized the factorized interactions to calculate nucleon–nucleon observables, such as scattering phase shifts and the binding energy of the deuteron. In Figure 5.9, we present the np phase shifts and mixing parameter of the SVD–SRG evolved EM interaction in the deuteron channel, as well as selected partial waves. The results indicate that about 5 to 10 interaction components are sufficient to achieve agreement with conventional matrix–based evolution, with only the 3S1 − 3D1 mixing angle ϵ1 requiring a few additional components. These findings align with our expectations based on the SVD of the initial interaction (cf. Figure 5.3), and the need to accommodate a slight increase in rank via “over-sampling” as we evolve to λ = 2.0 fm−1. 95 Figure 5.9: Selected neutron–proton phase shifts and mixing angles of the EM interaction resolution λ = 2.0 fm−1. The SVD–SRG evolution for different ranks is compared to the matrix–based evolution, which exactly preserves the phase shifts of the unevolved EM inter- action by construction. 96 0100200300400Elab [MeV]200204060 [deg]1S0Fullr=5r=10r=15r=200100200300400Elab [MeV]020406080100120140 [deg]3S1Fullr=5r=10r=15r=200100200300400Elab [MeV]20151050 [deg]3D1Fullr=5r=10r=15r=200100200300400Elab [MeV]101234561 [deg]3S13D1Fullr=5r=10r=15r=200100200300400Elab [MeV]2520151050510 [deg]3P0Fullr=5r=10r=15r=200100200300400Elab [MeV]25201510505 [deg]3P1Fullr=5r=10r=15r=200100200300400Elab [MeV]25201510505 [deg]3P1Fullr=5r=10r=15r=200100200300400Elab [MeV]25201510505 [deg]3P1Fullr=5r=10r=15r=20 Figure 5.10: Selected proton-proton phase shifts and mixing-angles of the EM interaction resolution λ = 2.0 fm−1. Note the higher ranks for the SVD-SRG evolution compared to Figure 5.9. 97 0100200300400Elab [MeV]200204060 [deg]1S0Fullr=25r=30r=35r=400100200300400Elab [MeV]0.02.55.07.510.012.515.017.520.0 [deg]3P2Fullr=25r=30r=35r=400100200300400Elab [MeV]0.00.51.01.52.0 [deg]3F2Fullr=25r=30r=35r=400100200300400Elab [MeV]0.00.51.01.52.02.53.02 [deg]3P23F2Fullr=25r=30r=35r=40 When dealing with the Coulomb interaction in the proton–proton channel, we have found it necessary to increase the number of components to 30−35, as indicated by the partial wave and mixing angle results in Figure 5.10. Initially, the convergence of the phase shift with the rank r may appear irregular. However, this is a result of gradually including components based on the size of the associated singular value alone, regardless of whether they are dominated by the Coulomb or the nuclear interaction. As we increase the rank from 25 to 30, we add components that are primarily influenced by VC, which significantly affects the phase shifts at low energy or long distances. While this effect is somewhat obscured by the overall size of δ in the 1S0 phase shift plot, it is quit clear in higher partial waves. Increasing the rank from 30 − 35 adds components that are dominated by nuclear interactions, which have a substantial impact on the Elab = 100 − 250 MeV region, and for r > 35, the phase shift is essentially converged. Next, we investigate the SVD–SRG evolution of the deuteron ground–state energy Ed as a function of resolution scale and rank per partial wave in Figure 5.11. The energy Ed should be unaffected by unitary transformations in the two–body system, and therefore, the curves for different resolution scales (λ) should collapse once sufficient components have been included. At r = 12, the deviations from the exact outcome are in the few–keV range, consistent with the slower convergence of the mixing angle ϵ1 with rank r than with other np scattering quantities discussed earlier. At lower ranks r, we observed that the artifacts due to the violation of the unitarity of the evolution get more pronounced as λ decreases. This is because the rank of the interaction starts to grow more rapidly once the SRG begins to reshuffle the long-range pion physics, as explained in Section 5.3.2. In general, the outcomes of SVD–SRG evolution for the EM interaction indicate that the method is effective for chiral N N interactions, especially for those that adopt nonlocal 98 Figure 5.11: Ground-state energy of the deuteron for the SVD-SRG-evolved EM interaction at different resolution scales λ. Here, r refers to the number of components per partial wave. regularization schemes. However, for a hard interaction such as AV18, the story is quite different. When we apply the SVD method based on the singular values of AV18, almost all components must be kept to reproduce observables in the two–nucleon system. We present the AV18 deuteron ground–state energy in Figure 5.12 as an example, where approximately 170 singular components are required to maintain the invariance of the energy. At lower resolution λ, there are plateaus, indicating that many interaction components no longer contribute to the energy due to the decoupling of momentum scales. As λ decreases, these plateaus become more prominent, indicating that low–rank approximations to the evolved AV18 interaction are more accurate, as observed in the Vlow−k approach by Bogner [12] et al., and in line with our discussion in Section 5.3.2. Unfortunately, the structure of unitary transformation (5.12) (cf. Figure 5.6) is too complex to allow us to restrict the evolution only to these components early in the flow. 99 5101520r4.54.03.53.02.52.01.51.0Ed [MeV]EM==10.0 fm1 =4.0 fm1 =3.0 fm1 =2.5 fm1 =2.0 fm1 =1.8 fm1 =1.5 fm1 Figure 5.12: Ground–state energy of the deuteron for the SRG–evolved rank–r approximation of AV18 at different resolution scales λ. 5.3.4 Harmonic Oscillator Basis In the process of preparing the nuclear interaction for many–body methods in config- uration space, at some point a basis change from relative momentum states to (spherical) harmonic oscillator (HO) states is performed. The HO basis is commonly used in nuclear many-body calculations because it allows the exact separation of the center–of–mass and in- trinsic degrees of freedom in many–body states, especially when working in Emax–complete Hilbert spaces. In such spaces, the total energy of the oscillator state is characterized by Emax = ∑i(2ni + li), as explained in previous works (e.g., [39]). In this study, we have two options to perform the SRG evolution: either change the basis of the singular vectors from momentum to HO states through a unitary transformation or directly implement the SVD– SRG in HO representation. The latter option is more practical when evolving three–nucleon forces in an Emax–complete Jacobi–HO representation, because the antisymmetrizer has a block–diagonal structure that makes it easier to implement [34, 35, 50, 75, 85]. 100 0255075100125150175200r543210Ed [MeV]AV18==10.0 fm1 =4.0 fm1 =3.0 fm1 =2.5 fm1 =2.0 fm1 =1.8 fm1 =1.5 fm1 In this study, we have implemented both the momentum–space and harmonic oscillator (HO) approaches and demonstrated that they produce consistent results for the deuteron ground–state energy. This is illustrated for the SVD–SRG evolved EM interaction at a res- olution scale of λ = 2.0 fm−1, as shown in Figure 5.13. Similarly to Figure 5.11, the behavior of Ed should become invariant under SVD–SRG evolution once the rank is sufficiently high. We observe that this occurs at r ≈ 15, which is independent of the HO energy scale ̵hω. For r < 15, the behavior of Ed is similar to that of the momentum–space SVD-SRG, which is also shown for comparison (as seen in Figure 5.11). It should be noted that the magni- tude of the deviations from the exact value weakly depends on ̵hω, and deviations from the momentum–space curve are more significant for the lowest and highest value of ̵hω. These choices correspond to adjusting the effective infrared and ultraviolet cutoffs of the finite HO basis to the scales of the problem (in this case, the deuteron wave function). This topic has received significant attention in recent years in the context of large-basis extrapolations [19, 28, 55, 70, 77, 106]. The key point that we can take away from our analysis is that the SVD–SRG method performs as well in the harmonic oscillator representation as it does in the momentum–space framework, and the conclusions regarding the rank of nuclear interactions hold true, whether for the better (EM) or worse (AV18) case. Furthermore, the representations of the nuclear potential in the momentum and harmonic oscillator bases are quite similar, and although the kinetic energy operator is tri-diagonal rather than diagonal in the harmonic oscillator representation, it can only connect basis states that are close in energy. Therefore, the structure and action of the generator will also be very similar. 101 Figure 5.13: Deuteron ground–state energy of the SVD–SRG evolved EM interaction at λ = 2 fm−1 as a function of the rank. The SVD–SRG and subsequent diagonalization are performed in relative HO bases with different ̵hω. 5.4 Transformation to the Laboratory Frame Moving on from investigating the SVD and SVD–SRG in the two–body system, our next objective is to employ factorized interactions in many-body computations. To that end, we need to analyze the Talmi-Moshinsky transformation of the interaction from the intrinsic frame, which is described by center-of-mass and Jacobi HO states, to the laboratory frame [52]). The singular vectors are coupled to the center–of–mass HO states ∣NcmLcm, uiα;JM ⟩ ≡ ∑ Mcmm ⟨LcmMcmjm ∣ JM ⟩ ∣NcmLcmMcm⟩ ⊗ ∣uiαjm⟩ (5.25) where we have defined the collective partial wave index α ≡ (lsT MT ) for brevity. The corresponding right singular vectors can be constructed analogously. All singular vectors can be transformed individually using the unitary Talmi-Moshinsky transformation matrix. 102 51015202530r2.82.62.42.22.01.8Ed [MeV]momentum=10 MeV=20 MeV=30 MeV=40 MeV=50 MeV Figure 5.14: Matrix elements in the J = 0 and J = 2 channels of the EM interaction at resolution scale λ = 2.0 fm−1, represented in the HO states ∣NcmLcm, n(ls)j; JT MT ⟩. The left and center panels are obtained with an Emax = 16 truncations, while the right panel uses (Ncm, Lcm) = (8, 16) (see text). It is evident from Eq. (5.25) that every singular vector and singular value will appear k times, where k is the number of center–of–mass states. In the extended basis (5.25), the matrix representation of V can be obtained by performing the Kronecker product of the identity matrix in the center–of–mass space with the factorized interaction in the relative space. For illustration, we present the matrices obtained for the EM interaction at λ = 2.0 fm−1 in the J, T, MT = (0, 1, 0) and (2, 1, 0) channels in Figure 5.14. The size and structure of the matrices are dependent on the truncation we apply to the oscillator basis: We can use an Emax truncation and require E = 2Ncm + Lcm + 2n + l = 2n1 + l1 + 2n2 + l2 ≤ Emax (5.26) (with single–particle oscillator quantum numbers ni, li in the laboratory frame), or we can introduce independent truncations Ncm, n ≤ Nmax and Lcm, l ≤ Lmax. In the Emax truncation, we decrease the size of the partial-wave copies as Ncm grows (as depicted in the left and center panels of Figure 5.14), leading to a change in the size of the singular values resulting 103 0255075100125150i020406080100120140160jEM, =2.0 fm1 J=0,EmaxV [fm2]10010110210310410501051041031021011000100200300400500i0100200300400500jEM, =2.0 fm1 J=2,EmaxV [fm2]1001011021031041050105104103102101100050010001500200025003000i050010001500200025003000jEM, =2.0 fm1 J=2,(Nmax,Lmax)V [fm2]1001011021031041050105104103102101100 Figure 5.15: Singular values of matrices shown in Figure 5.14 (plus the J = 0 channel in (Nmax, Lmax) = (8, 16) truncation). The green curves are obtained by compressing duplicate singular values from the full (Nmax, Lmax) = (8, 16) sets. from the projection into a smaller space. In contrast, an (Nmax, Lmax) truncation yields exact replicas of the partial waves (as shown in the right panel). Figure 5.15 illustrates the singular value spectra of the aforementioned matrices. It is not surprising that the factorized matrix with truncation (Nmax, Lmax) has more significant singular vectors than the Emax truncation. However, the former is easier to compress because we would only need to store one representative for each group of identical copies of a given partial wave. On the other hand, the factorized matrix in Emax truncation is not readily compressible due to the projection of the relative partial wave into smaller bases, and the 104 0255075100125150175200104103102101100i [fm2]050010001500200025003000r106105104103102101100i [fm2]EM, =2.0 fm1J=0,(Nmax,Lmax)J=0,(Nmax,Lmax)J=0,EmaxJ=2,(Nmax,Lmax)J=2,(Nmax,Lmax)J=2,Emax resulting change of the singular values. Although the ranks of the Emax and (Nmax, Lmax) matrices are roughly similar overall, a closer examination of the dominant singular vectors indicates a slight advantage for the latter (see top panel of Fig. 5.15). In the case of (Nmax, Lmax), the rank of the interaction in each channel is determined by adding the ranks of the partial waves that meet the selection rules for angular momentum and parity. For instance, in the (J, T, MT ) = (0, 1, 0) channel, the rank is the summation of the ranks of all T = 1 neutron-proton partial waves, where each relative angular momentum j is coupled with the corresponding Lcm to yield a total angular momentum J = 0. For the (J, T, MT ) = (2, 1, 0) channel, all partial waves with ∣Lcm − 2∣ ≤ j ≤ Lcm + 2 are permitted, and the number of allowed couplings will increase with the total J. This observation aligns with a recent investigation that applied tensor factorizations to nuclear interactions and found that their rank increases with J [95]. In summary, we observed that the embedding of partial-wave factorized interactions into a larger space ahead of the Talmi–Moshinsky transformation to the laboratory frame, creates duplicates of the singular values that formally increase the rank of the interaction. Our analysis suggests that the best approach to address this issue in a matrix formulation is to perform the transformation in (Nmax, Lmax) truncation, as this ensures that the duplicates are identical and eliminates the need for additional storage. A more significant version of this problem emerges in the SVD-SRG implementation for three–body forces, where embedding the two–body relative partial waves into the three–body relative partial waves is necessary to track induced forces [35]. We will come back to this problem in chapter 6. 105 Figure 5.16: IMSRG(2) ground–state energies of selected closed–shell nuclei as a function of the flow parameter s for the EM1.8/2.0 N N + 3N interaction (see [33] and text). The SVD– SRG with different ranks (per partial wave) is used to construct the evolved N N component of the interaction. The results are obtained for a HO basis with emax = 8, E3max = 14 and ̵hω = 20 MeV, which is sufficiently close to convergence with respect to basis truncations for the EM1.8/2.0 interaction and these nuclei. 5.5 Many–Body Calculations The final topic we want to address in this chapter is the application of the factorized interactions in many–body calculations. Here, we focus on testing the accuracy of rank–r SVDs by reconstructing the interactions from the factors as we feed them into our many-body methods. In chapter 7, we will attempt to reformulate the IMSRG approach to leverage the factorization directly. In Figure 5.16, we present results of ground–state energy computations for closed–shell 106 104103102101100101102s3029282726252423E [MeV]4Her=10,E=-30.0r=20,E=-29.8r=30,E=-29.8r=40,E=-29.2untruncatedE=-29.2104103102101100101s13012011010090E [MeV]16Or=10,E=-133.2r=20,E=-130.0r=30,E=-129.0r=40,E=-127.0untruncatedE=-127.0104103102101100101s360340320300280260240E [MeV]40Car=10,E=-371.9r=20,E=-352.1r=30,E=-345.8r=40,E=-341.2untruncatedE=-341.2104103102101100101s1000900800700600E [MeV]132Snr=10,E=-1007.2r=20,E=-862.0r=30,E=-831.5r=40,E=-821.4untruncatedE=-821.4 nuclei in the IMSRG(2) approach [37, 39]. They were obtained using the EM1.8/2.0 inter- action, which encompasses the EM interaction progressed to λ = 1.8 fm−1 and an NNLO 3N interaction with cutoff λ = 2.0 fm−1. The low–energy constants of this interaction have been calibrated to the triton binding energy and 4He charge radius [33, 76]. Although this interaction is not entirely consistent from the chiral EFT perspective, it has been empiri- cally successful in describing the ground–state energies of a large number of nuclei, albeit while underestimating radii by a few percent (see [40] and the cited literature, particularly [92]). Thus, we view calculations with EM1.8/2.0 as a benchmark for the performance of the factorized interactions under “realistic” conditions. If we solely rely on the SVD–SRG evolved N N interaction, it is well known that nuclei will be overbound and much too small (see [39] and the behavior of the N N -only triton energy in Fig. 2.4, for example). We find that we need 30–40 components per partial wave in the SVD–SRG interaction to recover the results without factorization accurately, with no attempts to fine tune the decomposition at this stage. This observation is in line with our findings in the two–nucleon system. This encompasses both the Hartree–Fock calculations used to prepare the reference states for normal ordering (cf. Section 1.3.7) and the full IMSRG(2) flow of the ground–state energy concerning as a function of the flow parameter s (cf. Chapter 3 and Refs. [37, 39]). As in our applications in the two-nucleon system, the need to include 30–40 components is primarily driven by the SVD of the Coulomb interaction VC between the protons — for calculations with the strong interaction alone, 5 to 10 components are sufficient for highly precise reproduction of the full flow, but the results do not describe realistic nuclei. Looking in more detail, it is evident that the convergence of the studied nuclei is not completely uniform. For example, the ground–state energy of 4He remains almost the same as the rank increases form r = 20 to 30, while 16O and 40Ca experience a significant gain 107 of several hundred keV and a few MeV, respectively. Additionally, we note that further increasing the rank to 40 results in more considerable changes in the ground–state energies of 4He and 16O, while the change in 40Ca is smaller. As noted in our discussion of the pp phase shifts in Section 5.3.3, these findings can be attributed to the inclusion of interaction components based on the size of their associated singular values alone, leading to an alterna- tion between components of the nuclear and Coulomb interactions that affects these nuclei differently, given the significantly different proton numbers. Furthermore, 4He, 16O and 40Ca probe the interaction at increasingly higher (laboratory–frame) angular momenta, hence it is reasonable to expect that specific interaction components may only become relevant for heavier masses. We have confirmed that the ground–state energies of the chosen nuclei reach convergence by gradually increasing the rank from 40 to the full rank. We have not included the graphs of these curves in Figure 5.16 to avoid confusion, but we refer the reader to Figure 5.17, which illustrates that the rank of 40 produces a precise and convergent replication not just of the ground-state energies but also of the root-mean-square radii of the studied nuclei. To explore the performance of the SVD–SRG for general observables, we construct the evolved mean-square radius operator by starting from the initial form R2 = 1 A2 A ∑ i=1 (ri − Rcm)2 , Rcm = 1 A A ∑ i=1 ri . (5.27) and using Eq. (5.12) to derive the unitary transformation in the two–body system in terms of the singular vectors of the truncated SVD. The expectation value of the evolved R2 is then computed in IMSRG(2). Figure 5.17 illustrates the relationship between the rank of the SVD and root–mean– 108 Figure 5.17: Accuracy of the SVD–SRG evolution for the root–mean–square radius of se- lected closed–shell nuclei: The curves show the ratio R/Rall as a function of the rank r, where Rall is obtained from a conventional matrix–based SRG evolution. The results for the bare operator are shown for comparison (dash–dotted lines). All calculations were performed with the EM1.8/2.0 interaction, using the IMSRG(2) truncation, eemax = 8, E3max = 14 and ̵hω = 20 MeV (cf. Figure 5.16). square radii R ≡ √ ⟨R2⟩ of 4He, 16O and 40Ca in the IMSRG(2). As in the case of the ground- state energy energy, a rank of 30 to 40 is adequate to achieve an accurate reconstruction of the unitary transformation. We also include the radii obtained with the unevolved R2 operator for comparison, which vary by about 5%, 0.5% and 0.1% in 4He,16 O and 40Ca, respectively. Since the SRG evolution aims for physics at high–momentum or short–range, its impact on a long–ranged operator like R2 is weak and negligible when compared to other sources of uncertainty at present. 109 20406080100120140r1.001.051.101.151.20R(r)/Rall4He16O40Ca==2.0 fm1 We conclude our discussion by stating although R2 is undoubtedly one of the simplest observables to examine besides the energy, we anticipate no problems in using the factorized unitary transformation for more complicated operators, like in applications for electroweak transitions [26, 79, 111]. The reason is that ultimately, U (s) is entirely determined by the characteristics of the Hamiltonian. 5.6 Conclusions and Directions for Future Research In this chapter, we used the SVD to perform principal-component analyses of nucleon- nucleon interactions, using the chiral N3LO interaction by Entem and Machleidt [22] and the Argonne V18 interaction [110] as representative examples. Our findings demonstrated that the EM interaction can be easily represented by a low–rank model that is constructed by truncating the SVD based on the magnitude of the singular values. Conversely, the AV18 interaction proved to be a much more challenging case because of its local nature and the resulting structure in momentum space, as well as its high implicit resolution scale. We have merged the SVD and SRG techniques and demonstrated that the low–rank representation can be effectively evolved to lower resolution scales. However, the efficiency of this approach compared to traditional SRG depends on the rank of the original interaction. While the evolution of two-nucleon interactions using SRG is a well–established process, the next step is to extend these techniques to three-nucleon forces, where much more substantial efficiency improvements can be achieved — this will be the subject of Chapter 6. We have carried the factorized representation of N N interactions through the various stages of the many-body theory workflows, and identified challenges, in particular when switching from the intrinsic/center-of-mass frame to the laboratory frame. At the final 110 stage, however, we resorted to reconstructing the input matrix elements from the factors, since formulating many–body methods capable of exploiting thes factorization more directly is a highly challenging task — we will come back to this issue in Chapter 7. While our primacy focus lies in the aforementioned developments, there are additional interesting directions for future research: 1. In light of the adverse impact of the Coulomb interaction on the achievable rank compression in pp channels, it is worth exploring alternative treatments of VC for addressing this issue. 2. The consistent ranks we observed across all partial waves of the interaction suggests that it may be a manifestation of a few essential operators projected into the different channels. Investigating such a possible connection may offer new ways of constructing compact operator bases for SRG and IMSRG flows. 3. The implementation of the Talmi-Moshinsky transformation on the SVD factors demon- strated how embedding an operator in larger product Hilbert spaces artifically increases the rank of its matrix representations through redundant copies. Given that many– body Hilbert spaces themselves possess a product structure, tensor representations appear to be a particularly suitable candidate for efficiently handling the physical in- formation encoded in an operator. 111 Chapter 6 SVD-SRG for Three-Nucleon Interactions Two-nucleon interactions alone are insufficient to accurately describe low-energy states in light nuclei, hence we are forced to include at least three-body (3N ) interactions in our calculations (and in the worst case, even higher many-body forces, although we expect them to only matter in systems with large A or infinite matter, if at all, as discussed in Sections 2.4). The computational cost and memory requirements for handling 3N forces are substantial, which is why we will explore their factorization and the development of an SVD-SRG for 3N interactions in this chapter. 6.1 Three-Body Basis States We will begin by focusing on the fundamental basis set, which pertains to the eigenstates of the harmonic oscillator (HO). These eigenstates are established in terms of Jacobi coor- dinates, which disentangle the relative and center-of-mass motion in the three-body system. Several options exists for constructing these coordinates — here, we adopt the same defini- tion and conventions as in previous references [16, 35, 72, 73]. The spatial Jacobi coordinates 112 are defined recursively as √ √ ξ0 = ξn = 1 A n n + 1 [r1 + r2 + ⋯ + rA] [ 1 n (r1 + r2 + ⋯ + rn) − rn+1] (6.1) and the corresponding Jacobi momenta are [p1 + p2⋯ + pA] √ √ π0 = πn = 1 A n n + 1 [ 1 n (p1 + p2 + ⋯ + pn) − pn+1] (6.2) where ri and pi represent the position and momentum of the i-th particle, respectively, within the system. The coordinates ξ0 and π0 are directly related to the center-of-mass properties of the entire A-body configuration. For 0 < n ≤ A − 1, the Jacobi coordinates ξn and πn capture the relative motion within the A-body system. From the definition, we see that the nth Jacobi coordinates describe the distance of the (n + 1)th particle from the center-of-mass of the first n particles. In Section 5.4, we previously encountered the HO basis states for the intrinsic frame in a two-nucleon system, which are given by ∣NcmLcmMLcm⟩ ⊗ ∣N1; (L1S1)J1MJ1; T1MT1⟩a . (6.3) Here, the subscript “cm” identifies the quantum numbers of the center-of-mass state encoding the collective motion of the system as a whole, rather than the individual particles within it. The quantum numbers N1 and L1 refer to the radial and orbital angular momentum of 113 an oscillator state in the first Jacobi coordinate. S1, T1 and MT1 denote the spin and isospin quantum numbers, as usual, and the angular momentum of the relative two-nucleon state and its projection are denoted as J1 and MJ1, respectively. We use the subscript a in Eq. (6.3) to indicate that the state describing the two-nucleon (sub)system has been antisymmetrized. This implies that its quantum numbers satisfy the condition (−1)L1+S1+T1 = 1. Generalizing (6.3) to the three-body system, we introduce the basis states (see Ref. [16, 35] for a detailed discussion) ∣NcmLcm⟩ ⊗ ∣N1N2; [(L1S1)J1, (L2 1 2 )J2]J12; (T1 1 2 )T12⟩ a1 . (6.4) Let us unpack the structure of these states in detail: • N2 and L2 are the radial and orbital angular momentum quantum numbers associated with the second Jacobi coordinate. • Based on the radial and orbital angular momentum quantum numbers, we can see that the total HO energy of (6.4) is given by E = Ecm + E12, where Ecm = 2Ncm + Lcm represents the center-of-mass contribution and E12 corresponds to the relative HO energy given by E12 = (2N1 + L1) + (2N2 + L2). • The total relative angular momenta associated with the two coordinates, J1 and J2, are coupled to J12. Since the interactions are rotationally invariant, we drop the angular momentum projection quantum numbers for brevity. • The third nucleon adds additional spin and isospin degrees of freedom to (6.4), but it is clear that we must always have S2 = 1 2 and T2 = 1 2, and we indicate so explicitly. The total isospin of the three-nucleon state is obtained by coupling T1 = 0, 1 and T2, which 114 means that we can have T12 = 1 2, 3 2. • The subscript a1 specifically indicates that this state is only partially antisymmetric, namely with respect to the exchange of the first two nucleons that are associated with the first Jacobi coordinate. The partially antisymmetric three-body basis is introduced because we need to embed the evolving two-body interaction into the three-body system in order to track the induced interactions discussed in Section 2.4. Of course, we must eventually antisymmetrize the three-body states completely if we want to obtain physical results in few- and many-body calculations. The most practical approach is to work with the antisymmetrizer A: We can represent A in the partially antisymmetric Jacobi-HO basis (as given in (6.4)) and diagonalize the corresponding matrix to obtain the eigenstates of this operator [16, 35, 75]. The antisymmetrizer is block-diagonal in the Jacobi-HO basis, as briefly mentioned in Chapter 5. Permutations of the particles must not change the intrinsic energy of the state, hence E12 is conserved. The antisymmetrization does not affect the rotational symmetry of the states 6.4 either, hence the 3N interaction will be diagonal in the angular momentum J12, and it is also diagonal in parity (−1)π = (−1)L1+L2. Last but not least, one usually also ignores isospin-breaking effects in the 3N interaction in many-body calculations because they are small compared to most other sources of uncertainty. Thus, the 3N interaction will also be diagonal in isospin T12. As a result, the diagonalization of A can be performed separately for each (E12, J12, π12, T12)-block. The fully antisymmetric eigenstates are associated with the eigenvalue (−1), and these eigenstates can be presented as linear combinations of the partially antisymmetric states (6.4) in each block. he introduction of the antisymmetric Jacobi-HO basis also has a practical benefit, because 115 the number of interaction matrix elements is reduced compared to the merely partially antisymmetric basis. This reduction is about an order of magnitude. As a final remark, we note that the antisymmetrization procedure discussed here cannot be directly applied in arbitrary basis sets, because the antisymmetrizer will not be block-diagonal in general and other approaches for ensuring the proper antisymmetry are needed [34, 35, 75, 105]. A particularly important example is the 3N partial wave momentum basis [34, 35]. 6.2 SRG Evolution in the Three-Nucleon System Let us now discuss the SRG evolution for a three-nucleon system. As discussed previously, we can use what we call the “α” basis1, which is partially anti- symmetric under exchanging particles 1 and 2, but has no definitive symmetry with respect to exchanges involving particle 3. Whenever we encounter an interaction of the third par- ticle with the other particles in the system, e.g., V23, we can rewrite it in terms of V12, the interaction between particles 1 and 2: V23 = P13V12P13 = P12P23P12V12P12P23P12 V13 = P23V12P23 (6.5a) (6.5b) Here the Pij are basic permutation operators that swap particles i and j. In this expression, we have used P13 = P12P23P12: To swap particles 1 and 3 in ∣abc⟩ is equivalent to swap- ping particles 1 and 2 to obtain ∣bac⟩, then swapping particles 2 and 3 (∣bca⟩), and finally 1 and 2 again, which ultimately yields ∣cba⟩), as desired. Note that the permutation oper- 1The collective index α is introduced to represent the “partial wave” quantum numbers in the states 6.4, i.e., α12 = {L1, S1, J1, L2, J2, J12, T1, T12}. 116 ators commute with their corresponding 2-body potentials, e.g., [P12, V12] = 0, etc. Due to antisymmetry, the matrix elements of V23 and V13 are the same: ⟨α′∣ V23 ∣α⟩ = ⟨α′∣ P12P23V12P23P12 ∣α⟩ = (−1)2 ⟨α′∣ P23V12P23 ∣α⟩ (6.6) The same works for the kinetic energy. The P23 operator is unitary and Hermitian, and thus P 2 23 = 1. This implies that its only eigenvalues are +1 and −1. Denoting the eigenvectors to eigenvalue -1 as ∣E, i, J12, π12, T12⟩. we can write P23 in spectral representation as P23 = 1J12,T12 − 2 ∑ ∣E, i, J12, π12, T12⟩ ⟨E, i, J12, π12, T12∣ (6.7) E,i Let us now construct the SRG flow equation for the 3-body part of the evolving Hamil- tonian (induced + initial). It is given by [34, 35, 50] dH ds = [η, H] = [[T, V ], H] =[[T12 + T13 + T23, V12 + V13 + V23 + V123], T12 + T13 + T23 + V12 + V13 + V23 + V123] (6.8) As before, T is explicitly kept constant, and induced terms are absorbed into the evolving interaction. We can split the flow equation into four classes of terms: dV ds ≡ dVa ds + dVb ds + dVc ds + dVd ds (6.9) 117 with and dVa ds dVb ds dVc ds dVd ds ≡[[T12, V12], T12 + V12] + [[T13, V13], T13 + V13] + [[T23, V23], T23 + V23] , (6.10a) ≡[[T12, V12], T13 + V13] + [[T12, V12], T23 + V23] + [[T13, V13], T12 + V12] + [[T13, V13], T23 + V23] + [[T23, V23], T12 + V12] + [[T23, V23], T13 + V13] , (6.10b) ≡[[T12, V13 + V23], T12 + T13 + T23 + V12 + V13 + V23] + [[T13, V12 + V23], T12 + T13 + T23 + V12 + V13 + V23] + [[T23, V12 + V13], T12 + T13 + T23 + V12 + V13 + V23] , (6.10c) ≡[[T12 + T13 + T23, V123], T12 + T13 + T23 + V12 + V13 + V23] + [[T12 + T13 + T23, V123], V123] + [[T12 + T13 + T23, V12 + V13 + V23], V123] (6.10d) The first class, dVa ds , described the evolution of the pairwise two-nucleon interactions, embedded in the three-nucleon system, in the presence of a spectator nucleon. The con- tributions dVb ds and dVc ds are given in terms of two-body operators, but all three nucleons are actively affected by the evolution, hence they generate induced three-body forces. Finally, dVd ds is the induced force due to the flowing three-body interaction. We can simplify the expression: In the end, we are only interested in the antisymmetric 118 matrix elements ⟨E′, i′, J12, π12, T12∣ dV ds ∣E, i, J12, π12, T12⟩ = ⟨E′, i′, J12, π12, T12∣ A dV ds A ∣E, i, J12, π12, T12⟩ . (6.11) For any permutation operator Pij, we have APij = PijA = −A, so whenever we can pull out one of the permutations, we can replace it by a minus sign. Applying this idea to the first class of terms, Eq. (6.10a), we get A dVa ds A = A[[T12, V12], T12 + V12]A + AP23[[T12, V12], T12 + V12]P23A + = AP13[[T23, V23], T23 + V23]P13A = 3A[[T12, V12], T12 + V12] . (6.12) When integrated, this will yield the evolved two-nucleon interaction embedded in the three- nucleon space, which will be subtracted from the evolved Hamiltonian along with the kinetic terms to obtain the evolved three-nucleon interaction (see below). For dVb ds (Eq. (6.10b)), we obtain A dVb ds A = A[[T12, V12], T13 + V13]A + A[[T12, V12], T23 + V23]A + = A[[T13, V13], T12 + V12]A + A[[T13, V13], T23 + V23]A + = A[[T23, V23], T12 + V12]A + A[[T23, V23], T13 + V13]A = AP12[[T12, V12], T23 + V23]P12A + A[[T12, V12], T23 + V23]A + = AP12[[T23, V23], T12 + V12]P12A + AP23[[T12, V12], T23 + V23]P23A + 119 = A[[T23, V23], T12 + V12]A + AP23[[T23, V23], T12 + V12]P23A + = 3A[[T12, V12], T23 + V23]A + 3A[[T23, V23], T12 + V12]A = 3A[[T12, V12], T23 + V23]A + 3AP12P13P23[[T12, V12], T23 + V23]P23P13P12A = 6A[[T12, V12], T23 + V23]A , (6.13) and dVc ds (Eq. (6.10c)) can be simplified to A dVc ds A = A[[T12, V13 + V23], T12 + T13 + T23 + V12 + V13 + V23]A + = A[[T13, V12 + V23], T12 + T13 + T23 + V12 + V13 + V23]A + = A[[T23, V12 + V13], T12 + T13 + T23 + V12 + V13 + V23]A = [A[T12, V13 + V23]A, A(T12 + T13 + T23 + V12 + V13 + V23)A] + = [A[T13, V12 + V23]A, A(T12 + T13 + T23 + V12 + V13 + V23)A] + = [A[T23, V12 + V13]A, A(T12 + T13 + T23 + V12 + V13 + V23)A] = [A[T12, V13 + V23]A, 3A(T12 + V12)A] + = [A[P23T13P23, P23(V12 + V23)P23]A, 3A(T12 + V12)A] + = [A[P13T23P13, P13(V12 + V13)P13]A, 3A(T12 + V12)A] = 3[A[T12, V13 + V23]A, 3A(T12 + V12)A] = 9[A[T12, P12V13P12 + V23]A, A(T12 + V12)A] = 18[A[T12, V23]A, A(T12 + V12)A] (6.14) 120 Here, we have used that for any unitary operator U we can write [B, C] = BC − CB = U U †BU U †CU U † − U U †CU U †BU U † = U [U †BU, U †CU ]U † (6.15) and for any idempotent operator A2 = A that commutes with C, we have A[B, C]A = ABCA − ACBA = ABA2CA − ACA2BA = [ABA, ACA] (6.16) Finally, the terms involving the flowing three-body interaction (Eq. (6.10d)) can be rewritten as A dVd ds A = A[[T12 + T13 + T23, V123], T12 + T13 + T23 + V12 + V13 + V23]A + = A[[T12 + T13 + T23, V123], V123]A + = A[[T12 + T13 + T23, V12 + V13 + V23], V123]A = 3[[AT12A, V123], 3A(T12 + V12)A + V123] + = 9[[AT12A, AV12A], V123] (6.17) Since dVa ds is the flowing two-nucleon part, the rest must be the ODE for the genuine three- body terms. Altogether, we have dH ds = d ds (V [2] + V [3]) =3A[[T12, V12], T12 + V12]A 121 + 6A[[T12, V12], T23 + V23]A + 18[A[T12, V23]A, A(T12 + V12)A] + 3[[AT12A, V123], 3A(T12 + V12)A + V123] + 9[[AT12A, AV12A], V123] (6.18) In practice, the simplest approach to implement the evolution is to construct the initial Hamiltonian in the antisymmetric three-body basis, and evolve it as a whole to the desired flow parameter s (or resolution scale λ)). Then we embed the two-body interaction V [2] at the same s or λ, which can be evolved separately, into the antisymmetric three-body basis, and subtract both T and V [2] from the Hamiltonian to obtain the evolved three-body interaction V [3]. Because of the symmetries discussed above, we can perform the evolution in (J12, π12, T12) blocks, which helps us alleviate the computational effort. Still, the dimensions of these blocks are significant if we need to allow large energy quantum numbers E12 — state-of the art implementations typically use E12 ≲ 40 in the lowest partial waves, and ramp down slightly as the angular momentum increases [35, 85]. In contrast, evolutions for two- body interactions in Jacobi-HO representation can be performed with energy cutoffs as high as E1 ∼ 200 − 300. This seems to be sufficient — if challenging — for currently used chiral interactions with soft, nonlocal regulators, but there is a need to accelerate the three-body SRG evolution to allow for other, harder types of interactions as we explore the chiral power counting, and to enable the exploration of the parameter space for the chiral LECs. 6.3 Singular Value Decomposition of 3N Interactions Let us now proceed and apply (randomized) SVDs to typical three-nucleon interactions. In Fig. 6.1, we show the singular spectra of a 3N interaction in various (J, π, T ) partial 122 (a) T = 1 2 J = 1 2 π = 0 (b) T = 1 2 J = 3 2 π = 1 (c) T = 1 2 J = 5 2 π = 0 (d) T = 3 2 J = 5 2 π = 0 Figure 6.1: Singular-value spectra of a chiral 3N interaction in the J, π, T partial waves indicated by the sub-captions. This is an N2LO three-nucleon interaction with cutoff Λ = 2.0 fm−1 that has been fit to accompany the SRG-evolved EM N3LO N N interaction (SRG λ = 2.0 fm−1) [33]. The red circular dot indicates that the running sum of singular values constitutes 99% of the total, whereas the red square dot indicates when we have reached 99.9% of the total. 123 01000200030004000N108106104102100iNi = 54, i = 0.007851066Nj = 153, j = 0.000404635Dimension = 4263Singular Values01000200030004000500060007000N108106104102100iNi = 120, i = 0.002901565Nj = 328, j = 0.000174424Dimension = 7511Singular Values0200040006000800010000N108106104102100iNi = 204, i = 0.001079518Nj = 523, j = 0.000070266Dimension = 11340Singular Values010002000300040005000N108106104102100iNi = 85, i = 0.000327486Nj = 241, j = 0.000018617Dimension = 5670Singular Values Figure 6.2: Singular values of the initial (solid black) and evolved chiral 3N interaction interactions (dashed red) are displayed. We mark the points when the running sum of singular values reaches 99% and 99.9% of the total, as in Fig. 6.1. waves, as indicated in the sub-captions (we drop the label “12” referring to the Jacobi coordinates for brevity). Specifically, we consider an N2LO interaction with a cutoff value of Λ = 2.0 fm−1 that has been tailored to complement the SRG-evolved EM N3LO nucleon- nucleon interaction with SRG resolution scale λ = 2.0 fm−1 (corresponding to s = λ−4 = 0.0625 fm4) — see Ref. [33] for details. It is evident from the figures that the singular values exhibit a rapid decay: With a few tens or hundreds of singular values, we exhaust 99% or 99.9% of the total sum of singular values in the shown partial waves, which are representative for the overall behavior. We expect that these modes should be sufficient to capture the essential physical information of the interaction. Unfortunately, the SRG evolution seems to spoil the low-rank structure of the evolved interactions, as illustrated for a representative partial wave in Fig. 6.2. The culprit is almost certainly the interplay of the two- and three-nucleon terms in the three-body flow equations: 124 01000200030004000500060007000108106104102100102iT = 1/2, J = 3/2, P = 1V[3](0)V[3](s) As we saw in chapter 5, the evolving N N interaction maintains a low rank, but when we embed it in the three-body system, we form a Kronecker product of this matrix with the identity matrix for the third particle / second Jacobi coordinate. This creates duplicates of the singular-value spectrum, just as in the case of the Talmi-Moshinsky transformation discussed in Sec. 5.4. The situation is then complicated further by the application of the antisymmetrizer, which leads to a “spreading” of the matrix elements to sectors that were previously zero, and the resulting matrix then appears in products with T and V [3]. Because the issues we encountered in the present application as well in the Talmi- Moshinsky transformation stem from the Kronecker product structure of our Hilbert spaces, we hope that tensor factorization techniques that explicitly take such structures into account can offer a viable path forward. They will be explored in future work. 6.4 Performance of the Randomized SVD We conclude this chapter with a discussion of the performance of the various RSVD algorithms in our applications for three-nucleon interactions. In Fig. 6.3, we compare the run times for the regular SVD with the RSVD variants described in Sec. 4.4. To establish a meaningful comparison with RSVD, we adopted a truncation criterion wherein 99.9% of the total singular values were retained, while maintaining a fixed oversampling parameter of 10. Figure 6.3(a) reveals that the RSVD algorithms is typically about two orders of magnitude faster than the regular SVD for sufficiently large matrices with dimensions greater than 2000- 3000, which are common in our applications for 3N forces (as well as the IMSRG, see Chapter 10). For smaller matrices, the speedup is still at least one order of magnitude. If needed, a few power iterations can be performed to refine the singular value spectrum at a small 125 (a) (b) Figure 6.3: Run time as a function of dimension of the 3N matrices for SVD and multiple RSVD variants (cf. Sec. 4.4.) Panel a: SVD vs. RSVD algorithm 1 with different numbers of power iterations q. Panel b: Comparison of the main RSVD algorithms without power iteration. computational overhead (a prefactor of 1-1.5 for up to q = 2 power iterations). In 6.3(b), we compare the performance of our main algorithms (without power iterations), and find that algorithm 2 is approximately a factor 2 more effective than algorithms 1 and 3 — we will see below that this seems to come at the cost of a slight loss of precision in the tail of the singular-value spectrum that is likely irrelevant in practical applications. In Figure 6.4, we compare the singular-value spectra of SVD and RSVD for the three- nucleon interaction in the (J π, T ) = ( 3 2 + , 1 2). To reduce computational complexity, we trun- cate the singular values to retain 99.9% of the total running sum, which translates to 328 components out of the original 7511 (cf. Figs. 6.1 and 6.2). The oversampling parameter is k = 10, as before. We see that all RSVD algorithms give an excellent reproduction of the full singular value spectrum with 1-2 orders of magnitude less computational cost. There is a slight loss of precision in the tail of the spectrum, which can be corrected by performing 126 200040006000800010000Dim101102103104T(s)SVDRSVD-1,q=0RSVD-1,q=1RSVD-1,q=2200040006000800010000Dim101102T(s)RSVD-1,q=0RSVD-2,q=0RSVD-3,q=0 (a) (b) (c) Figure 6.4: Comparison of singular value spectra for the (J π, T ) = ( 3 , 1 2) partial wave of 2 the EM1.8/2.0 3N interaction from SVD and RSVD with (q = 1) and without (q = 0) power iterations. We only show the first 328 components — the full dimension of the matrix is 7511 (see text). + 127 050100150200250300104103102101100SVDRSVD-1, q=0RSVD-1, q=1050100150200250300104103102101100SVDRSVD-2, q=0RSVD-2, q=1050100150200250300104103102101100SVDRSVD-3, q=0RSVD-3, q=1 a single power iteration to refine the spectrum. As mentioned before, there is some noisy behavior for RSVD algorithm 2. 128 Part III Factorization in Many-Body Systems 129 Chapter 7 SVD and IMSRG We now move from the few-body systems discussed in the previous chapters to many- body systems, with the goal of applying SVDs to factorize the quantities that appear in the IMSRG flow equations. In Section 5.4, we saw that it is not straightforward to “transfer” low- rank structures that were identified in the intrinsic frame to the laboratory frame, because embeddings of the factors into larger Hilbert spaces and basis transformations can inflate the rank of the interaction again. However, there are compelling reasons for exploring SVDs of the normal-ordered Hamiltonian and generator appearing in the IMSRG: As discussed in Section 3.4, the latter has the structure η = ∑ ai ηai{a† aai} + 1 4 ∑ abij ηabij{a† aa† bajai} − H.c. , (7.1) and the N 2 p × N 2 h matrix ηabij is flat and wide in realistic applications because Nh ≪ Np. Thus, its numerical rank should satisfy k ≤ N 2 h, which is much smaller than the naive rank N 4 for normal-ordered two-body operators. In the present chapter, we reformulate the IMSRG(2) (cf. Section 3) in a way that allows us exploit a rank reduction, if revealed through SVD. In subsequent chapters, we will apply this new formulation to a schematic model, as well as infinite nuclear matter. We note that similar ideas have been explored in Coupled Cluster calculations in quantum chemistry; after 130 first attempts in Refs. [44, 54], significant progress has been made in just the last couple of years, in parallel to the work described here [45, 46, 62, 63, 78]. 7.1 The IMSRG(2) in Rank-Reduced Form To start, we assume that fpq = ∑ µν p f µνuν uµ q Γpqrs = ∑ αβ U α pqΓαβU β rs (7.2) (7.3) where f and Γ are the one- and two-body parts of the normal-ordered Hamiltonian, as dis- cussed in Section 1.3.8. The coefficients uµ p and U α pq are the elements of basis transformations in the one- and two-body spaces — the idea is that these coefficients could be obtained as singular vectors of a truncated SVD of f , Γ or other appropriate one- and two-body operators. Denoting the ranges of the transformed one- and two-body bases by d and D, respectively, we ideally would like that d ≪ N , D ≪ N 2 , (7.4) so that there will be a significant compression in comparison to the working basis. While we have touched upon the tensor structure of our operators and many-body wave functions in Chapters 1 and 5, we work with matrix representations in the following — the exploration of tensor representations will be left for future work. 131 In analogy to Eqs. (7.2) and (7.3), we define the transformed generator as ηpq = ∑ µν p ηµνuν uµ q , ηpqrs = ∑ αβ U α pqηαβU β rs . (7.5) (7.6) In our first attempt at merging SVD and IMSRG, we take a slightly different approach than in the SVD-SRG: Instead of evolving both the singular value matrices and left and right singular vectors, We assume that the transformations uµ p and U α pq defined by the singular vectors remain fixed, and we evolve the (reduced) quantities f µν and Γαβ. Let us now consider the zero-body flow equation (3.9) and plug in Eqs. (7.2)–(7.6). First, we split the equation into two contributions: dE ds (np − nq)ηpqfqp = ∑ pq ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ 1○ + Defining the transformed density matrix ∑ pqrs ηpqrsΓrspqnrns¯np¯nq 1 2 ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ 2○ . (7.7) ρµν ≡ ∑ p p npuν uµ p , (7.8) and assuming that the reduced basis uµ is still orthonormalized (which would be true for a basis of singular vectors), uµ p uν p = δµν , ∑ p (7.9) 132 we obtain 1○ = ∑ ab (na − nb) ∑ µν ∑ µ′ν′ aηµνuν uµ b uµ′ b f µ′ν′ uν′ a = ∑ µν ∑ µ′ν′ (∑ a anauν′ uµ a ∑ b b uµ′ uν b − ∑ b b nbuµ′ uν b ∑ a auν′ uµ a ) ηµνf µ′ν′ = ∑ µν ∑ µ′ν′ (ρµν′ δνµ′ − ρνµ′ δµν′) ηµνf µ′ν′ = ∑ µµ′ν′ (ηµµ′ f µ′ν′ ρν′µ) − ∑ µµ′ν (ρνµ′ f µ′µηµν) (ηλµf µν − f λµηµν) ρνλ = ∑ λµν = Tr (ηf ρ − f ηρ) . Similarly, we can introduce the two-body particle and hole density matrices Pαβ ≡ ∑ pq αβ P ≡ ∑ pq and rewrite the second term as U α pqnpnqU β pq = ∑ ij U α pq ¯np¯nqU β pq = ∑ ab ijninjU β U α ij , ab¯na¯nbU β U α ab , 2○ = = = = 1 2 1 2 1 2 1 2 ∑ pqrs ∑ αβ ∑ α′β′ U α pqηαβU β rsU α′ rs Γα′β′ U β′ pq npnq ¯nr ¯ns ∑ αβ ∑ α′β′ (∑ pq pqnpnqU β′ U α pq ) (∑ rs rs¯nr ¯nsU α′ U β rs ) ηαβΓα′β′ βα′ ηαβP Γα′β′ Pβ′α ∑ αβ ∑ α′β′ Tr (ηPΓP) . 133 (7.10) (7.11a) (7.11b) (7.12) Putting the two terms together, we obtain dE ds = Tr (ηf ρ − f ηρ) + 1 2 Tr (ηPΓP) , (7.13) and we see that the flow of the energy can be easily computed in terms of the reduced quantities. Ultimately, this is not surprising, because the RHS of the flow equation is given in terms of traces that must be invariant under basis changes. Next, we consider the one-body flow equation (3.10). First, we note that na − nb = na − nb + nanb − nanb = na(1 − nb) − nb(1 − na) = na¯nb − ¯nanb , (7.14) hence we can rewrite the flow equation as dfpq ds = ∑ r (1 + Ppq)ηprfrq + ∑ rs (nr ¯ns − ¯nrns)(ηrsΓsprq − frsηsprq) + 1 2 (nrns¯nt + ¯nr ¯nsnt)(1 + Ppq)ηtprsΓrstq . ∑ rst (7.15) Applying the transformation to the reduced representation and summing, we obtain df µν ds uµ p (∑ r (1 + Ppq)ηprfrq) uν q = ∑ pq ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ 1○ uµ p (∑ rs (nr ¯ns − ¯nrns)(ηrsΓsprq − frsηsprq)) uν q + ∑ pq ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ 2○ 1 2 uµ p ( (nrns¯nt + ¯nr ¯nsnt)(1 + Ppq)ηtprsΓrstq) uν q + ∑ pq ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ 3○ ∑ rst = (ηf − f η)µν + (ΥΥηΓ − ΥΥηΓ − ΥΥf η + ΥΥf η) µν ′ (PΞηΓ − PΞηΓ + PΞ ηΓ − PΞ′ηΓ) µν + 1 2 (7.16) 134 where we have defined Υµνα ≡ ∑ pq p uν uµ q U α qpnq , µνα Υ ≡ ∑ pq p uν uµ q Uqp¯nq , (7.17) and In detail, Ξαβµν ≡ ∑ pqr U α rpU β rquµ p uν q nr , αβµν Ξ ≡ ∑ pqr U α rpU β rquµ p uν q ¯nr . 1○ = ∑ ij uµ i ∑ a (ηiafaj + ηjafai) uν j = ∑ ij ∑ a ∑ µ′ν′ ∑ µ′′ν′′ i (uµ′ uµ i ηµ′ν′ a uµ′′ uν′ a f µ′′ν′′ j + uµ′ uν′′ j ηµ′ν′ a uµ′′ uν′ a f µ′′ν′′ uν′′ i ) uν j = ∑ ij ∑ µ′ν′ ∑ µ′′ν′′ i (uµ′ uµ i ηµ′ν′ δν′µ′′ f µ′′ν′′ j + uµ′ uν′′ j ηµ′ν′ δν′µ′′ f µ′′ν′′ uν′′ i ) uν j = ∑ ij ∑ µ′ν′ν′′ i (uµ′ uµ i ηµ′ν′ f ν′ν′′ j + uµ′ uν′′ j ηµ′ν′ f ν′ν′′ uν′′ i ) uν j = ∑ ij ∑ µ′ν′′ i (uµ′ uµ i (ηf )µ′ν′′ j + uµ′ uν′′ j (ηf )µ′ν′′ uν′′ i ) uν j δµµ′(ηf )µ′ν′′ δν′′ν + δνµ′(ηf )µ′ν′′ δν′′µ = ∑ µ′ν′′ = (ηf )µν + (ηf )νµ = (ηf − f η)µν where in the last equality, we used the symmetry of f and antisymmetry of η, 2○ = ∑ ij uµ i (∑ ab (na − nb)(ηabΓbiaj − fabηbiaj)) uν j 135 (7.18) (7.19) = ∑ ij (na − nb) (uµ i uν j uµ′ a ηµ′ν′ ∑ ab uν′ b U α biΓαβU β aj − uµ i uν j uµ′′ a f µ′′ν′′ b U α′ uν′′ bi ηα′β′ U β′ aj ) (Υνµ′βΥ µν′α νµ′β − Υ Υµν′α) ηµ′ν′ Γαβ− (Υνµ′′β′ µν′′α′ Υ νµ′′β′ − Υ Υµν′′α′) f µ′′ν′′ ηα′β′ ∑ µ′′ν′′ ∑ α′β′ (Υνµ′′β′ µν′′α′ Υ νµ′′β′ − Υ Υµν′′α′) (ηµ′ν′ Γαβ − f µ′ν′ ηαβ) = ∑ µ′ν′ ∑ αβ = ∑ µ′ν′ ∑ αβ = (ΥΥηΓ) µν − (ΥΥηΓ) µν − (ΥΥf η) µν + (ΥΥf η) µν , (7.20) and 3○ = ∑ ij uµ i ( 1 2 (nanb¯nc + ¯na¯nbnc)(1 + Pij)ηciabΓabcj) uν j ∑ abc = = = = 1 2 1 2 1 2 1 2 (nanb¯nc + ¯na¯nbnc) (U α ciU β abU α′ ab U β′ cj uµ i uν j ) (1 + Pij) ηαβΓα′β′ ∑ ij ∑ abc (∑ ab abU α′ U β ab nanb ∑ ijc ciU β′ U α cj uµ i uν j ¯nc − ∑ ab abU α′ U β ab ¯na¯nb ∑ ijc ciU β′ U α cj uµ i uν j nc) × (1 + Pij) ηαβΓα′β′ (Pβα′ αµβ′ν Ξ βα′ − P Ξαµβ′ν + Pβα′ ανβ′µ Ξ βα′ − P Ξανβ′µ) ηαβΓα′β′ ′ (PΞηΓ − PΞηΓ + PΞ ηΓ − PΞ′ηΓ) µν (7.21) with P and P as defined in 7.11. 136 For the two-body flow equation (3.11), dΓijkl ds U α ij = ∑ αβ dΓαβ ds U β kl (1 − Pij)(ηiaΓajkl − fiaηajkl) − (1 − Pkl)(ηakΓijal − fakηijal) = ∑ a + 1 2 (1 − na − nb)(ηijabΓabkl − Γijabηabkl) ∑ ab (na − nb)(1 − Pij)(1 − Pkl)ηaibkΓbjal + ∑ ab dΓαβ ds = ∑ ijkl U α ij ∑ a {(1 − Pij)(ηiaΓajkl − fiaηajkl) − (1 − Pkl)(ηakΓijal − fakηijal)}U β kl ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ 1○ + U α ∑ ijkl (1 − na − nb)(ηijabΓabkl − Γijabηabkl)U β kl 1 2 ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ 2○ ij ∑ ab U α ij ∑ ab (na − nb)(1 − Pij)(1 − Pkl)ηaibkΓbjalU β kl + ∑ ijkl ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ 3○ = (̃ΥηΓ − ̃Υf η − ̃Υ′ηΓ − ̃Υ′f η) αβ + 1 2 (η(P − P)Γ − Γ(P − P)η) αβ + df adf asf a where we have 1○ = ∑ ijkl ∑ a (1 − Pij) (U α ijU β kluµ i uν aU α′ aj U β′ kl ) (ηµνΓα′β′ − f µνηα′β′) − ∑ ijkl ∑ a (1 − Pkl) (U α ijU β kluµ auν kU α′ ij U β′ al ) (ηµνΓα′β′ − f µνηα′β′) = ̃Υνα′µα (ηµνΓα′β − f µνηα′β) − ̃Υµβ′νβ (ηµνΓαβ′ − f µνηαβ′) = (̃ΥηΓ − ̃Υf η − ̃Υ′ηΓ − ̃Υ′f η) αβ 137 (7.22) (7.23) (7.24) (7.25) (7.26) where Υ is defined as ̃Υ = ∑ ijkl ∑ a UijuiuaUaj (7.27) ̃Υ demonstrates marginal dissimilarity compared to the definition of Υ or Υ provided in equation (7.17), where the former encompasses summation across the entirety of the single- particle axis. ̃Υ and ̃Υ′ exhibit permutation-induced variation in indices. 2○ = = = = 1 2 1 2 1 2 1 2 ∑ ijkl U α ij ∑ ab (¯na¯nb − nanb) (U α′ ij U β′ ab U α′′ ab U β′′ kl − U α′′ ij U β′′ ab U α′ ab U β′ kl ) U β klηα′β′ Γα′′β′′ (δαα′ δββ′′ β′α′′ P − δαα′′ δββ′ P β′α′′ − δαα′ δββ′′ Pβ′α′′ + δαα′′ δββ′ Pβ′α′′) ηα′β′ Γα′′β′′ (ηαβ′ (P β′α′′ − Pβ′α′′) Γα′′β − Γαβ′′ (P β′α′′ − Pβ′α′′) ηα′β) (η(P − P)Γ − Γ(P − P)η) αβ (7.28) where the final term involves a straightforward matrix multiplication operation. 3○ = ∑ ijkl ijU β U α kl ∑ ab (na¯nb − ¯nanb) (1 − Pij − Pkl + PijPkl) ∑ α′β′ ∑ α′′β′′ = ∑ ijkl ijU β U α kl (1 − Pij − Pkl + PijPkl) ∑ α′β′ ∑ α′′β′′ ai U β′ U α′ bk U α′′ bj U β′′ al [(U α′ ai U β′′ al na)(U β′ bk U α′′ bj ¯nb) − (U α′ ai U β′′ al ¯na)(U β′ bk U α′′ bj nb)] ηα′β′ Γα′′β′′ ∑ ab = (1 − Pij − Pkl + PijPkl) ∑ α′β′ ∑ α′′β′′ ∑ ijkl ijU β U α kl (Tα′β′′ il T β′α′′ kj − T α′β′′ il Tβ′α′′ kj ) ηα′β′ Γα′′β′′ = (1 − Pij − Pkl + PijPkl) ∑ ijkl ∑ β′β′′ ijU β U α kl ((Tη)β′′β′ il = (1 − Pij − Pkl + PijPkl) ∑ ijkl ijU β U α kl (◻ilkj − ◻kjil) 138 (TΓ)β′′β′ kj − (Tη)β′′β′ il (TΓ)β′′β′ kj ) (7.29) (7.30) where we can choose to retain the indices (i, j, k, l) in the preceding expression at this inter- mediate stage, opting against further contraction. Subsequent contraction would exacerbate the computational intricacies involved. It is worth emphasizing that this particular term holds sifnificant importance in the context of the SVD-IMSRG development, which will be elaborated upon in the subsequent section. Here, we have formally defined (Tη)ββ′(TΓ)ββ′′ , ◻ = ∑ ββ′ Tij = ∑ i UaiUajna , Tij = ∑ i UaiUaj ¯na (7.31) 7.2 Computational Scaling Let us briefly discuss the computational scaling of the SVD-IMSRG(2) approach com- pared with IMSRG(2) before delving into its implementation in the following section. Since the effort is primarily driven by the two-body operators, we focus on those terms in the flow equation. During a single integration step, the two-body contribution to the flow equa- tion (7.13) requires O(D2) compared to O(N 4) in IMSRG(2), where N indicates the size of the single particle basis and D < N 2 — ideally D ≪ N 2 — in rank-reduced space. The computational effort is dominated by the two-body flow equation (3.11), which naively re- quired O(N 6) operations IMSRG(2). In the context of SVD-IMSRG(2), a notable achieve- ment is the successful compression of the ladder term (7.28) from its original complexity of O(D3), again assuming that D ≪ N 2. It is imperative to recognize that the intricacy in the particle-hole term (7.29) remains unaltered, thereby necessitating a comprehensive exploration of alternative contraction methodologies or diverse decomposition techniques to effectively tackle this aspect. Moreover, there arises a demand for the retention of specific matrices or tensors, essential for matrices/tensors precomputed and stored before the flow 139 equations. Notably, the tensorial representation of these flow equations introduces auxiliary tensors, thereby requiring their temporary storage to ensure seamless implementation. 140 Chapter 8 SVD-IMSRG for Schematic Models In the following section, we will explore the application of the SVD-IMSRG to a family of schematic models based on pairing and pair-breaking interactions. We will see that the en- hanced computational scheme can lead to significant speedups compared to the conventional IMSRG approach, provided the pair-breaking interaction is weak. 8.1 The Schematic Models 8.1.1 The Pairing Model The so-called pairing model Hamiltonian describes a quantum system of A fermions in a set of uniformly spaced energy levels, each having a spin degeneracy of two. The Hamiltonian is given by H = δ ∑ pσ (p − 1)a† pσapσ ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ ≡H0 1 2 − g ∑ pq p−aq−aq+ p+a† a† ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ ≡V . (8.1) Here, p is the energy quantum number of a single-particle state, and σ denotes the spin projection. In the one-body operator H0, the constant δ controls the spacing of the energy levels. V is a two-body interaction term with a coupling constant g, which describes the interactions of pairs of particles with the same energy but opposite spin projection. As a concrete example, let us consider A = 4 particles in N = 8 single-particle states. Since 141 Figure 8.1: Illustration of a pairing model configuration with 4 particles in 8 single-particle states. The first two levels with shading represent (occupied) hole states, while the upper two levels depict (unoccupied) particle states. H0 is a diagonal one-body operator, we know that its eigenstates are Slater determinants (cf. Sec. 1.3.9) that correspond to all the possible ways in which we can distribute 4 particles into 8 levels — (8 4) = 70 in total. We can then use these states as a basis for our many-body Hilbert space. In Fig. 8.1, we illustrate the basis configuration corresponding to the ground state of H0, with the lowest two energy levels occupied by pairs. Upon further inspection, we note that H does conserve the total spin projection MS: Since H0 is diagonal, it merely counts whether a given single-particle state ∣p, σ⟩ is occupied or not, and V only acts on pairs with spin projection 0, either counting them (for p = q) or scattering them to different energy (p ≠ q). Thus, the Hamiltonian matrix consists of blocks with constant MS (see Fig. 8.2). From Fig. 8.1 it is apparent that the possible values of MS range from −2 (all 4 particles with in states with spin − 1 2) to +2 (all 4 particles in states with spin + 1 2 ). Since V only acts on pairs of nucleons, as explained above, the number of pairs NP is an additional conserved quantity, so there will be sub-blocks with constant NP within each MS block. We can easily find their dimensions: For MS = ±2, all four particles will be in different energy levels, and there will be no pairs. There is just one possible arrangement (since the fermions are undistinguishable), hence the dimension of these blocks is just 1. For MS = ±1, we have need to have three particles with the same spin projection, and one with 142 𝑝=0𝑝=1𝑝=2𝑝=3𝜎+𝜎−Fermi Level𝛿 Figure 8.2: Illustration of the Hamiltonian matrix for 4 particles in 8 single-particle states. Since the pairing interaction conserves MS, we obtain blocks ranging from MS = −2 to 2. These blocks have additional structure because the number of pairs is conserved as well. When a pair-breaking is added (cf. Sec. 8.1.2), the off-diagonal blocks shaded in gray are populated. the opposite projection. Thus, the total dimension of each of these blocks is ( N+ 3 ) × ( N− 1 ) = ( N+ 1 ) × ( N− 3 ) = 4 × 4 = 16 . Within these blocks, there are 12 states for which there is one pair, i.e., when one spin-up and one spin-down fermion have the same energy, and four states in which the single-particle energies are all distinct. Finally, we need to consider the MS = 0 block, for which we have ( N+ 2 ) × ( N− 2 ) = 36 . 143 𝑀!=+2𝑀!=+1𝑀!=0𝑀!=−1𝑀!=−2𝑁!=0𝑁!=1𝑁!=0𝑁!=1𝑁!=2𝑁!=1𝑁!=0𝑁!=0𝑁!=0 Among these states, there are three possible numbers of pairs: NP = 2, 1, 0. For NP = 2, the spin-up and spin-down fermions must occupy the same energy levels, so there are (4 2) = 6 options in total. For NP = 1, we place one pair in one of the energy levels, and the spin- up and spin-down fermions must be assigned to levels with different energies. There are four possible choices for the pair, three for the first unpaired fermion, and then two for the remaining one, which implies that the dimension of the block is 4 × 3 × 2 = 24. Finally, for NP = 0, there are six configurations in which the particles are completely unpaired. Since the overall dimension of the Hamiltonian matrix in our example is rather small, we can diagonalize H exactly, which will be very useful for benchmarking in the following. We can reduce this effort even further by diagonalizing the blocks separately. Moreover, if we’re only interested in comparing ground-state energies, then for an attractive pairing force (g > 0), we can assume that the ground state will lie in the six-dimensional MS = 0, NP = 2 block of the matrix (cf. Fig. 8.2, and we will only need to diagonalize H0,2 ≡ −g/2 −g/2 −g/2 ⎡ ⎢ 2δ − g −g/2 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ −g/2 −g/2 −g/2 −g/2 0 0 −g/2 −g/2 −g/2 0 4δ − g −g/2 −g/2 0 −g/2 6δ − g 0 −g/2 −g/2 0 6δ − g −g/2 −g/2 −g/2 −g/2 8δ − g −g/2 −g/2 −g/2 −g/2 10δ − g ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (8.2) Note that for weak pairing g ≪ 1, the matrix is heavily dominated by the diagonal, and the diagonal entries are the eigenvalues of the one-body part of H. Thus, the configuration shown in Fig. 8.1 will be the dominant component of the exact ground state. Now let us consider performing an IMSRG(2) calculation for the chosen model. Our 144 first task is to choose a reference state. As argued in previous chapters, it is beneficial to use a reference that is a good approximation to the ground state of the system, so the aforementioned basis configuration will be a suitable choice for ∣Φ⟩. Previously, we have argued for the use of Hartree-Fock optimized Slater determinants as references. In the present case, it turns out that the HF determinant is identical to ∣Φ⟩, and only the single- particle energies get modified: ⟨pσiσ′∣ v ∣pσiσ′⟩ = δ(p − 1) + npσ ⟨pσp − σ∣ v ∣pσp − σ⟩ ϵpσ = δ(p − 1) + ∑ iσ′ g 2 = δ(p − 1) − npσ , (8.3) where we have used the properties of the pairing interaction, as well as npσ = 1 for hole and 0 for particle states. This is not a surprising result: The optimal determinant produced by the HF procedure is decoupled from 1p−1h excitations by construction, but as mentioned above, the interaction term can only excite fermions pairwise. Thus, the MS = 0, NP = 0 block of the Hamiltonian matrix is spanned by ∣Φ⟩ as well as its 2p − 2h and 4p − 4h excitations. In Fig. 8.3, we show the IMSRG(2) flow for the pairing model with 4 particles in 8 levels, using δ = 1.0 and a pairing strength g = 0.5. The curves illustrate how correlations that can be captured by second- and third-order MBPT are absorbed into the flowing Hamiltonian, since the MBPT corrections to E(s) are practically zero for s ≫ 1. The final final ground- state energy E(∞) = 1.4133 is slightly lower than the exact ground-state energy value of 1.4168, indicated by the dashed line, which we obtained by exact diagonalization of the Hamiltonian. This differences is due to the truncation of the IMSRG to two-body operators. 145 Figure 8.3: IMSRG(2) ground-state energy E(s) (i.e., the zero-body part of H(s)) for 4 particles in 8 levels described by a pairing Hamiltonian with δ = 1.0 and pairing strength g = 0.5. The White generator (3.18) was used for the calculation. The figure also displays the impact of the perturbative second (∆E(2)) and third-order energy corrections (∆E(3)) as a function of the flow parameter (s). The dashed line in the plot corresponds to the exact ground-state energy value. Adapted from [41]. 8.1.2 The Pairing plus Pair-Breaking Model As explained in the previous section, the pairing model exhibits symmetries that can be leveraged to simplify calculations, but which ultimately restrict the types of configurations that can mix with the ground-state. While spin and angular momentum are usually con- served by realistic nuclear interactions, the number of pairs is not, although pairing does play an important role in nuclei. To enhance the complexity of our model, we can add a pair-breaking interaction of the form H = Hpair + 1 2 b ( p≠p′ ∑ pp′q p+a† a† p′+aq−aq+ + H.c.) (8.4) 146 10−410−310−210−1100101s1.421.441.461.481.5E[a.u.]E+∆E(2)+∆E(3) where H.c. indicates the Hermitian conjugate of the previous term in the parenthesis — this term would create a pair out of two fermions with opposite spins and different energy quantum numbers. Referring to the Hamiltonian matrix as shown in Fig. 8.2, the pair-breaking interaction will populate the blocks adjacent to the diagonal, which couple the blocks with NP to blocks with NP ± 1. The dimensions are still sufficiently small to allow an exact diagonalization. 8.2 SVD-IMSRG(2) for the Pure Pairing Model Let us now consider the SVD-IMSRG(2) for the pure pairing model. Since the one-body operator’s dimension, N , is significantly smaller than that of the two-body operator, N 2, our main focus is on decomposing constructing a low-rank approximation of the two-body part of the flowing Hamiltonian, Γ, that retains critical information and characteristics while allowing us to substantially reduce the computational cost of the IMSRG evolution. We start by examining the singular values of the initial Gamma (Γ0) and final Gamma (Γ∞) coefficient matrices, which are shown in Fig. 8.4 for our example with N = 8 levels (i.e., Γ can be represented by a 64×64 matrix). Notably, Γ0 exhibits a rank of one, indicating that it contains only one dominant feature — this is as expected, because the pairing interaction (8.1) is characterized by the single constant g. However, when considering the singular values of Γ∞, we see that it has a rank of approximately 20, although most of the singular values are quite small compared to the dominant ones. Thus, the evolution induces contributions to Γ that effectively increase the rank, analogous to what we found in the momentum-space SVD-SRG applications discussed in Chapter 5, although the objects we are considering here are somewhat different. 147 Figure 8.4: Singular-value spectra of the initial Γ0 and the final Γ∞ for a pairing model with A = 4, N = 8. Note that Γ0 matrix essentially consists of a single nonzero singular value, while the evolved Γ∞ has a higher rank of 20. Clearly, it would be a bad idea to use the one dominant singular vector of Γ0 as a projector for the SVD-IMSRG flow equation, and it is impossible to anticipate which of the other remaining singular vectors that are associated with vanishing singular values may become important during the evolution. Naturally, we also want to avoid using the singular vectors of Γ∞ as projectors, because this would force us to perform the evolution first, and defeat the purpose of the SVD-IMSRG. However, we do have quantities available at s = 0 that may indicate how the system is going to evolve, namely the generator η and the derivative dΓ ds ∣ s=0 matrix — we’ll refer to it as dΓ for the sake of brevity in the following. It turns out that the singular vectors of dΓ are capable of capturing all the essential features of the Γ∞, as illustrated in Figure 8.5b. In panel (a), we show Γ∞ in the singular vector basis of Γ0. While the resulting matrix is sparse and structured, non-zero matrix elements appear in its far off-diagonal sectors. In contrast, panel (b) shows Γ∞ in the basis of singular vectors of dΓ, and we find that the matrix elements are confined to a 20 × 20 block in the upper left corner that suggests that 148 0102030405060rank101610131010107104101iSingular values of 0 (a) (b) (c) Figure 8.5: These graphs showcase the transformation of Γ matrices into singular-vector spaces derived from various sources. Panel (a): Representation of Γ∞ in the singular-vector space of Γ0. Panel (b): Representation of Γ∞ in the singular-vector space of dΓ. Panel (c): Representation of Γ∞ in the singular-vector space of dΓ with basis permuted. the evolution can be captured in a space of dimension 20 instead of 64. Moreover, there are only two matrix elements located at the off-diagonal corners of the matrix in this proposed effective space. It turns out that the “emergence” of these off-diagonal matrix elements can be anticipated based on the structure of the generator η at s = 0, which is shown in the singular-vector basis of dΓ in Fig. 8.6. We compute the matrix elements ⟨ui∣ η(s = 0) ∣uj⟩ and permute the basis vectors such that the non-zero matrix elements are all concentrated in a single corner, as shown in Fig. 8.6b. In our specific example, we obtain a 4 × 4 block. In the rearranged basis, Γ∞ also has a prominent 4 × 4 block in the upper left corner, and it is diagonal in the space spanned by the 16 vectors that remain from the initial truncation to dimension 20. We will see below that this diagonality is a particular feature of the pure pairing model, but that characteristic block structures within the singular space will also appear for the pairing plus pair-breaking model (Sec. 8.3), as well as for realistic nuclear interactions (Sec. 9.2). 149 010203040506001020304050600102030405060010203040506001020304050600102030405060 (a) (b) Figure 8.6: Representation of η in the singular-vector basis of dΓ before (panel a) and after permuting the basis (panel b). Let us now consider the ground-state energy of our model system in the SVD-IMSRG(2) as a function of the rank, i.e., the number of singular vectors of dΓ that are used to project the flow equations after performing the basis rearrangement. We find that even with just 4 components, the SVD-IMSRG(2) calculation produces highly accurate results, yielding an error of approximately 0.1%. Given the prominence of the 4 × 4 block in Figs. 8.5c and 8.6b, it is perhaps not surprising that it is sufficient to capture the essential features of our model. As we include more components, the error diminishes exponentially, and once we include the 20 relevant singular vectors of dΓ, we recover the IMSRG(2) ground-state energy exactly. Including additional singular vectors in the basis does not lead to further improvement, as shown in Figure 8.7a. We can think of r = 20 as the effective rank of Γ∞, or our overall model. It corresponds to about 30% of the full dimension of the two-body interaction. In realistic calculations, we also need to explore the convergence of the IMSRG(2) with respect to truncation of the single-particle basis truncation, hence we also compare SVD- IMSRG(2) and IMSRG(2) for a pairing model with 4 particles in 20 single-particle states in Fig. 8.7b. We find that the overall behavior of the SVD-IMSRG(2) and IMSRG(2) ground- state energies as a function of the rank is very similar to the N = 8 case. We start from 150 0204060010203040506002040600102030405060 (a) (b) Figure 8.7: Ratio of SVD-IMSRG(2) and IMSRG(2) ground-state energies as a function of the truncation rank for pairing models with 4 particles in 8 single-particle states (panel a) and 20 single-particle states (panel b), respectively. a slightly higher relative error of 0.3% when we include only 4 components, which is still highly accurate. We observe the same steep drop in the relative error of the g.s. energy as we go from 4 to 8 components, which means that this feature must be entirely driven by the number of particles, or equivalently, the number of hole states. The exponential tail of the relative error decays a little bit more slowly than for 8 single-particle states. The SVD-IMSRG(2) and IMSRG(2) ground-state energies become identical once we include 69 components. Compared to the dimension 400 of the two-body basis, this corresponds to an effective rank of 17%. Overall, our present findings suggest that the IMSRG(2) can enable substantial computa- tional savings, especially as we approach realistic basis sizes where the number of hole states Nh, which is identical to the particle number, is much smaller than the number of particle states Np and the overall basis size, Nh ≪ Np ∼ N . 151 4812162024rank1.00001.00021.00041.00061.00081.00101.00121.0014Er/EfullGround-state energy as a function of rankSVD-IMSRGIMSRG4812162024rank1.00001.00051.00101.00151.00201.00251.0030Er/EfullGround-state energy as a function of rankSVD-IMSRGIMSRG Figure 8.8: Singular-value spectra of the initial Γ0 and the final Γ∞ for a pairing plus pair- breaking model with A = 4, N = 8. In contrast to Fig. 8.4, the initial Γ0 has two features, which we associate with the coupling constants of the pairing and pair-breaking interactions. The evolved Γ∞, has a higher rank of 28, but we note that the singular values still rather quickly drop off by orders of magnitude. 8.3 SVD-IMSRG(2) for a Pairing plus Pair-Breaking Model Let us now enhance the model by adding a pair-breaking term to the Hamiltonian (cf. Sec. 8.1.2) and analyze the singular values of the Γ0 and Γ∞ once more. As we can see in 8.8, Γ0 exhibits a rank of 2, which signifies the presence of only two dominant features, the strength of the pairing and pair-breaking interactions, respectively. After IMSRG evolution, Γ∞ has a higher rank of 28. This finding once again emphasizes that the singular vectors of Γ0 are a poor choice for capturing the features that are induced by the IMSRG evolution in a compact fashion. This is further illustrated by Fig. 8.9a, which shows that Γ∞ is dense in the basis of Γ0 singular vectors. Figure 8.9b reveals that the singular vectors of dΓ allow a more compact representation for Γ∞. We obtain a 28 × 28 block, with dimensions matching the rank of Γ∞ that was 152 0102030405060rank10151012109106103100iSingular values of 0 (a) (b) (c) (d) Figure 8.9: Representation of Γ0 and Γ∞ in bases of select singular vectors. Panel a: Γ∞ represented in the singular-vector basis of Γ0. Panel b: Γ∞ in the singular-vector basis of dΓ, without basis rearrangement. Panels c,d: Γ0 and Γ∞ in the singular-vector basis of dΓ after basis permutation. 153 01020304050600102030405060010203040506001020304050600102030405060010203040506001020304050600102030405060 (a) (b) Figure 8.10: Left panel: Singular values of Γ∞ associated with the two blocks identified in Fig. 8.9. These results were obtained for b = 0.1. Right panel: SVD-IMSRG(2) ground-state energy using the primary block only, compared to the full IMSRG(2) ground-state energy, as a function of the pair-break strength b. revealed by the SVD (cf. Fig. 8.8). By implementing our basis rearrangement procedure from the previous section, we find two prominent blocks: A 10 × 10 block that is present in both Γ0 (Fig. 8.9c) and Γ∞ (Fig. 8.9d), and an 18 × 18 block that is induced during the IMSRG evolution. The greater dimension of the “primary” block and the expansion of the “secondary” block from a diagonal to a dense matrix can be attributed to the inclusion of the pair-breaking interaction in the Hamiltonian. In Fig. 8.10a, we show the singular values associated with the primary and secondary blocks of Γ∞ (also cf. Fig. 8.8). The singular values of the primary block constitute approximately 84% of the trace of Γ∞, those of the secondary block the remaining 16%. Naively, we might expect this to give us an estimate for the error in the ground-state energy if we omit the secondary block altogether. However, Fig. 8.10b reveals that the impact on the ground-state energy is significantly smaller. For a pair-breaking strength b = 0.1, for example, the difference between the IMSRG(2) and SVD-IMSRG(2) with the primary block alone is a mere 0.1%. This is roughly comparable to what we found when evolving the pure 154 0510152025rank102101100iSingular values of 1st block2nd block0.100.150.200.250.30b1.241.261.281.301.321.341.361.381.40EgsGround-state energy as a function of pair-breaking strengthIMSRGSVD-IMSRG pairing model with only 4 components, and we note that there are three singular values in Fig. 8.10a that are significantly larger than the others. Given that the majority of the elements of the secondary block are tied to the presence of the pair-breaking interaction in our model Hamiltonian, it is not surprising that the error from its omission grows with the strength of the pair-breaking interaction. However, increasing the pair-breaking strength three-fold to b = 0.3, we can still reproduce the IMSRG(2) ground-state energy within about 1% at a fraction of the computational cost. 155 Chapter 9 SVD-IMSRG for Infinite Matter 9.1 Preliminaries In this section, we briefly recapitulate the construction of single-particle basis and the representation of two-nucleon interactions in infinite matter calculations, following Ref. [64]. 9.1.1 Single-Particle Basis for Infinite Matter The description of infinite nuclear matter relies on plane-wave functions that are confined to a box with volume Ω and side length L, which quantizes the spectrum. (Eventually, one needs to take the limit L → ∞ in computed expectation values.) The wave function is given by Ψkσ(r) = 1 √ 2Ω exp(ikr)ξσ (9.1) where k denotes the wave number, and ξσ are basic spinors for nucleons with either spin up or spin down: ξσ=+1/2 = ( 1 ) 0 ξσ=−1/2 = ( 0 ) . 1 (9.2) We can think of an infinitely extended system in terms of periodic copies of our “elemen- tary” box, which implies the use of periodic boundary conditions on the single-particle wave functions. These conditions impose constraints on the allowed wave numbers, which can be 156 represented as kni = 2πni L i = x, y, z ni = 0, ±1, ±2, ⋯ (9.3) The kinetic energy operator can then be expressed as ̵h2k2 p 2m T = ∑ pσp a† pσpapσp (9.4) where a† pσp and a† pσpapσp are the creation and annihilation operators, respectively, associated with the momentum state represented by the wave number p. When applying periodic boundary conditions, the discretization of the single-particle momenta lead to the following expression for the single-particle energy: εnx,ny,nz = ̵h2 2m (k2 nx + k2 ny + k2 nz ) = ̵h2 2m ( 2π L 2 ) (n2 x + n2 y + n2 z) . (9.5) To avoid singling out specific spatial directions, we impose truncations on the single- particle basis by restricting the single-particle energy εnx,ny,nz ≤ εmax (or, equivalently, the total single-particle momenta or quantum numbers ni). It is clear from Eq. (9.5) that for each ε allowed by the truncation, we may have several energetically degenerate realizations through different combinations of the quantum numbers nx, ny, nz. Thus, the single-particle basis will have a structure analogous to a shell model. In our calculations, we ensure a “closed-shell” structure for both occupied and unoccupied single-particle states by carefully including all single-particle states with energies below our predetermined cutoff. Consider, for example, the case (n2 x + n2 y + n2 z) ≤ 3 . (9.6) 157 y + n2 y + n2 x + n2 n2 0 1 1 1 1 1 1 2 2 2 2 2 2 2 z nx ny nz N(cid:153) 0 2 0 0 -1 1 0 0 -1 1 -1 1 0 0 0 0 -1 1 0 0 0 0 -1 -1 1 1 -1 -1 1 0 0 0 0 0 -1 1 0 0 0 0 -1 1 -1 14 x + n2 n2 2 2 2 2 2 3 3 3 3 3 3 3 3 38 z nx ny nz N(cid:153) 0 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 1 0 0 0 0 -1 -1 -1 -1 1 1 1 1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 54 Table 9.1: Total number of particles N(cid:153) for various n2 fermion species, accounting for a spin degeneracy of 2 for each energy state. x + n2 y + n2 z values for one spin-1/2 The resulting list of energy levels is presented in Table 9.1, considering one species of identical fermions and accounting for both spin-up and spin-down solutions. If we were to consider n2 x + n2 y + n2 z = 4, we obtain 6 additional energy levels (with one ni = ±2 and the others vanishing) and 12 additional states due to spin degeneracy, resulting in a new magic number of 66. Proceeding in this manner, one obtains the magic numbers of 2, 14, 38, 54, 66, 114, . . .. For the study of infinite nuclear matter with both protons and neutrons, the magic numbers become 4, 28, 76, 108, 132, 228, and so on. Once the number of particles in the simulation is determined, and a specific density ρ is chosen, the Fermi momentum kF of the system is determined by the equation ρ = g k3 F 6π2 (9.7) where g represents the degeneracy. For a system of one type of spin-1/2 particles, the degeneracy is 2, while for nuclear matter with protons and neutrons, it is 4. 158 Since we are typically interested in the energy of neutron or nuclear matter at a fixed density, we can use ρ and the particle number A that we choose according to our “closed-shell criterion” to define the size of the simulation box via Ω = L3 = A ρ . (9.8) 9.1.2 Two-Nucleon Operators Next, we discuss the form of two-body operators in the plane-wave basis defined in the previous section. The matrix elements of the initial two-nucleon interaction have the form ⟨kpσpτpkqσqτq∣ ˆv ∣krσrτrksσsτs⟩ , (9.9) where the quantum numbers are the components of the single-particle momenta kp, the spin σp, and the isospin τp. The interaction conserves the spin projection MS = σp+σq = σr +σs and isospin projection (or charge) MT = τp + τq = τr + τs. Then MT = 0 is the pn (proton-neutron) channel, MT = +1 the pp (proton-proton) channel, and MT = −1 the nn (neutron-neutron) channel (cf. Chapters 2 and 5). Furthermore, translational and Galilean invariance imply that it must be diagonal in the center-of-mass momentum K = kp + kq = kp + kq and the matrix elements may not explicitly depend on K. Switching to plane waves in K and the relative momentum k = 1 2 (kp − kq) (9.10) 159 we have ⟨kpσpτpkqσqτq∣ˆv ∣krσrτrksσsτs⟩ = ⟨Kkσpτpσqτq∣ ˆv ∣K′k′σrτrσsτs⟩ = δMS ,M ′ S δMT M ′ T δ(K − K′) ⟨kMSMT ∣ ˆv ∣k′MSMT ⟩ . (9.11) (9.12) We can couple the spins of the individual nucleons to total spin S and introduce ⟨kMSMT ∣ˆv(k, k′) ∣k′MSMT ⟩ 1 2 σp 1 2 ⟨ = ∑ S σq∣SMS⟩ ⟨ 1 2 σr 1 2 σs∣SMS⟩ ⟨kSMSMT ∣ ˆv(k, k′) ∣kSMSMT ⟩ , (9.13) where ⟨ 1 2σp 1 2σq∣SMS⟩ and ⟨ 1 2σr 1 2σs∣SMS⟩ are Clebsch-Gordan coefficients. The two-nucleon interaction is frequently analyzed in a partial wave representation, be- cause it is symmetric under rotations and parity transformations. A plane wave state with linear Cartesian momentum k, irrespective of whether it characterizes a single nucleon or the center-of-mass or relative momenta of a pair of nucleons, can be expanded in partial wave momentum states as ∣k⟩ = 4π ∞ ∑ l=0 l ∑ m=−l ilY (ˆk) ∣klm⟩ , (9.14) where Ylm are spherical harmonics, and ˆk is a unit vector characterizing the angular coordi- nates in momentum space. Then ⟨kSMST MT ∣ ˆv ∣k′SMST MT ⟩ = (4π2) ∑ JM ∑ lm ∑ l′m′ il′−lY ∗ lm(ˆk)Ylm(ˆk′) ⟨lmSMS∣JM ⟩ ⟨l′m′SMS∣JM ⟩ ⟨klST J∣ ˆv ∣k′l′ST J⟩ , (9.15) 160 Figure 9.1: Plane-wave matrix representation of Γ0 derived from an NNLO chiral interaction. We used an Nmax = 4 basis consisting of 66 single-particle states, and A = 38 particles. The block structure results from the conservation of K, MS and MT . Matrix elements were provided by O. Udiani and K. Yu. where ⟨klST J∣ ˆv ∣k′l′ST J⟩ are the partial-wave matrix elements used throughout chapters 2 and 5. General two-body operators we must handle during the IMSRG and SVD-IMSRG evo- lution of an infinite-matter system will have a structure that is similar to what we dis- cussed for the two-nucleon interaction, with one important difference: When we normal order the infinite-matter Hamiltonian with respect to some single-particle momentum distri- bution np(k),the center-of-mass momentum K will still be conserved, but Γ0 will in general become dependent on its value. Thus, Γ0 will assume a block-diagonal structure, as illus- trated in Fig. 9.1. Other relevant operators like η, dΓ, and Ω will inherit this structure as well, hence it can be leveraged in the implementation of the IMSRG evolution. 161 0500100015002000250030003500400005001000150020002500300035004000 9.2 Towards SVD-IMSRG(2) for Infinite Neutron Mat- ter In this section, we analyze the singular-value spectra of Γ0, Γ∞ as well as the Magnus operator Ω (cf. Chapter 3) in a similar way as those of the schematic models in the previous chapter. The development of an SVD-IMSRG(2) code for infinite-matter problems is still in progress, primarily because of the challenge of finding an efficient treatment of the particle- hole terms in the IMSRG(2) flow equation in a suitably chosen singular-vector basis, as discussed in Chapter 7. Since the treatment of particle-hole type correlations in infinite-matter calculations is difficult and they have been estimated to give only small contributions to the neutron matter and nuclear matter equations of state, they have historically been omitted in many beyond- mean field approaches, including the IMSRG(2) calculations in Ref. [10]. In the case of the SVD-IMSRG(2), this omission allows the flow to be implemented at O(r3) cost, with r ≪ N 2, yielding order-of-magnitude speedups, especially for the high-dimensional plane- wave single-particle bases that need to be considered here. However, the validity of this omission deserves ongoing scrutiny. As a concrete example, we consider pure neutron matter with A = 38 neutrons in an Nmax = 4 basis consisting of 66 states — i.e., the setup is rather similar to our schematic model calculations in the previous chapter. Since the matrix dimension of Γ and other two- body matrices is aready substantial at 662 = 4356, testing convergence with respect to the single-particle basis truncation is numerically challenging. As our interaction, we choose a chiral two-nucleon force at next-to-next-to-leading order (NNLO), as briefly introduced in Chapter 1. 162 Figure 9.2: Singular-value spectra of Γ0, Γ∞ and dΓ for an NNLO chiral interaction. We used an Nmax = 4 basis with N = 66 and A = 38, as in Fig. 9.1. In Fig. 9.2, we show the singular-value spectra of Γ0, Γ∞ and dΓ. Contrary to what we found for the schematic models in Chapter 8, the distribution of singular values for Γ0 and Γ∞ is very similar, and very extended — there is a quick dropoff after just a few components, but there are more than 2000 singular values that fall in a range covering of 1-10% of the largest ones, and it is difficult to tell from this figure which of these operators may be relevant for the evolution or not. The similarity of the singular values actually suggests that many components of the interaction are not significantly evolving at all, and we should be subtracting of a constant “background” part of Γ — this would be analogous to our approach in the free-space SVD-SRG. The structure of dΓ’s singular-value spectrum supports this, because it exhibits a much more pronounced drop off. However, this raises the question whether the projection of that hypothetical background part into the truncated singular-vector basis of dΓ will yield a sufficiently accurate evolution. Let us now proceed as in our analysis of schematic models and project Γ∞ into the (re- arranged) basis of singular vectors of dΓ. The resulting matrix is shown in Fig. 9.3a. We recognize the same two-block structure, and in the secondary block a pronounced diagonal. 163 01000200030004000101710141011108105102101Singular Value Spectrum0d (a) (b) (c) (d) Figure 9.3: Representation of Γ∞ (panel a) and Ω∞ (panel b) in the singular-vector basis of dΓ. Panels c and d show the singular value spectra associated with the primary and secondary blocks. 164 0500100015002000250030003500400005001000150020002500300035004000050010001500200025003000350040000500100015002000250030003500400001000200030004000101710141011108105102101Singular value spectrum of Singular Spectrum for 1st blockSingular Spectrum for 2nd block010002000300040001021101810151012109106103100Singular value spectrum of Singular Spectrum for 1st blockSingular Spectrum for 2nd block Figure 9.3c attributes the singular values to the primary and secondary blocks. The cumula- tive sum of the singular values in Γ∞’s primary block yields about 76% of the cumulative sum of the singular values, somewhat similar to our findings for the pairing plus pair-breaking Hamiltonian. In that case, we were able to omit the secondary block while only causing a small loss of accuracy compared to the full IMSRG(2) flow. For the schematic model, we also saw that the secondary block is strongly tied to the particle-hole term in the flow equa- tions, this could support the idea that the associated correlations are not very important in infinite-matter calculations, as discussed earlier. Instead of considering Γ∞, we can also consider the IMSRG evolution from the perspective of the Magnus approach discussed in Sec. 3.5. In Fig. 9.3b, we show the representation of Ω∞ in the basis of singular vectors of dΓ, and Fig. 9.3d is the associated singular-value spectrum. We see that Ω∞ has a similar block structure as Γ∞, but importantly, it turns out that the primary block contains 99.7% of the cumulative sum of singular values. This suggests that we may be more successful in using the SVD to compress the flow equations for Ω (cf. Sec. 3.5) than those for H. While these flow equations also include particle-hole terms, the present result suggests that it may be much safer to omit them in the evolution of Ω than in the evolution of H. In our present example, this would mean evolving a block whose dimension is approximately 500 instead of the the size of the overall dimension 4356. Moreover, with Ω∞ at hand, we can construct not just H∞ but general transformed observables, so this approach may be preferable over direct evolution of H, anyway. This is a key direction for future research. 165 Chapter 10 Further Developments 10.1 Exploiting Block Structures in the IMSRG Evo- lution In the previous chapters, we have repeatedly considered omitting parts of the Hamiltonian that may not change significantly throughout the IMSRG flow, or contribute to the evolution of the parts of the Hamiltonian that are changing. This has motivated a more detailed investigation of the (SVD-)IMSRG evolution of specific blocks of H, and potential efficiency gains that could be realized by leveraging these structures more explicitly. 10.1.1 General Considerations In Fig. 10.1, we present a schematic illustration of the structure of Γ. In essence, we need to consider six sectors of the matrix: hhhh, hhhp, hhpp, hphp, hppp, and pppp. The other sectors are related to these either because of the antisymmetry of the two-body basis, or Hermiticity. In Table 10.1, we summarize the time complexities associated with the two-body flow equations for the six essential blocks. The flow equation (3.11) contains three classes of terms: the interplay between one-body and two-body components of η and H (1B-2B), a ladder term and a particle-hole term. We assume generators with a ph / pphh structure, which 166 Figure 10.1: Schematic illustration of the block structure of the Γ matrix. In realistic applications, we typically have Nh ≪ Np ∼ N , hence the pppp block is by far the largest block of the matrix. encompasses the White, White-ArcTan and imaginary time ans¨atze, the Wegner generator has a more complex structure (cf. Sec. 3.4). For example, the 1B-2B term contribution to the pppp block is given by dΓ(1) abcd ds = ∑ p (ηapΓpbcd + ηbpΓapcd − Γabpdηpc − Γabcpηpd) (fapηpbcd + fbpηapcd − ηabpdfpc − ηabcpfpd) − ∑ p (ηaiΓibcd + ηbiΓaicd − Γabidηic − Γabciηid) , (10.1) = ∑ i where we have used the single-particle index conventions of Chapter 3. As we can see, the overall time complexity of the IMSRG(2) is dominated by evaluating the particle-hole term in the evolution of the hppp block, and the ladder term in the pppp and hhpp sectors, which both scale as O(N 2 hN 4 p ). 167 hhhpphpphhhpphpp index type 1B-2B ladder term p-h term 1 2 3 4 5 6 hhhh O(N 4 hNp) O(N 4 p ) hN 2 None hhhp O(N 3 hN 2 p ) O(N 3 hN 3 p ) O(N 4 p ) hN 2 hhpp O(N 2 hN 3 p ) O(N 2 hN 4 p ) O(N 3 p ) hN 3 hphp O(N 2 p ) hN 3 None O(N 3 p ) hN 3 hppp O(NhN 4 p ) O(N 3 hN 3 p ) O(N 2 p ) hN 4 pppp O(NhN 4 p ) O(N 2 p ) hN 4 None Table 10.1: Time complexity analysis of each term in the evolution of the 6 essential blocks. All estimates are based on the generators with an explicit ph / pphh structure (e.g., White, White-ArcTan, and imaginary time). 10.1.2 Numerical Applications Let us now leverage the block structures we identified in the previous section in numerical applications. We consider a schematic Hamiltonian that consists of a pairing plus pair- breaking interaction (cf. Sec. 8.1.2, but we enhance the complexity by adding an additional contribution ijkl = Rijkl exp−(fii+fjj −fkk−fll)2/λ2 ΓR (10.2) where R is random but Hermitized and antisymmetrized, and the Gaussian factor imposes a smooth cutoff not unlike what is found in nuclear interactions. In the right panel of Fig. 10.2 we illustrate the computational time required for perform- ing the full and blocked IMSRG(2) (referred to as Block-IMSRG(2) henceforth) for increasing single-particle basis size N , while Nh = A = 4 is held fixed. We see that when the ratio of particle and hole states, Np/Nh, reaches approximately 10, which starts to reflect realistic applications, the computing time for a Block-IMSRG evolution is reduced tenfold. For larger bases, the savings will be even more substantial. Remarkably, this acceleration is achieved 168 Figure 10.2: Right panel: Computational time for IMSRG(2) and Block-IMSRG(2) as a function of the single-particle basis size N . Left panel, top: IMSRG(2) and Block-IMSRG(2) ground-state energies, along with an assessment of omitted contributions. Left panel, bottom: Relative error in the ground-state energy due to the omission of specific contributions (see text). purely through a more streamlined manipulation of reduced sets of matrix elements in the RHS of the flow equations. In the left panels of Fig. 10.2, we present ground-state energies obtained through IM- SRG(2), Block-IMSRG(2), and a modified Block-IMSRG(2) approach that omits specific contributions in the flow equations: In the context of hhhp and hppp blocks, the O(N 2 p ) hN 4 particle-hole term is omitted, while for pppp, we disregard the interplay between one-body and two-body components, which scales as O(NhN 4 p ) for most of the common generators but O(N 5 p ) for the Wegner generator. The omission of the particle-hole term is considered with an eye on merging block and SVD techniques, because of the implementation difficulties described in previous chapters. The IMSRG(2) flow equation will still scale as O(N 2 p ) hN 4 due to the ladder term in the pppp and hhpp blocks, but these contributions are readily 169 0.40.60.81.01.21.4E1020304050# of s.p. states0.0000.0050.0100.0150.0200.0250.0300.0350.040E/EFullBlockBlock, no hhhp/hpppBlock, no hhhp/hppp/pppp1020304050# of s.p. states100101102103TFullBlock compressible in an SVD basis. The top left panel of Fig. 10.2 shows that the effect of omitting the aforementioned blocks is difficult to discern. In the bottom left panel, we consider the relative error in the ground- state energy instead. We see that as the number of single-particle states N (or, equivalently, Np) increases, the influence of the omitted contributions becomes increasingly pronounced. As we reach Np/Nh ≈ 10, as in the runtime comparison, the error due to the omission of the 1B-2B coupling in the pppp block causes an error of a few percent, while the error due to the omission of the particle-hole term is still merely 0.5%. Next, we consider how the relative errors in the ground-state energy depend on the model parameters. In Figure 10.3, we present a comprehensive analysis of the relative ground-state energy errors that stem from the selective omission of the particle-hole term in blocks 2–5, and the 1B-2B term in the pppp block (#6) — we refer to Tab. 10.1 for the labeling. Since the effects of the particle-hole term are almost negligibile for a pure pairing interaction, we fix the level spacing δ = 1.0 and the pairing strength g = 0.5 of our model Hamiltonian (cf. Chapter 8), but vary the strengths of the pair-breaking interaction b and the random perturbation k. For larger basis sizes, we encountered stability issues when the pair-breaking interaction was increased beyond b ≈ 0.16, hence the range of b is more limited in the center and right columns of the figure. In general, we see that the errors due to the omission of the discussed terms increase with b and k, which is not surprising. The most promising finding is that it seems to be safe to omit the particle-hole term in the hpppp block (#5) without compromising accuracy: For most of the parameter sets studied here, the resulting error in the IMSRG(2) ground-state energy is below 0.1%. There is a jump in that error for the N = 16 basis with pair-breaking strengths b = 0.14 and 0.16, but this is likely due to the onset of instability for the IMSRG(2) 170 Figure 10.3: Relative error in the IMSRG(2) ground-state energy due to the omission of the particle-hole term in the six unique channels (cf. Tab. 10.1) as function of the pair-breaking strength b and the strength of the perturbation k (g = 0.5, δ = 1). The particle number is fixed at A = 4, while the single particle basis size is N = 8 (left), 16 (center) and 30 (right column). Note the more restricted range of the pair-breaking strength in calculations with larger single-particle basis size — this is due convergence problems with the IMSRG(2) flow for greater values of b. flow (as evidenced by the simultaneous growth of all the errors.) Importantly, this is the single particle-hole term that scales as O(N 2 hN 4 p ), in the other blocks, the scaling is milder (cf. Tab. 10.1). While it is not time-critical for most generator choices, we also see that the error due to omitting the 1B-2B term in the pppp channel is in the sub-percent range in most cases. The relative importance does seem to increase in the larger single-particle bases, but as the absolute energy differences shown in Fig. 10.4 reveal, this is owed more to the decrease of 171 105104103102k=1/2023452/52/3/52/4/52/3/4/562/5/6105104103102k=1/1023452/52/3/52/4/52/3/4/562/5/6105104103102k=1/5104103102k=1/2023452/52/3/52/4/52/3/4/562/5/6104103102101k=1/1023452/52/3/52/4/52/3/4/562/5/6105104103102101k=1/523452/52/3/52/4/52/3/4/562/5/6105104103102101k=1/2023452/52/3/52/4/52/3/4/562/5/6104103102101k=1/1023452/52/3/52/4/52/3/4/562/5/6104103102101k=1/5excluded blocksE/Eb=0.1b=0.2b=0.3b=0.4b=0.10b=0.12b=0.14b=0.16b=0.09b=0.11b=0.13b=0.15 Figure 10.4: Absolute error in the IMSRG(2) ground-state energy due to the omission of the particle-hole term in the six unique channels (cf. Tab. 10.1) for δ = 1, g = 0.5, b = 0.1 and k = 1/20. The particle number is fixed at A = 4, while the single particle basis size is varied from N = 8 to 30 (Also see Fig. 10.3). the total ground-state energy than to the absolute magnitude of this contribution. Because of the inclusion of the random perturbation in our Hamiltonian, there is reason to expect that our present findings can be generalized to realistic Hamiltonians, hence block techniques and the strategic omission of the particle-hole terms in certain blocks should be explored in IMSRG calculations with chiral nuclear interactions in future work. This is particularly true because of the difficulties with implementing SVD-based rank reductions for these terms that were discussed in chapters 7 through 9. 10.2 Low-Rank Structures in Magnus-IMSRG Calcu- lations for Finite Nuclei In Sec. 9.2, we found evidence of potentially exploitable low-rank structures in the operator Ω that defines the unitary evolution of infinite neutron matter in the Magnus 172 23452/52/3/52/4/52/3/4/562/5/6excluded blocks0.00000.00250.00500.00750.01000.01250.01500.01750.0200E8 s.p. states16 s.p. states30 s.p. states (a) (b) (c) Figure 10.5: Singular-value spectra of selected (J, π, MT ) channels of the Magnus operator Ω at various stages of the MAGNUS(2) evolution for 16O. Calculations were performed with the EM 1.8/2.0 interaction in a single-particle basis with emax = 10, E3max = 14 and ̵hω = 20 MeV (see Sec. 5.5). formulation of the IMSRG. It turns out that such structures also exist in finite nuclei. As an example, we consider MAGNUS(2) calculations for 16O based on the chiral EM1.8/2.0 two- plus three-nucleon interaction. In Fig. 10.5, we show the singular-value spectra of a few selected (J, π, MT ) channels of the flowing Magnus operator. Even with a still modest basis size of emax = 10, a large number of channels of the operator have dimensions in excess of 10000, hence we employ our RSVD implementation (cf. Sec. 4.2 and Chapter 6). We see that the singular-value spectrum of the Magnus operator is very sharply peaked. 173 The selected partial waves are representative for all the channels of the operator, and likewise, the results for other closed-shell nuclei are similar to those for 16O — there is merely some variation in the number of components due to the changing number of hole states Nh, which will be identical to either the proton or neutron number of the nucleus. Since Ω(s = 0) = 0, it is not surprising that the rank of the operator grows during the Magnus-IMSRG evolution. Once we have evolved past the initial stages of the flow, e.g., around s = 1 for a White generator, the change in singular values happens only in the tail of the singular-value distribution, and the induced singular values in the tail are still surpressed by multiple orders of magnitude. This probably explains the results of Hoppe et al. in Ref. [48], who successfully applied ad-hoc importance measures to reduce the computational effort for MAGNUS(2) calculations. While one could think that it is a problem that Ω(s = 0) = 0, which effectively prevents us from computing its SVD, our success with the bases constructed from η and dΓ suggests that this issue can be overcome. The generator seems a particular promising candidate for defining a singular-vector basis, since it is the leading contribution to the RHS of the Magnus flow equation, dΩ ds = η − 1 2 [Ω, η] + 1 12 [Ω, [Ω, η]] + . . . (10.3) (cf. Sec. 3.5). As in the case of infinite-matter wcalculations, the investigation of Magnus operator decompositions is an important direction for future research. 174 Part IV Conclusions 175 Chapter 11 Conclusions and Outlook In this dissertaion, we have explored the application of mathematical factorization tech- niques — in particular, singular value decompositions (SVDs) — to nuclear interactions and many-body methods, with the primary objectivesof enhancing the computational efficiency and reducing memory demands in state-of-the-art nuclear many-body calculations. We suc- cessfully used the SVD to construct precise low-rank models for two-nucleon interactions, and developed a new implementation of the Similarity Renormalization Group (SRG), a widely used RG method in nuclear theory, to evolve these interactions directly in the SVD-factorized form — this work was published in Ref. [113]. Next, we sought to extend these investigation to three-nucleon interactions, where there is arguably a greater need for model order reduction because of the sized of the matrices that need to be manipulated. We employed randomized SVD algorithms to address the challenges posed by significantly larger basis dimensions. While evidence of low-rank structures in three- nucleon interactions emerges, we identify impediments to sustained rank reduction, which are primarily the result of the interplay of SRG-evolving two- and three-nucleon forces in a three-body basis. Furthermore, we explored SVD-based model order reduction within the framework of the In-Medium SRG (IMSRG) approach for many-body calculations. We derived a represen- tation of the IMSRG flow equations in the singular-vector basis of general operators, and 176 identified the derivative dH/ds and the IMSRG generator η at the initial stages of the flow as suitable sources for generating the singular vectors. However, we also identified persistent computational hurdles linked to the treatment of the particle-hole terms of the IMSRG flow equation in singular-vector basis representations. We began implementing these model-order reduction techniques in a variety of systems, from schematic models of growing complexity to infinite matter and finite nuclei in calcu- lations based on realistic interactions from chiral Effective Field Theory (EFT). While we found emergent low-rank structures in the Hamiltonian and/or the Magnus operator that can be used to parameterize the IMSRG transformation, leveraging them for speedup is made difficult by the aforementioned problems with the handling of the particle-hole term. This challenge has motivated a more detailed investigation of the scaling and significance of all contributions within the IMSRG flow. For Hamiltonians that consist of a schematic pairing plus pair-breaking interaction as well as a randomized interaction with smooth cutoff but tunable strength, we found evidence that the particle-hole term can be safely neglected in several channels of the IMSRG evolution, which would not only eliminate one of the domi- nant contributions to the methods computational scaling, but also pave the way for further computational savings through the application of SVD-based model-order reduction to the remaining terms. In summary, this dissertation advances our understanding of the potential utility of fac- torization techniques in the context of nuclear many-body calculations, offering valuable insights into the potential for model-order reduction and the computational savings it may unlock. It identifies both promising directions for future research as well as a few key obsta- cles that need to be overcome to advance this work to the next stage. 177 BIBLIOGRAPHY [1] E Anderson et al. “Block diagonalization using similarity renormalization group flow equations”. In: Physical Review C 77.3 (2008), p. 037001. [2] ER Anderson et al. “Operator evolution via the similarity renormalization group: The deuteron”. In: Physical Review C 82.5 (2010), p. 054001. [3] W. E. Arnoldi. “The principle of minimized iterations in the solution of the matrix eigenvalue problem”. In: Quarterly of Applied Mathematics 9.1 (1951), pp. 17–29. [4] C. Barbieri and A. Carbone. “Self-Consistent Green’s Function Approaches”. In: An Advanced Course in Computational Nuclear Physics. Ed. by M. Hjorth-Jensen, M. Lombardo, and U. van Kolck. Lecture Notes in Physics 936. Springer, 2017. Chap. 11. [5] Bruce R. Barrett, Petr Navr´atil, and James P. Vary. “Ab initio no core shell model”. In: Prog. Part. Nucl. Phys. 69.0 (2013), pp. 131 –181. [6] Sven Binder et al. “Ab initio path to heavy nuclei”. In: Physics Letters B 736 (2014), pp. 119–123. [7] Sven Binder et al. “Extension of coupled-cluster theory with a noniterative treatment of connected triply excited clusters to three-body Hamiltonians”. In: Physical Review C 88.5 (2013), p. 054319. [8] Sergio Blanes et al. “The Magnus expansion and some of its applications”. In: Physics reports 470.5-6 (2009), pp. 151–238. [9] SK Bogner, Richard J Furnstahl, and RJ Perry. “Similarity renormalization group for nucleon-nucleon interactions”. In: Physical Review C 75.6 (2007), p. 061001. [10] SK Bogner, RJ Furnstahl, and Achim Schwenk. “From low-momentum interactions to nuclear structure”. In: Progress in Particle and Nuclear Physics 65.1 (2010), pp. 94– 147. [11] SK Bogner, Thomas Tzu Szu Kuo, and A Schwenk. “Model-independent low mo- mentum nucleon interaction from phase shift equivalence”. In: Physics Reports 386.1 (2003), pp. 1–27. [12] SK Bogner et al. “Convergence of the Born series with low-momentum interactions”. In: Nuclear Physics A 773.3-4 (2006), pp. 203–220. [13] SK Bogner et al. “Nonperturbative shell-model interactions from the in-medium sim- ilarity renormalization group”. In: Physical Review Letters 113.14 (2014), p. 142501. 178 [14] Peter N Brown and Alan C Hindmarsh. “Reduced storage matrix methods in stiff ODE systems”. In: Applied Mathematics and Computation 31 (1989), pp. 40–91. [15] Gerhard Buchalla, Oscar Cat`a, and Claudius Krause. “On the power counting in effective field theories”. In: Physics Letters B 731 (2014), pp. 80–86. [16] Angelo Calci. “Evolved chiral Hamiltonians at the three-body level and beyond”. In: (2014). [17] James Chadwick. “The existence of a neutron”. In: Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 136.830 (1932), pp. 692–708. [18] Moody T Chu. “Scaled Toda-like flows”. In: Linear algebra and its applications 215 (1995), pp. 261–273. [19] S. A. Coon et al. “Convergence properties of ab initio calculations of light nuclei in a harmonic oscillator basis”. In: Phys. Rev. C 86 (5 Nov. 2012), p. 054002. [20] Carl Eckart and Gale Young. “The approximation of one matrix by another of lower rank”. In: Psychometrika 1.3 (1936), pp. 211–218. [21] A Ekstr¨om et al. “What is ab initio in nuclear theory?” In: Frontiers in Physics 11 (2023), p. 125. [22] DR Entem and R Machleidt. “Accurate charge-dependent nucleon-nucleon potential at fourth order of chiral perturbation theory”. In: Physical Review C 68.4 (2003), p. 041001. [23] Evgeny Epelbaum, Hermann Krebs, and Patrick Reinert. “High-precision nuclear forces from chiral EFT: state-of-the-art, challenges, and outlook”. In: Frontiers in Physics 8 (2020), p. 98. [24] EB Fel’dman. “On the convergence of the Magnus expansion for spin systems in periodic magnetic fields”. In: Physics Letters A 104.9 (1984), pp. 479–481. [25] Alexander L. Fetter and John D. Walecka. Quantum Theory of Many-Particle Sys- tems. Dover Publications, June 2003. isbn: 0486428273. [26] Paolo Finelli, Matteo Vorabbi, and Carlotta Giusti. “Microscopic Optical Potentials: recent achievements and future perspectives”. In: Journal of Physics: Conference Series. Vol. 2453. 1. IOP Publishing. 2023, p. 012026. [27] RJ Furnstahl. “The renormalization group in nuclear physics”. In: Nuclear Physics B-Proceedings Supplements 228 (2012), pp. 139–175. 179 [28] RJ Furnstahl et al. “Infrared extrapolations for atomic nuclei”. In: Journal of Physics G: Nuclear and Particle Physics 42.3 (2015), p. 034032. [29] Eskendr Gebrerufael, Angelo Calci, and Robert Roth. “Open-shell nuclei and excited states from multireference normal-ordered Hamiltonians”. In: Physical Review C 93.3 (2016), p. 031301. [30] G. Hagen et al. “Coupled-cluster computations of atomic nuclei”. In: Rept. Prog. Phys. 77.9 (2014), p. 096302. [31] G Hagen et al. “M. W loch, and P. Piecuch”. In: Phys. Rev. C 76.3 (2007). [32] Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp. “Finding structure with randomness: Stochastic algorithms for constructing approximate matrix decomposi- tions”. In: (2009). [33] K. Hebeler et al. “Improved nuclear matter calculations from chiral low-momentum interactions”. In: Phys. Rev. C 83 (3 Mar. 2011), p. 031301. [34] Kai Hebeler. “Momentum-space evolution of chiral three-nucleon forces”. In: Physical Review C 85.2 (2012), p. 021002. [35] Kai Hebeler. “Three-nucleon forces: Implementation and applications to atomic nuclei and dense matter”. In: Physics Reports 890 (2021), pp. 1–116. [36] Werner Heisenberg. “On the structure of atomic nuclei”. In: Z. Phys 77.1 (1932). [37] H Hergert. “In-medium similarity renormalization group for closed and open-shell nuclei”. In: Physica Scripta 92.2 (2016), p. 023002. [38] H Hergert and R Roth. “Treatment of the intrinsic Hamiltonian in particle-number nonconserving theories”. In: Physics Letters B 682.1 (2009), pp. 27–32. [39] H Hergert et al. “The in-medium similarity renormalization group: A novel ab initio method for nuclei”. In: Physics reports 621 (2016), pp. 165–222. [40] Heiko Hergert. “A guided tour of ab initio nuclear many-body theory”. In: Frontiers in Physics 8 (2020), p. 379. [41] Heiko Hergert et al. “In-medium similarity renormalization group approach to the nuclear many-body problem”. In: An Advanced Course in Computational Nuclear Physics. Springer, 2017, pp. 477–570. [42] Alan C Hindmarsh. “ODEPACK, a systemized collection of ODE solvers”. In: Scien- tific computing (1983). 180 [43] Alan C Hindmarsh et al. “SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers”. In: ACM Transactions on Mathematical Software (TOMS) 31.3 (2005), pp. 363–396. [44] Osamu Hino, Tomoko Kinoshita, and Rodney J. Bartlett. “Singular value decompo- sition applied to the compression of T3 amplitude for the coupled cluster method”. In: J. Chem. Phys. 121.3 (Jan. 2004), pp. 1206–1213. [45] Edward G. Hohenstein et al. “Rank-reduced coupled-cluster. III. Tensor hypercon- traction of the doubles amplitudes”. In: J. Chem. Phys. 156.5 (Feb. 2022), p. 054102. [46] Edward G. Hohenstein et al. “Rank reduced coupled cluster theory. II. Equation-of- motion coupled-cluster singles and doubles”. In: J. Chem. Phys. 151.16 (Apr. 2019), p. 164121. [47] Jason D Holt et al. “Counter terms for low momentum nucleon–nucleon interactions”. In: Nuclear Physics A 733.1-2 (2004), pp. 153–165. [48] J. Hoppe et al. “Importance truncation for the in-medium similarity renormalization group”. In: Phys. Rev. C 105 (3 Mar. 2022), p. 034324. [49] Calvin W Johnson. “Unmixing symmetries”. In: Physical Review Letters 124.17 (2020), p. 172502. [50] ED Jurgenson, Petr Navratil, and RJ Furnstahl. “Evolving nuclear many-body forces with the similarity renormalization group”. In: Physical Review C 83.3 (2011), p. 034301. [51] Eric Donald Jurgenson. “Applications of the similarity renormalization group to the nuclear interaction”. In: arXiv preprint arXiv:0912.2937 (2009). [52] G.P. Kamuntaviˇcius et al. “The general harmonic-oscillator brackets: compact ex- pression, symmetries, sums and Fortran code”. In: Nuclear Physics A 695.1 (2001), pp. 191–201. [53] Stefan Kehrein. The flow equation approach to many-particle systems. Vol. 217. Springer, 2007. [54] Tomoko Kinoshita, Osamu Hino, and Rodney J. Bartlett. “Singular value decom- position approach for the approximate coupled-cluster method”. In: J. Chem. Phys. 119.15 (Oct. 2003), pp. 7756–7762. [55] S. K¨onig et al. “Ultraviolet extrapolations in finite oscillator bases”. In: Phys. Rev. C 90 (6 Dec. 2014), p. 064007. 181 [56] Werner Kutzelnigg. “Quantum chemistry in Fock space. I. The universal wave and energy operators”. In: The Journal of Chemical Physics 77.6 (1982), pp. 3081–3097. [57] Werner Kutzelnigg. “Quantum chemistry in Fock space. III. Particle-hole formalism”. In: The Journal of chemical physics 80.2 (1984), pp. 822–830. [58] Werner Kutzelnigg and Sigurd Koch. “Quantum chemistry in Fock space. II. Effec- tive Hamiltonians in Fock space”. In: The Journal of chemical physics 79.9 (1983), pp. 4315–4335. [59] C. Lanczos. “An iteration method for the solution of the eigenvalue problem of linear differential and integral operators”. In: J. Res. Nat’l Bur. Std. 45 (1950), p. 255. [60] G Peter Lepage. “What is renormalization?” In: arXiv preprint hep-ph/0506330 (2005). [61] GP Lepage. “How to renormalize the Schrodinger equation”. In: Nuclear Physics (ed. by CA Bertulani et al.) (1997), p. 135. [62] Micha(cid:32)l Lesiuk. “Near-Exact CCSDT Energetics from Rank-Reduced Formalism Sup- plemented by Non-iterative Corrections”. In: J. Chem. Theory Comput. 17.12 (Dec. 2021), pp. 7632–7647. [63] Micha(cid:32)l Lesiuk. “Quintic-scaling rank-reduced coupled cluster theory with single and double excitations”. In: J. Chem. Phys. 156.6 (Oct. 2022), p. 064103. [64] Justin G Lietz et al. “Computational nuclear physics and post hartree-fock methods”. In: An Advanced Course in Computational Nuclear Physics: Bridging the Scales from Quarks to Neutron Stars (2017), pp. 293–399. [65] R Machleidt. “High-precision, charge-dependent Bonn nucleon-nucleon potential”. In: Physical Review C 63.2 (2001), p. 024001. [66] R Machleidt. “The theory of nuclear forces: Is the never-ending story coming to an end?” In: arXiv preprint nucl-th/0701077 (2007). [67] Re Machleidt. “The meson theory of nuclear forces and nuclear structure”. In: Ad- vances in nuclear physics. Springer, 1989, pp. 189–376. [68] Ruprecht Machleidt and David Rodrignez Entem. “Chiral effective field theory and nuclear forces”. In: Physics Reports 503.1 (2011), pp. 1–75. [69] Wilhelm Magnus. “On the exponential solution of differential equations for a linear operator”. In: Communications on pure and applied mathematics 7.4 (1954), pp. 649– 673. 182 [70] S. N. More et al. “Universal properties of infrared oscillator basis extrapolations”. In: Phys. Rev. C 87 (4 Apr. 2013), p. 044326. [71] TD Morris, NM Parzuchowski, and SK Bogner. “Magnus expansion and in-medium similarity renormalization group”. In: Physical Review C 92.3 (2015), p. 034331. [72] P Navr´atil, Gintautas Pranciˇskus Kamuntaviˇcius, and Bruce R Barrett. “Few-nucleon systems in a translationally invariant harmonic oscillator basis”. In: Physical Review C 61.4 (2000), p. 044001. [73] Petr Navratil. “Local three-nucleon interaction from chiral effective field theory”. In: Few-Body Systems 41.3-4 (2007), pp. 117–140. [74] Petr Navr´atil et al. “Unified ab initio approaches to nuclear structure and reactions”. In: Phys. Scripta 91.5 (2016), p. 053002. [75] A Nogga et al. “Spectra and binding energy predictions of chiral interactions for Li 7”. In: Physical Review C 73.6 (2006), p. 064002. [76] Andreas Nogga, Scott K. Bogner, and Achim Schwenk. “Low-momentum interaction in few-nucleon systems”. In: Phys. Rev. C 70 (6 Dec. 2004), p. 061002. [77] D. Odell, T. Papenbrock, and L. Platter. “Infrared extrapolations of quadrupole mo- ments and transitions”. In: Phys. Rev. C 93 (4 Apr. 2016), p. 044331. [78] Robert M. Parrish et al. “Rank reduced coupled cluster theory. I. Ground state en- ergies and wavefunctions”. In: J. Chem. Phys. 150.16 (Apr. 2019), p. 164118. [79] N. M. Parzuchowski et al. “Ab initio electromagnetic observables with the in-medium similarity renormalization group”. In: Phys. Rev. C 96 (3 Sept. 2017), p. 034324. [80] Maria Piarulli and Ingo Tews. “Local nucleon-nucleon and three-nucleon interactions within chiral effective field theory”. In: Frontiers in Physics 7 (2020), p. 245. [81] William H Press et al. Numerical recipes 3rd edition: The art of scientific computing. Cambridge university press, 2007. [82] P. Ring and P. Schuck. The nuclear many-body problem. New York: Springer-Verlag, 1980. [83] David Rodriguez Entem, Ruprecht Machleidt, and Yevgen Nosyk. “Nucleon-nucleon scattering up to N5LO in chiral effective field theory”. In: Frontiers in Physics 8 (2020), p. 57. 183 [84] Robert Roth and Petr Navr´atil. “Ab Initio Study of Ca 40 with an Importance- Truncated No-Core Shell Model”. In: Physical review letters 99.9 (2007), p. 092501. [85] Robert Roth et al. “Evolved chiral N N+ 3 N Hamiltonians for ab initio nuclear structure calculations”. In: Physical Review C 90.2 (2014), p. 024325. [86] Robert Roth et al. “Medium-mass nuclei with normal-ordered chiral N N+ 3 N in- teractions”. In: Physical Review Letters 109.5 (2012), p. 052501. [87] Robert Roth et al. “Similarity-Transformed Chiral N N+ 3 N Interactions for the Ab Initio Description of C 12 and O 16”. In: Physical Review Letters 107.7 (2011), p. 072501. [88] Micah D Schuster et al. “Operator evolution for ab initio theory of light nuclei”. In: Physical Review C 90.1 (2014), p. 011301. [89] Isaiah Shavitt and Rodney J Bartlett. Many-body methods in chemistry and physics: MBPT and coupled-cluster theory. Cambridge university press, 2009. [90] Vittorio Som`a. “Self-Consistent Green’s Function Theory for Atomic Nuclei”. In: Frontiers in Physics 8 (2020), p. 340. [91] H. P. Stapp, T. J. Ypsilantis, and N. Metropolis. “Phase-Shift Analysis of 310-Mev Proton-Proton Scattering Experiments”. In: Phys. Rev. 105 (1 Jan. 1957), pp. 302– 310. [92] SR Stroberg et al. “Ab initio limits of atomic nuclei”. In: Physical Review Letters 126.2 (2021), p. 022501. [93] SR Stroberg et al. “Ground and excited states of doubly open-shell nuclei from ab initio valence-space Hamiltonians”. In: Physical Review C 93.5 (2016), p. 051301. [94] SR Stroberg et al. “Nucleus-dependent valence-space approach to nuclear structure”. In: Physical Review Letters 118.3 (2017), p. 032502. [95] Alexander Tichai et al. “Tensor-decomposition techniques for ab initio nuclear struc- ture calculations: From chiral nuclear potentials to ground-state energies”. In: Physical Review C 99.3 (2019), p. 034320. [96] AJ Tropiano, SK Bogner, and RJ Furnstahl. “Operator evolution from the similar- ity renormalization group and the Magnus expansion”. In: Physical Review C 102.3 (2020), p. 034005. [97] CM Vincent and SC Phatak. “Accurate momentum-space method for scattering by nuclear and Coulomb potentials”. In: Physical Review C 10.1 (1974), p. 391. 184 [98] Sergey Voronin and Per-Gunnar Martinsson. “RSVDPACK: An implementation of randomized algorithms for computing the singular value, interpolative, and CUR decompositions of matrices on multi-core and GPU architectures”. In: arXiv preprint arXiv:1502.05366 (2015). [99] Meng Wang et al. “The AME 2020 atomic mass evaluation (II). Tables, graphs and references”. In: Chinese Physics C 45.3 (2021), p. 030003. [100] F. Wegner. “Flow equations for Hamiltonians”. In: Ann. Phys. (Leipzig) 3 (1994), p. 77. [101] Franz J Wegner. “Flow equations for Hamiltonians”. In: Physics Reports 348.1-2 (2001), pp. 77–89. [102] Steven Weinberg. “Effective chiral Lagrangians for nucleon-pion interactions and nu- clear forces”. In: Nuclear Physics B 363.1 (1991), pp. 3–18. [103] Steven Weinberg. “Phenomenological lagrangians”. In: Physica a 96.1-2 (1979), pp. 327– 340. [104] Steven Weinberg. “What is quantum field theory, and what did we think it is?” In: arXiv preprint hep-th/9702027 (1997). [105] K. A. Wendt. “Similarity renormalization group evolution of three-nucleon forces in a hyperspherical momentum representation”. In: Phys. Rev. C 87 (6 June 2013), p. 061001. [106] K. A. Wendt et al. “Infrared length scale and extrapolations for the no-core shell model”. In: Phys. Rev. C 91 (6 June 2015), p. 061301. [107] Steven R. White. “Numerical canonical transformation approach to quantum many- body problems”. In: J. Chem. Phys. 117.16 (2002), pp. 7472–7482. [108] G. C. Wick. “The Evaluation of the Collision Matrix”. In: Phys. Rev. 80 (2 Oct. 1950), pp. 268–272. [109] Kenneth G Wilson. “Problems in physics with many scales of length”. In: Scientific American 241.2 (1979), pp. 158–179. [110] Robert B Wiringa, VGJ Stoks, and R Schiavilla. “Accurate nucleon-nucleon potential with charge-independence breaking”. In: Physical Review C 51.1 (1995), p. 38. [111] J. M. Yao et al. “Ab Initio Treatment of Collective Correlations and the Neutrinoless Double Beta Decay of 48Ca”. In: Phys. Rev. Lett. 124 (23 June 2020), p. 232501. 185 [112] Hideki Yukawa and Shoichi Sakata. “On the Theory of the β-Disintegration and the Allied Phenomenon”. In: Proceedings of the Physico-Mathematical Society of Japan. 3rd Series 17 (1935), pp. 467–479. [113] B Zhu, R Wirth, and H Hergert. “Singular value decomposition and similarity renor- malization group evolution of nuclear interactions”. In: Physical Review C 104.4 (2021), p. 044002. 186