THEORETICAL AND COMPUTATIONAL IMPROVEMENTS TO THE IN-MEDIUM SIMILARITY RENORMALIZATION GROUP By Jacob Davison A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Physics—Doctor of Philosophy Computational Mathematics, Science and Engineering—Dual Major 2023 ABSTRACT There has been much progress in ab initio nuclear many-body theory in recent years, including a first attempt on an ab initio based mass table and a calculation of the neutron skin of 208 Pb based on two- plus three-nucleon interactions from chiral Effective Field Theory, with a sophisticated statistical uncertainty quantification. Going forward, there is a need to continue extending the capabilities of ab initio nuclear many-body theory to observables in medium-mass and heavy open-shell nuclei, enhance the description of exotic nuclei, and to build on theoretical and computational advances in nuclear structure to tackle nuclear dynamics. These advances in ab initio nuclear many-body theory, along with new methods for carrying controlled uncertainties into the relevant calculations, will provide new insight into nuclear forces and phenomena, as well as important input for searches for Beyond-Standard Model physics or nuclear astrophysics. All of these goals requires theoretical and computational improvements of current nuclear many-body methods in order to enhance the stability of predictions, determine proper theoretical error bars, and to provide greater computational efficiency. In this body of work, we pursue these directions for the In-Medium Similarity Renormalization Group (IMSRG), which has become an important tool in ab initio nuclear many-body theory over the past decade. We will discuss three major contributions in these directions. First, we present a model for the future infrastructure of IMSRG production codes using tensor network architectures, and show that it can soften the memory requirements and favorably improve the efficiency of IMSRG calculations. Second, we present a methodological improvement, the so-called reference state ensemble, that allows us to mitigate truncation errors in the IMSRG flow without performing calculations at a higher truncation rank. We show that the reference state ensemble improves the stability and accuracy of the IMSRG flow by “informing” the underlying operator basis about the features of the many-body system’s excitations. Last but not least, we use a data-driven technique called Dynamic Mode Decomposition (DMD) for emulating the IMSRG solution which reduces the cost of complete integration to convergence from hundreds or thousands of iterations to only tens of iterations. We also present a parametric emulation technique, powered by DMD and trained on IMSRG data, which we show can robustly predict the IMSRG results for Hamiltonians, including chiral two- plus three-nucleon Hamiltonians with more than 20 parameters. We show that the parametric emulator can produce a global sensitivity analysis of statistical significance in only a few minutes, where full IMSRG(2) calculations would take thousands of years—improving the feasibility of uncertainty quantification for the IMSRG. Copyright by JACOB DAVISON 2023 To my wife, Sophie. To my family. To my friends. To my advisors, my mentors, and my teachers throughout my career as a student. “Serious students are not commonly plunged into fits of despair on turning the pages of a textbook and discovering that some further topic is known to the author but not to the student. Usually the students struggle a little, acquire the new knowledge, and, following an ancient human tradition, continue to turn the pages.” – Carl Sagan, Pale Blue Dot: A Vision of the Human Future in Space v ACKNOWLEDGEMENTS Foremost, I would like to acknowledge my advisor, Heiko Hergert. Without his support, encour- agement, patience, and gracious flexibility, I would have been quite lost on my foray into graduate physics. Heiko is a veritable font of knowledge of all things many-body physics, and more. Thank you, Heiko, for guiding me on this journey. I would also like to thank Scott Bogner, Alexei Bazavov, Longxiu Huang, and Jaideep Singh, for their commitment to serve on my guidance committee. I want to thank the collaborating members of what I affectionately refer to as the “Heiko group,” to whom I presented my research progress every week and received thought-provoking advice and fruitful discussions. These members are Jiangming Yao, Boyao Zhu, and additionally Roland Wirth. When I first arrived at MSU in the summer of 2018, I joined a rather large, exciting collaboration of three professors, Heiko Hergert (who, of course, would go on to be my direct advisor), Scott Bogner, and Morten Hjorth-Jensen, who were running a kind of “summer school” for incoming graduate students interested in their research. Their excitement for this field was infectious, and this introduction was absolutely critical and foundational for my success in the program, and eventually my own research. I want to thank Heiko, Scott, and Morten for being great advisors together. I also want to thank Justin Lietz for his mentorship in this initial summer introduction, and for being my WaMPS peer mentor when I first arrived as a new graduate student. Finally, I thank fellow graduate students Jane Kim, Julie Butler, and Yani Udiani for fruitful discussions and help in classes. The main reason I pursued this doctoral program was my experience as an undergraduate in physics at Central Michigan University. I want to thank the entire CMU Department of Physics for providing an enriching undergraduate experience. I especially want to thank Georgios Perdikakis for his mentorship, as well as Alan Jackson, Christopher Tycner, Aaron LaCluyzé and Marco Fornari for helping me decide how to proceed in my research career. I want to thank the entire MSU Department of Physics, although much larger than CMU, for providing an enriching graduate experience. Thanks to Kim Crosslan for navigating the vi bureaucracies of the graduate school. I also want to thank the Theory Alliance here at the Facility for Rare Isotope Beams for their resources and support for aspiring nuclear theorists. There is a truly amazing collection of minds here at MSU; what a privilege to have learned from them. My family is important to me. Without them, I certainly would not be where I am. This list is not exhaustive, but thanks to my mom Mona and my sister Jenna, my dad Jamie, Grandpa Ray and Grandma Linda, Uncle Stein and Aunt Emilie, Uncle Dave and Aunt Lina, Uncle Joe, my cousins Madelyn, Linnea, Gibson, and Sam. Thanks to my father-in-law Joe and my mother-in-law Tracy. Thanks to my dog Bean. Your unconditional encouragement has been a sail in the wind and a rudder in the sea. Love you all! Thanks to my best friends, Austin, Austin, Alan, Ben, Brady, Brian, Jake, Peter, and Zach, all of whom I have known, in some capacity, since kindergarten. We did it boys. Thanks to my wife Sophie, who I met at Central Michigan University. She was a chemistry major and I was a physics major. A love story for the ages. Since we met, we have completed every progression of our professional careers together. We pursued and accomplished an undergraduate research internship together, we finished our Bachelor’s degrees together, we attended MSU together for our graduate degrees, and now I like to think we are finishing this dissertation together (although, I assure you that I wrote this dissertation myself). Love you, Sophie. vii TABLE OF CONTENTS LIST OF ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x CHAPTER 1 OVERVIEW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 CHAPTER 2 MANY-BODY THEORY IN SECOND QUANTIZATION . . . . . . . . 8 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 Many-Body Basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Annihilation and Creation Operators . . . . . . . . . . . . . . . . . . . . . . . 10 2.4 Normal-Ordering and Wick Contractions . . . . . . . . . . . . . . . . . . . . . 12 2.5 Wick’s Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.6 Correlated Reference States and Irreducible Density Matrices . . . . . . . . . . 16 2.7 Normal-ordering the Two-body Hamiltonian . . . . . . . . . . . . . . . . . . . 17 2.8 Density Matrix Normal-ordering . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.9 Slater Determinant Occupation Number Representation . . . . . . . . . . . . . 20 2.10 Full Configuration Interaction Theory to Solve Model Hamiltonians . . . . . . 21 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 CHAPTER 3 IN-MEDIUM SIMILARITY RENORMALIZATION GROUP . . . . . . 24 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Similarity Renormalization Group . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3 IMSRG Formalism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.4 Generators for the IMSRG . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 CHAPTER 4 MODEL HAMILTONIANS . . . . . . . . . . . . . . . . . . . . . . . . 38 4.1 The Pairing Model Plus Particle-Hole Model . . . . . . . . . . . . . . . . . . . 38 4.2 Chiral Effective Field Theory for Nuclear Interactions . . . . . . . . . . . . . . 42 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 CHAPTER 5 BUILDING AN IMSRG MANY-BODY CODE USING TENSOR NETWORKS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.2 Tensor Network Diagrams . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.3 Factorization with Tensor Train Decomposition . . . . . . . . . . . . . . . . . 52 5.4 Tensor Networks in the IMSRG . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.5 Tensor-factors IMSRG: Python Implementation for the P3H Model . . . . . . 61 5.6 Tensorized C++ IMSRG: C++ Implementation for the P3H Model . . . . . . . . 64 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 CHAPTER 6 OPTIMAL REFERENCE STATE ENSEMBLES FOR THE IN-MEDIUM SIMILARITY RENORMALIZATION GROUP . . . . . . . . . . . . . . 69 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 viii 6.2 Truncation Error in the IMSRG . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.3 Formulation of the Reference State Ensemble . . . . . . . . . . . . . . . . . . 70 6.4 Ensemble Optimization Scheme . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.5 Regularization of Generator Divergences . . . . . . . . . . . . . . . . . . . . . 73 6.6 Improving the Eigenspectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.7 Ensemble Consistency with Non-Targeted Observables . . . . . . . . . . . . . 77 6.8 Using Exact Eigenstates as References . . . . . . . . . . . . . . . . . . . . . . 79 6.9 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 CHAPTER 7 KOOPMAN OPERATOR FOR THE IN-MEDIUM SIMILARITY RENORMALIZATION GROUP . . . . . . . . . . . . . . . . . . . . . . 86 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 7.2 The Koopman Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 7.3 Dynamic Mode Decomposition for IMSRG . . . . . . . . . . . . . . . . . . . 87 7.4 Emulating the IMSRG Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 7.5 Emulating the Parametric IMSRG Flow . . . . . . . . . . . . . . . . . . . . . 97 7.6 Interpolating DMD using RBFs . . . . . . . . . . . . . . . . . . . . . . . . . . 99 7.7 Optimized Training for Large Bases using Streaming SVD . . . . . . . . . . . 106 7.8 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 CHAPTER 8 GLOBAL SENSITIVITY ANALYSIS USING EMULATION OF THE IMSRG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 8.2 GSA with a Chiral N3LO Nucleon-Nucleon Interaction . . . . . . . . . . . . . 115 8.3 GSA with Chiral NNLO Two- Plus Three-Nucleon Interactions . . . . . . . . . 116 8.4 Discussion and Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 CHAPTER 9 CONCLUSIONS AND OUTLOOK . . . . . . . . . . . . . . . . . . . . 127 APPENDIX A TENSOR TRAIN DECOMPOSITION OF THE OCCUPATION FAC- TORS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 APPENDIX B MULTI-REFERENCE IMSRG(2) FLOW EQUATIONS . . . . . . . . . 132 APPENDIX C CONTRIBUTED CODE PROJECTS . . . . . . . . . . . . . . . . . . . 135 ix LIST OF ABBREVIATIONS IMSRG In-Medium Similarity Renormalization Group DMD Dynamic Mode Decomposition IE Interpolation Engine IEM Interpolation Engine Model rIEM Reduced Interpolation Engine Model UQ Uncertainty Quantification FCI Full Configuration Interaction IBCs Irreducible Brillouin Conditions SD Slater determinant LECs Low-energy constants QCD Quantum chromodynamics LO Leading order NLO Next-to leading order NNLO Next-to-next-to leading order N3 LO Next-to-next-to-next-to leading order TT Tensor Train IDMs Irreducible Density Matrices SVD Singular Value Decomposition rKOI Reduced Koopman Operator Interpolation RBF Radial Basis Function SKL Sequential Karhunen-Loéve GSA Global Sensitivity Analysis PCA Principal Component Analysis x CHAPTER 1 OVERVIEW There has been much progress in ab initio [3] nuclear many-body theory in the past decade, which is illustrated in Figure 1.1. For example, there has been a first attempt on an ab initio based mass table [11] and the neutron skin of 208 Pb has been computed based on two- plus three-nucleon interactions from chiral Effective Field Theory [8, 9], with a sophisticated statistical uncertainty quantification [6]. The next tasks of ab initio nuclear theory are (i) to extend its capabilities into the domain of medium-mass and heavy open-shell nuclei, which exhibit strong collective correlations associated with intrinsic deformation, (ii) to enhance the description of exotic nuclei, where weak-binding and the coupling to continuum degrees of freedom lead to subtle new phenomena, and (iii) to build on the theoretical and computational advances in nuclear structure in order to tackle nuclear dynamics. In this way, ab initio nuclear many-body calculations with controlled uncertainties will provide crucial new input to tackle grand scientific challenges, like the search for Beyond-Standard Model physics [2], or the improvement of our understanding of supernovae and neutron star mergers, and their detailed role in nucleosynthesis. All of this requires theoretical and computational improvements of current nuclear many-body methods, in order to enhance the stability of predictions, determine proper theoretical error bars, and to provide much greater computational efficiency. In the present work, we pursue next steps in these directions for the In-Medium Similarity Renormalization Group (IMSRG), which has become an important tool in ab initio nuclear many-body theory over the past decade. Chapter 2 introduces some essential aspects of nuclear many-body theory theory in second quantization that are needed for the discussion in this work; following that is Chapter 3, which is dedicated to a detailed presentation of the IMSRG. The IMSRG is an ab initio, or first principles, method for solving the nuclear many-body problem. Importantly, the IMSRG can also be used to pre-process interactions and operators for use in other methods (Valence-Space IMSRG [12], 1 IM-GCM for Double Beta Decay [13], IN-NCSM [4]). Functionally, the IMSRG decouples a target reference state, e.g. the ground state, from all other excitations in the many-body basis via a continuous, unitary transformation that effectively block-diagonalizes the Hamiltonian (Chapter 4 discusses the model Hamiltonians used for the analyses in this work). Excitations are gradually suppressed according to a particular decoupling scheme which separates diagonal elements from off-diagonal elements in a chosen operator basis. Typically, the most convenient operator basis is generated through a process called normal-ordering in the language of second quantization (Chapter 2). Then, instead of calculating the IMSRG transformation explicitly, we solve for the transformation implicitly through integration of the IMSRG flow equations (Chapter 3.3), which involves a commutator of the Hamiltonian and an object called the generator; this object generates the desired behavior of off-diagonal elements of the Hamiltonian over the integration range. The commutator itself generates terms of increasingly higher-body interactions, and tracking all of them is not computationally feasible. Thus, we must introduce a truncation scheme. In the IMSRG(2), for example, we truncate all operators in the chosen operator basis at the two-body level. The IMSRG is favorable due its polynomial scaling in basis size. Additionally, the IMSRG is built in a framework which allows systematic improvement toward exact solutions via the increase in truncation rank, which introduces increasingly higher-fidelity physics information. Chapter 5 introduces a model for the future infrastructure of the IMSRG production codes, the motivation for which is built on the exchange of tensor libraries on the backend. We also show how the Tensor Train decomposition (Chapter 5.3) can be used to rewrite particular objects in the flow equations, which reduces the memory demand of those objects by several orders of magnitude in extreme cases. This technique holds the key for efficiently leveraging the ever-evolving hardware landscape, and lets us benefit from the work of library vendors who know best how to adapt tensor calculations to new platforms. Chapter 6 discusses the use reference state ensembles for the IMSRG, which is a technique for potentially mitigating truncation errors while avoiding the computational cost of performing IMSRG at a higher truncation rank. We show that error due to truncation at two-body operators 2 �� �� ���� �=�� �� �=��� �� �� �=�� �� � �� �=�� �� �=�� �=�� �� �=�� ���� ���� �� �=�� �=�� ���� ���� � ���� �=�� ���� �=�� ���� � �� �� �� �� �� �� �� �� �� �� �� ��� ��� ��� ��� ��� ��� � Figure 1.1 Figured updated from [5]. Each filled square represents a nucleus that has been successfully calculated using IMSRG. manifests as a loss of unitarity in the IMSRG(2) transformation. The severity of this loss depends on the choice of operator basis for the particular problem we study, and is related to the importance of induced terms which have been truncated. The goal of the reference state ensemble is to introduce additional physics information to the IMSRG(2) via correlations in the operator basis. In the reference state ensemble scheme, we construct a reference state from statistically-mixed many- body basis configurations (e.g. approximate ground- and excited-state configurations), and then use that reference as a starting point for the IMSRG(2). In this way, we introduce more descriptive physical information to the standard operator basis of the IMSRG. We show that the reference state ensemble is a controllable method for improving divergences in the IMSRG(2) related to truncation error, as well as improving accuracy in the evolved energy spectrum. 3 Chapter 7 presents a method for emulating the IMSRG solution using an approach called Dy- namic Mode Decomposition (DMD), which is based on Koopman operator theory. The Koopman operator is a linear operator which advances measurement functions of a nonlinear dynamical system forward in time [7, 1]. The eigendecomposition of the Koopman operator characterizes the dynamics of the system, as expressed through measurements of observables that depend on the system’s state variables. Thus, it can be understood as shifting the description of the system to coordinates where the nonlinear dynamics look linear, analogous to canonical transformations in the phase-space formulation of classical mechanics. DMD is a finite-measurement, discrete-time method for approximating the Koopman operator. The DMD operator is formulated as the best-fit linear operator which propagates the dynamical system forward one step in time [10]. The DMD is informed by direct measurements of the evolving system. We show that the DMD operator, built from IMSRG(2) observations and informed by simple physical constraints on the predicted dynamics, offers a method by which we may construct a true emulator for the IMSRG(2) solution. DMD emulation allows us to forecast the IMSRG(2) evolution for a particular Hamiltonian with high accuracy and short computational wall-time, at the cost of a few polynomial-time iterations of the IMSRG(2) solution which serve as “training points” for the emulator. We also present an interpolation engine (IE) for the IMSRG, constructed from DMD, which we use to train an interpolative model of the parametric IMSRG, i.e. a model where the inputs include Hamilto- nian coupling parameters as well as the dynamical variable. We show that fitted IE model (IEM) produces predictions with robust accuracy for up to 24 coupling parameters in the Hamiltonian. The parametric emulation of the IMSRG(2) offer a strategy for undertaking the large volumes of calculations required for uncertainty quantification (UQ) of the IMSRG, so that we may accurately assign theoretical constraints on the output. Chapter 9 demonstrates this capability by performing sensitivity analyses of observables computed in IMSRG(2) to the parameters of the underlying chiral interactions. We show that a statistically significant number of calculations can be performed in just a few minutes using the parametric emulator, whereas the same number of full-scale IMSRG(2) calculations would require thousands of years. The feasibility of these sensitivity analyses, using 4 parametric emulation, is a major step toward feasibility of UQ for the IMSRG. 5 BIBLIOGRAPHY [1] Steven L. Brunton et al. Modern Koopman Theory for Dynamical Systems. 2021. arXiv: 2102.12086 [math.DS]. [2] Vincenzo Cirigliano et al. Neutrinoless Double-Beta Decay: A Roadmap for Matching Theory to Experiment. 2022. arXiv: 2203.12169 [hep-ph]. [3] A. Ekström et al. “What is ab initio in nuclear theory?” In: Frontiers in Physics 11 (2023). issn: 2296-424X. doi: 10.3389/fphy.2023.1129094. url: https://www.frontiersin.org/ articles/10.3389/fphy.2023.1129094. [4] Eskendr Gebrerufael et al. “Ab Initio Description of Open-Shell Nuclei: Merging No- Core Shell Model and In-Medium Similarity Renormalization Group”. In: Phys. Rev. Lett. 118 (15 Apr. 2017), p. 152503. doi: 10.1103/PhysRevLett.118.152503. url: https: //link.aps.org/doi/10.1103/PhysRevLett.118.152503. [5] Heiko Hergert. “A Guided Tour of ab initio Nuclear Many-Body Theory”. In: Frontiers in Physics 8 (2020), p. 379. [6] Baishan Hu et al. “Ab initio predictions link the neutron skin of 208Pb to nuclear forces”. In: Nature Physics 18.10 (Aug. 2022), pp. 1196–1200. doi: 10.1038/s41567-022-01715-8. url: https://doi.org/10.1038%5C%2Fs41567-022-01715-8. [7] B. O. Koopman. “Hamiltonian Systems and Transformation in Hilbert Space”. In: Proceedings of the National Academy of Sciences 17.5 (1931), pp. 315–318. doi: 10.1073/pnas.17.5.315. eprint: https://www.pnas.org/doi/pdf/10.1073/pnas.17.5.315. url: https://www.pnas.org/doi/abs/10.1073/pnas.17.5.315. [8] R. Machleidt and D.R. Entem. “Chiral effective field theory and nuclear forces”. In: Physics Reports 503.1 (June 2011), pp. 1–75. doi: 10.1016/j.physrep.2011.02.001. url: https://doi.org/10.1016%5C%2Fj.physrep.2011.02.001. [9] Maria Piarulli and Ingo Tews. “Local Nucleon-Nucleon and Three-Nucleon Interactions Within Chiral Effective Field Theory”. In: Frontiers in Physics 7 (2020). issn: 2296- 424X. doi: 10.3389/fphy.2019.00245. url: https://www.frontiersin.org/articles/10.3389/ fphy.2019.00245. [10] Peter J. Schmid. “Dynamic mode decomposition of numerical and experimental data”. In: Journal of Fluid Mechanics 656 (2010), pp. 5–28. doi: 10.1017/S0022112010001217. [11] S. R. Stroberg et al. “Ab Initio Limits of Atomic Nuclei”. In: Phys. Rev. Lett. 126 (2 Jan. 2021), p. 022501. doi: 10.1103/PhysRevLett.126.022501. url: https://link.aps.org/doi/10. 1103/PhysRevLett.126.022501. 6 [12] S. Ragnar Stroberg et al. “Nonempirical Interactions for the Nuclear Shell Model: An Update”. In: Annual Review of Nuclear and Particle Science 69.1 (2019), pp. 307– 362. doi: 10.1146/annurev-nucl-101917-021120. eprint: https://doi.org/10.1146/ annurev-nucl-101917-021120. url: https://doi.org/10.1146/annurev-nucl-101917-021120. [13] J. M. Yao et al. “Ab Initio Treatment of Collective Correlations and the Neutrinoless Double Beta Decay of 48 Ca”. In: Phys. Rev. Lett. 124 (23 June 2020), p. 232501. doi: 10.1103/ PhysRevLett.124.232501. url: https://link.aps.org/doi/10.1103/PhysRevLett.124.232501. 7 CHAPTER 2 MANY-BODY THEORY IN SECOND QUANTIZATION 2.1 Introduction Much of the nuclear theory toolbox leveraged by this body of work is constructed within the “second quantization” formalism for many-body quantum theory. Second quantization is concerned with modeling a space which contains multiple quantum systems, e.g. the nucleus of an atom for which the number of particles is not constant. Second quantization, which exists inside of Fock space, provides tools which allow us to navigate between individual Hilbert spaces, and is the most convenient framework for quantum many-body theory. The purpose of this chapter is to provide context and definitions for these tools which are used in following chapters. Much of the information here is compiled from Shavitt and Bartlett [6] and may take on the flavor of quantum chemistry as a result; however, the notation and techniques are general enough to be understood in the context of nuclear theory. 2.2 Many-Body Basis The “many-body basis” may be understood as a chosen collection of states, varying particle number, which exist inside the Fock space. Formally, the fermionic Fock space is the direct sum over all antisymmetrized tensor products of Hilbert spaces of quantum states up to 𝑁 particle numbers,   F (𝐻) = C ⊕ 𝐻 ⊕ A (𝐻 ⊗ 𝐻) ⊕ . . . A 𝐻 ⊗𝑁 , (2.1) where A is an operator which generates the antisymmetrization. Equivalent notation which compresses the definition reads Ê𝑁 F𝐴 (𝐻) = A𝐻 ⊗𝑛 . (2.2) 𝑛=0 The Fock states are acted upon by operators which also exist inside the Fock space. For the purposes of this work, we are interested in fermionic Fock space because the techniques we apply are meant to model nuclear systems. Fermionic Hilbert spaces are antisymmetric under the exchange of particles. The complete single particle basis spans the Hilbert space of a single 8 particle, and is uniquely characterized by a set of quantum numbers that are associated with the static and dynamical properties of that particle. [6]. The Hilbert space for an 𝐴-body system is constructed as a product of single-particle Hilbert spaces, which is spanned by tensor products of the single-particle states. We notate the single particle basis 𝜙 𝑝 , which encodes the set of all possible single particle states indexed by 𝑝 = 1, 2, 3, . . . . We assume that the single particle basis is orthonormal. Schematically, we might write the Hilbert space for an 𝐴-body system as, Ì 𝐴 𝐻 ( 𝐴) = 𝜙𝑛 . (2.3) 𝑛=0 The antisymmetrizer A enforces that condition that product states, elements of the 𝐴-body Hilbert space in Equation 2.3, be antisymmetric under the exchange of particle, e.g. for a two- particle product state, 𝜙 𝑞(1) ⊗ 𝜙 (2) (1) 𝑝 = −𝜙 𝑝 ⊗ 𝜙 𝑞 . (2) (2.4) The antisymmetrizer A generates the so-called Slater determinant, such that   √︂ 1   (1) (2) (1) (2) Φ = A 𝜙 (1) ⊗ 𝜙 (2) = 𝜙 ⊗ 𝜙 − 𝜙 ⊗ 𝜙 𝑝 𝑞 2 𝑝 𝑞 𝑞 𝑝 √︂ (1) (2) (2.5) 1 𝜙𝑝 𝜙𝑝 = . 2 𝜙 (1) 𝜙 (2) 𝑞 𝑞 For brevity, we will often write a Slater determinant Φ as a ket of single particle states, or maybe single particle indices, Φ = |𝜙 𝑝 𝜙 𝑞 𝜙𝑟 . . .⟩ , (2.6) where, intuitively, the single particle states collected in this string are occupied within this particular state. The antisymmetrizer has been implicitly absorbed into the notation. Note that e.g. |𝜙 𝑝 𝜙 𝑞 ⟩ = − |𝜙 𝑞 𝜙 𝑝 ⟩. In principle, there an infinite number of unoccupied states which we could write for this particular state, but there is no reason to write all of them. Additionally, we might choose to drop the 𝜙 for brevity, settling on raw indices 𝑝, 𝑞, 𝑟, . . . to represent these single particle states which are occupied in the Slater determinant. 9 Of course, the order in which the single particle states appear in this string is deliberate, and the fermionic nature of the Slater determinant guarantees that permuting the order of these states results in a change of phase. Written succinctly [6], ˆ 𝑃ˆ | 𝑝𝑞𝑟 · · ·⟩ = (−1) 𝜎( 𝑃) | 𝑝𝑞𝑟 · · ·⟩ . (2.7) Permuting the indices of Slater determinant | 𝑝𝑞𝑟 · · ·⟩ via 𝑃ˆ results in a phase change of order 𝜎( 𝑃), ˆ which we must track whenever we perform operations on the Slater determinants in our relevant bases. This signature of the permutation counts the minimum number of two particle exchanges that is requried to produce the final ordering. The Slater determinants span a complete orthonormal basis of the fermionic subspace of 𝐻 𝐴 , where ⟨𝑝𝑞𝑟 · · · | 𝑝𝑞𝑟 · · ·⟩ = 1, (2.8) when the indices of the bra are in the same ordering as the indices of the ket. Note that indices in bras or kets can always be permuted to match the ordering, at the cost of a phase. 2.3 Annihilation and Creation Operators In the context of second quantization, the annihilation and creation operators allow us to add or remove (respectively) occupied states in a particular Slater determinant. Interpreted a different way, the annihilation and creation operators connect states of differing number of particles, within the Fock space. We denote the annihilation operator by 𝑎 𝑝 and its Hermitian adjoint, the creation operator, as 𝑎 †𝑝 [6]. The operators are defined by their action on the Slater determinant, in the following way: 𝑎 𝑝 | 𝑝𝑞𝑟 · · ·⟩ = |𝑞𝑟 · · ·⟩ , 𝑎 𝑝 |𝑞𝑟 · · ·⟩ = 0, (2.9) 𝑎 †𝑝 |𝑞𝑟 · · ·⟩ = | 𝑝𝑞𝑟 · · ·⟩ , 𝑎 †𝑝 | 𝑝𝑞𝑟 · · ·⟩ = 0. (2.10) Here we show that the annihilation operator 𝑎 𝑝 removes a particle from the occupied single-particle state associated with 𝑝 in the target Slater determinant. Similarly, the creation operator 𝑎 †𝑝 adds a new particle to the single-particle state 𝑝. Attempting to annihilate a particle in a state which 10 is already unoccupied, or attempting to create a particle in a state which is already occupied, is forbidden. Note that this zero result is different from the vacuum state |0⟩, which is defined as the Slater determinant with no occupied single-particle states. It is clear that any application of an annihilator to this vacuum leads to a zero result, i.e. 𝑎 𝑝 |0⟩ = 0 = ⟨0| 𝑎 †𝑝 . (2.11) The fermionic creation and annihilation operators obey the following anticommutation relations [6], which will be important for our subsequent discussion: {𝑎 𝑝 , 𝑎 𝑞 } = {𝑎 †𝑝 , 𝑎 †𝑞 } = 0; {𝑎 †𝑝 , 𝑎 𝑞 } = 𝛿 𝑝𝑞 . (2.12) Note that any 𝐴-body Slater determinant can be constructed by a string of creation operators acting on the vacuum [2], Ö 𝐴 𝑎 †𝑝 𝑘 |0⟩ = | 𝑝 1 𝑝 2 · · · 𝑝 𝐴 ⟩ . (2.13) 𝑘=1 In this representation, we could compute the matrix element associated with two Slater determinants (sometimes called the overlap) via Ö 𝐵 Ö 𝐴 ⟨𝑞 1 𝑞 2 · · · 𝑞 𝐵 | 𝑝 1 𝑝 2 · · · 𝑝 𝐴 ⟩ = ⟨0| 𝑎 𝑞𝑙 𝑎 †𝑝 𝑘 |0⟩ . (2.14) 𝑙=1 𝑘=1 The result depends on the number of operators 𝐴 and 𝐵, as well as the ordering of indices. If 𝐵 = 𝐴 and indices 𝑝 𝑘 , 𝑞 𝑙 have matching order, then the result is 1. In this work, we will only concern ourselves with particle number conserving operators, so we can assume that 𝐵 = 𝐴 unless stated otherwise. In practical applications within nuclear theory, the vacuum state is an inefficient reference point from which to define the many-body basis; this is because of the resolution scale of the nuclear interaction [2]. More importantly, the vacuum state contains no information about the energy scales and occupation density of the nuclear matter we wish to describe. As a result, we may also define the creation and annihilation operators with respect to a particular reference state |Φ⟩ (called the Fermi vacuum for a Slater determinant reference state), which is often chosen as an approximation 11 to the ground state. In this case, we would define particle and hole operators, which are distinct in the kind of single particle state index on which they operate. The hole operators, usually indexed by 𝑖, 𝑗, 𝑘, . . . , act on states which are occupied in |Φ⟩. The particle operators, usually indexed by 𝑎, 𝑏, 𝑐, . . . , act on states which are unoccupied in |Φ⟩. We can construct excited Slater determinants by acting on the reference state with combinations of creation and annihilation operators, e.g. 𝑎 †𝑎 𝑎𝑖 |Φ⟩ = |Φ𝑖𝑎 ⟩ , (2.15) where the operator pair 𝑎 †𝑎 𝑎𝑖 annihilates a particle in the occupied (hole) state 𝑖 from the reference and creates a particle in the unoccupied (particle) state 𝑎. This action can be interpreted as a one-particle excitation relative to the reference state, and we may define a complete many-body basis as the set of all possible one-, two-, up to 𝐴-particle excitations for an 𝐴-body reference state. 2.4 Normal-Ordering and Wick Contractions The process of subtracting off contractions in the product of operators is called normal-ordering. The normal-ordered product of two operators 𝐴ˆ 𝐵ˆ reads   𝑛 𝐴ˆ 𝐵ˆ = 𝐴ˆ 𝐵ˆ − 𝐴𝐵, (2.16) where the contraction, in this context also called the Wick contraction, 𝐴𝐵 is defined as the expectation value of 𝐴ˆ 𝐵ˆ in the chosen reference state (e.g. the vacuum, in this case), 𝐴𝐵 ≡ ⟨0| 𝐴ˆ 𝐵|0⟩ ˆ . (2.17) By construction, this definition of the contraction leads to the consequence that the expectation   value of the normal-ordered operator 𝑛 𝐴ˆ 𝐵ˆ is zero, in the reference state for which the contraction was defined:   𝐴𝐵 = ⟨0| 𝐴ˆ 𝐵|0⟩ ˆ − ⟨0|𝑛 𝐴ˆ 𝐵ˆ |0⟩ (2.18)   =⇒ ⟨0|𝑛 𝐴ˆ 𝐵ˆ |0⟩ = 0 by Eq. 2.17. (2.19) These definitions extend to any length string of operator products, but will not be derived here (refer to Shavitt and Bartlett [6] and Hergert et al. [2]). These definitions also extend to any general reference state in which we define the Wick contraction. 12 We can exchange any two operators inside the normal-ordered product at the cost of a sign change, e.g.     𝑛 𝐴ˆ 𝐵ˆ = −𝑛 𝐵ˆ 𝐴ˆ . (2.20) Note that the contractions induced by permuting operators in the operator string have been taken care of by the expansion in terms generated by the normal-ordering process. Nuclear many-body theory primarily deals with pairs of annihilator/creator operators, e.g. 𝑎 †𝑝 𝑎 𝑞 or 𝑎 𝑝 𝑎 †𝑞 . The Wick contraction of operator pair 𝑎 𝑝 𝑎 †𝑞 may be written, by virtue of Equation 2.16, as 𝑎 𝑝 𝑎 †𝑞 = 𝑎 𝑝 𝑎 †𝑞 − 𝑛 𝑎 𝑝 𝑎 †𝑞 = {𝑎 𝑝 , 𝑎 †𝑞 } = 𝛿 𝑝𝑞 ,   (2.21) where we have written the only operator pair which gives a nonzero contraction, with respect to the vacuum. In this case, a particle in state 𝑞 is created in the vacuum and a particle in state 𝑝 is annihilated — thus, the contraction is 1 for 𝑝 = 𝑞 and zero otherwise. Interpreted a different way, this contraction is the expectation value of the operator 𝑎 𝑝 𝑎 †𝑞 in the vacuum. In a similar fashion, the Wick contraction of operator pairs with respect to the general reference is defined as 𝑎 𝑎 𝑎 †𝑏 = 𝛿𝑎𝑏 (2.22) 𝑎𝑖† 𝑎 𝑗 = 𝛿𝑖 𝑗 . We can interpret the results from the perspective of occupied and unoccupied states in the reference. The particle operator contraction is the expectation value ⟨Φ|𝑎 𝑎 𝑎 †𝑏 |Φ⟩, and is only nonzero when 𝑏 = 𝑎. The hole operator contraction is the expectation value ⟨Φ|𝑎𝑖† 𝑎 𝑗 |Φ⟩, and is only nonzero when 𝑖 = 𝑗. For a general reference state, the expectation value ⟨Φ|𝑎 †𝑝 𝑎 𝑞 |Φ⟩ defines the one-body density matrix: 𝜌 𝑞 𝑝 = ⟨Φ|𝑎 †𝑝 𝑎 𝑞 |Φ⟩ = 𝑎 †𝑝 𝑎 𝑞 ⟨Φ|Φ⟩ + ⟨Φ|𝑁 𝑎 †𝑝 𝑎 𝑞 |Φ⟩ .   (2.23) Note that the reference state is normalized so that ⟨Φ|Φ⟩ = 1. The definition of the Wick contraction h i implies the expectation value of the normal-product ⟨Φ|𝑁 𝑎 †𝑝 𝑎 𝑞 |Φ⟩ must vanish [2]. It follows 13 that the reference expectation value of any normal-product, with respect to that reference state, must vanish (analogous to Equation 2.19): ⟨Φ|𝑁 [· · ·] |Φ⟩ = 0. (2.24) Finally, note that we can drop the index distinction in Equation 2.22 altogether using the anticommutation relations in Equation 2.12. Starting with the definition of the contraction between hole operators and allowing the indices to span the entire single particle basis, we can define a new contraction 𝑎 𝑝 𝑎 †𝑞 which is related to 𝑎 †𝑝 𝑎 𝑞 in the following way: ⟨Φ|𝑎 †𝑝 𝑎 𝑞 |Φ⟩ = 𝑎 †𝑝 𝑎 𝑞 = 𝜌 𝑞 𝑝 (2.25) ⟨Φ|𝑎 𝑞 𝑎 †𝑝 |Φ⟩ = 𝑎 𝑞 𝑎 †𝑝 = ⟨Φ|(𝛿 𝑝𝑞 − 𝑎 †𝑝 𝑎 𝑞 )|Φ⟩ = 𝛿 𝑝𝑞 − 𝜌 𝑞 𝑝 . These two contractions automatically encode the behavior of expectation values of hole/particle pairs of operators in the reference state. For a single Slater determinant reference state, the elements of 𝑎 †𝑝 𝑎 𝑞 are nonzero only where both 𝑝, 𝑞 indices are occupied states in the reference. If either 𝑝 or 𝑞 (or both) are unoccupied states, the element is zero for one of two reasons: 1) the element vanishes due to Slater determinant orthogonality or 2) the action of a particle-index annihilator on the reference yields a zero result. Therefore, for a single Slater determinant reference state, the contraction 𝑎 𝑞 𝑎 †𝑝 yields elements which are 0 where the one-body density matrix is nonzero, and nonzero where the one-body density matrix is 0. For brevity and convenience, we usually work in the eigenbasis of the one-body density matrix: 𝜌 𝑝𝑞 = 𝑛( 𝑝)𝛿 𝑝𝑞 , (2.26) where the factor 𝑛( 𝑝) encodes the occupation of the single particle state 𝑝, i.e.,   if p is unoccupied in ref. (particle state)  0   𝑛( 𝑝) = (2.27)   1 if p in occupied in ref. (hole state)    Analogously to Equation 2.25, we define a particle occupation factor 𝑛( 𝑝), such that 𝑛( 𝑝) = 1 − 𝑛( 𝑝). (2.28) 14 These occupation factors are convenient for expressing occupation information in a general sin- gle particle basis because the diagonality of the one-body density 𝜌 eliminates summations in expansions of operators (e.g., refer to the IMSRG(2) flow equations later in Chapter 3, Equation 3.18). 2.5 Wick’s Theorem Whether with respect to the vacuum or a reference state, Wick’s theorem states that any (particle- number conserving) string of creation and annihilation operators can be written as the sum of all possible normal-product contractions [2, 1, 6]: WT 𝑎 †𝑝1 · · · 𝑎 †𝑝 𝐴 𝑎 𝑞 𝐴 · · · 𝑎 𝑞1 = 𝑁 𝑎 †𝑝1 · · · 𝑎 †𝑝 𝐴 𝑎 𝑞 𝐴 · · · 𝑎 𝑞1   + 𝑎 †𝑝1 𝑎 𝑞1 𝑁 𝑎 †𝑝2 · · · 𝑎 †𝑝 𝐴 𝑎 𝑞 𝐴 · · · 𝑎 𝑞2 − 𝑎 †𝑝1 𝑎 𝑞2 𝑁 𝑎 †𝑝2 · · · 𝑎 †𝑝 𝐴 𝑎 𝑞 𝐴 · · · 𝑎 𝑞1 + single contractions       † † † † + 𝑎 𝑝1 𝑎 𝑞1 𝑎 𝑝2 𝑎 𝑞2 − 𝑎 𝑝1 𝑎 𝑞2 𝑎 𝑝2 𝑎 𝑞1 𝑁 [· · ·] + double contractions + triple contractions + · · · + full contractions. (2.29) Notice that the process of removing the contraction from the normal-product results in an appropriate phase change on the normal-product operator string. The resulting sum of operators assumes the general form [2], 𝑀+𝑁 ∑︁ 𝐴𝑀 𝐵𝑁 = 𝐶 (𝑘) , (2.30) 𝑘=|𝑀−𝑁 | where the lower bounds exist because we are not able to fully contract all indices unless 𝑀 = 𝑁. The power of Wick’s theorem is that we can take any many-body operator product and “decom- pose” the product into a sum of normal-ordered products. In this expansion, we need only worry about phases when permuting the operators in the normal-ordered string, because contributions from the 𝛿 in the anticommutator have vanished according to the contractions. Additionally, when computing an expectation value with this theorem, only full contractions survive by Equation 2.24. We may leverage Wick’s theorem to define the two-body density matrix. Analogously to the one-body hole density matrix in Equation 2.23, the two-body hole density encode two-body 15 correlations in the reference state. We write, WT 𝜌𝑟 𝑠𝑝𝑞 = ⟨Φ|𝑎𝑟† 𝑎 †𝑠 𝑎 𝑞 𝑎 𝑝 |Φ⟩ = 𝜌𝑟 𝑝 𝜌 𝑠𝑞 − 𝜌𝑟𝑞 𝜌 𝑠𝑝 , (2.31) where we have used Wick’s theorem to further decompose the two-body operator string. Using Wick’s theorem, we have shown that the two-body hole density matrix, for a single Slater deter- minant reference state, can be reduced to a product of one-body density matrices. The three-body density matrix factorizes analogously, 𝜌 𝑠𝑡𝑢 𝑝𝑞𝑟 = ⟨Φ|𝑎 †𝑠 𝑎 𝑡† 𝑎 𝑢† 𝑎𝑟 𝑎 𝑞 𝑎 𝑝 |Φ⟩ WT = 𝜌 𝑠𝑝 𝜌𝑡𝑞 𝜌𝑢𝑟 − 𝜌 𝑠𝑞 𝜌𝑡 𝑝 𝜌𝑢𝑟 + 𝜌 𝑠𝑞 𝜌𝑡𝑟 𝜌𝑢 𝑝 (2.32) − 𝜌 𝑠𝑟 𝜌𝑡𝑞 𝜌𝑢 𝑝 − 𝜌 𝑠𝑝 𝜌𝑡𝑟 𝜌𝑢𝑞 + 𝜌 𝑠𝑟 𝜌𝑡 𝑝 𝜌𝑢𝑞 . 2.6 Correlated Reference States and Irreducible Density Matrices In Chapter 2.5, we discussed the one-body and two-body density matrices for a single, uncor- related Slater determinant reference state. For a correlated reference state, the density matrices include irreducible terms which encode “pure” correlations in the reference state. In order to define contractions with respect to a general reference state, Wick’s theorem must be augmented with additional contractions that reflect the present correlations, using a formalism developed by Kutzelnigg and Mukherjee [4] in the context of quantum chemistry. The new contractions are given by the so-called irreducible density matrices of the reference state, for reasons that will become clear soon. Following the notation of Hergert [1], we will denote these irreducible density matrices with the symbol 𝜆. In the correlated reference state, the one-body density matrix itself is irreducible: ⟨Φ|𝑎 †𝑝 𝑎 𝑞 |Φ⟩ = 𝜆 𝑞 𝑝 = 𝜌 𝑞 𝑝 . (2.33) The irreducible two-body density matrix is defined as:  𝜆𝑟 𝑠𝑝𝑞 = 𝜌𝑟 𝑠𝑝𝑞 − 𝜆𝑟 𝑝 𝜆 𝑠𝑞 − 𝜆𝑟𝑞 𝜆 𝑠𝑟 . (2.34) 16 Comparison with Equation 2.34 shows that for an uncorrelated reference state, such as a single Slater determinant, the irreducible two-body density matrix vanishes. In other words, 𝜆𝑖 𝑗 𝑘𝑙 represents pure two-body correlation information, after contributions from independent particles have been removed from the full two-body density matrix [1]. The irreducible three-body density matrix is similarly defined, 𝜆 𝑠𝑡𝑢 𝑝𝑞𝑟 = 𝜌 𝑠𝑡𝑢 𝑝𝑞𝑟 − A (𝜆 𝑠𝑝 𝜆𝑡𝑢𝑞𝑟 ) − A (𝜆 𝑠𝑝 𝜆𝑡𝑞 𝜆𝑢𝑟 ) (2.35) where A (·) represents all antisymmetric combinations of its argument (refer to [4, 5] for their complete derivation). In this case, the three-body irreducible term is the pure three-body correlation, after contributions from independent particles, as well as a correlated pair in presence of an independent spectator, have been subtracted from the full three-body density matrix. 2.7 Normal-ordering the Two-body Hamiltonian Consider a many-body vacuum Hamiltonian which consists of a one-body operator term, corresponding to single-particle kinetic energy, and a two-body operator term, ∑︁ 1 ∑︁ 𝐻= ⟨𝑝|𝑡|𝑞⟩ 𝑎 †𝑝 𝑎 𝑞 + ⟨𝑝𝑞|𝑣|𝑟 𝑠⟩ 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 , (2.36) 𝑝𝑞 4 𝑝𝑞𝑟 𝑠 where the vacuum two-body coefficients ⟨𝑝𝑞|𝑣|𝑟 𝑠⟩ are antisymmetrized. Wick’s theorem can be applied to 𝐻 to extract the normal-ordered 𝐻 𝑁 . We work in the eigenbasis of the one-body operator and begin by applying Wick’s theorem to the one-body operator: WT 𝑎 †𝑝 𝑎 𝑞 = 𝑁 𝑎 †𝑝 𝑎 𝑞 + 𝑁 𝑎 †𝑝 𝑎 𝑞 = 𝑁 𝑎 †𝑝 𝑎 𝑞 + 𝑛( 𝑝)𝛿 𝑝𝑞 .       (2.37) The Wick’s theorem two-body operator is written, WT 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 = 𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 +   𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 + 𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 +     (2.38) 𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 + 𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 +     𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 + 𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 ,     17 and evaluating the Wick contractions reveals the two-body dependence on the occupation structure in the reference state, WT 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 = 𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟   − 𝑛( 𝑝)𝛿 𝑝𝑠 𝑁 𝑎 †𝑞 𝑎𝑟 + 𝑛( 𝑝)𝛿 𝑝𝑟 𝑁 𝑎 †𝑞 𝑎 𝑠     (2.39) + 𝑛(𝑞)𝛿 𝑞𝑠 𝑁 𝑎 †𝑝 𝑎𝑟 − 𝑛(𝑞)𝛿 𝑞𝑟 𝑁 𝑎 †𝑝 𝑎 𝑠     + 𝑛( 𝑝)𝑛(𝑞)𝛿 𝑝𝑟 𝛿 𝑞𝑠 − 𝑛( 𝑝)𝑛(𝑞)𝛿 𝑝𝑠 𝛿 𝑞𝑟 . Note that the occupation factors 𝑛( 𝑝) and 𝑛(𝑞) restrict the single particle indices to only occupied states, or hole indices. Consequently, the summation indices in Equation 2.36 are also restricted where relevant. We proceed by collecting zero-, one-, and two-body terms into 0B, 1B, and 2B normal-ordered operators which appear on the right-hand side of Wick’s theorem applied to the vacuum Hamiltonian in Equation 2.36. Beginning with the 0B normal-ordered operator (fully contracted), ∑︁ 1 ∑︁ 0B = ⟨𝑖|𝑡|𝑖⟩ + ⟨𝑖 𝑗 |𝑣|𝑖 𝑗⟩ ≡ 𝐸, (2.40) 𝑖 2 𝑖𝑗 where the sum runs over the occupied (hole) states. This term is referred to as the reference energy 𝐸. The 1B normal-ordered operator (one contraction) reads ∑︁ ⟨𝑝|𝑡|𝑞⟩ 𝑁 𝑎 †𝑝 𝑎 𝑞   1B = 𝑝𝑞 1 ∑︁  †  + 𝑁 𝑎 𝑝 𝑎 𝑞 (⟨𝑖 𝑝|𝑣|𝑖𝑞⟩ + ⟨𝑝𝑖|𝑣|𝑞𝑖⟩ 4 𝑖 𝑝𝑞 − ⟨𝑖 𝑝|𝑣|𝑞𝑖⟩ − ⟨𝑝𝑖|𝑣|𝑖𝑞⟩) ∑︁ ∑︁  ⟨𝑝|𝑡|𝑞⟩ 𝑁 𝑎 †𝑝 𝑎 𝑞 𝑁 𝑎 †𝑝 𝑎 𝑞 ⟨𝑖 𝑝|𝑣|𝑖𝑞⟩    = + 𝑝𝑞 𝑖 𝑝𝑞 ∑︁ ⟨𝑝| 𝑓 |𝑞⟩ 𝑁 𝑎 †𝑝 𝑎 𝑞 ,   ≡ 𝑝𝑞 where the Fock operator 𝑓 contains the kinetic energy of the particle plus a mean-field contribution from the interaction. 18 The 2B normal-ordered operator is simply the last remaining two-body term which contains no contractions, 1 ∑︁ ⟨𝑝𝑞|𝑣|𝑟 𝑠⟩ 𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 .   2B = (2.41) 4 𝑝𝑞𝑟 𝑠 Finally, the normal-ordered two-body Hamiltonian reads ∑︁  1 ∑︁ ⟨𝑝| 𝑓 |𝑞⟩ 𝑁 𝑎 †𝑝 𝑎 𝑞 + ⟨𝑝𝑞|𝑣|𝑟 𝑠⟩ 𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟    𝐻𝑁 = 𝐸 + (2.42) 𝑝𝑞 4 𝑝𝑞𝑟 𝑠 Wick’s theorem shuffles operator contributions into terms of fewer-body character by contracting with the density matrix of the independent particle system. There are situations where we may want to reverse the normal-ordering process in order to recover the true vacuum operators. To see how we can undo the normal-ordering of the Hamiltonian (or any other operator), we compare coefficients in Equation 2.42 (the normal-ordered Hamiltonian) and Equation 2.36 (the two-body vacuum Hamiltonian). Comparing coefficients where vac corresponds to vacuum operators and NO corresponds to normal-ordered operators, ⟨𝑝𝑞|𝑣|𝑟 𝑠⟩ vac = ⟨𝑝𝑞|𝑣|𝑟 𝑠⟩ NO ∑︁ ⟨𝑝|𝑡|𝑞⟩ vac = ⟨𝑝| 𝑓 |𝑞⟩ NO − ⟨𝑝𝑖|𝑣|𝑞𝑖⟩ vac (2.43) 𝑖 ∑︁ 1 ∑︁ 𝐻0 = 𝐸 − ⟨𝑖|𝑡|𝑖⟩ vac − ⟨𝑖 𝑗 |𝑣|𝑖 𝑗⟩ vac , 𝑖 4 𝑖𝑗 where the 𝐻0 is necessary to recover information lost due to truncation of higher-body operators generated by the normal-ordering process. In the case where no operators were truncated in the normal-ordered operator, 𝐻0 vanishes by Equation 2.40. If normal-ordered contributions which would have been shuffled into the zero-body term in Equation 2.40 are instead truncated off, then the unnormal-ordered 𝐻0 piece will not vanish, and must be included in the total vacuum Hamiltonian, written 𝐻 vac = 𝐻0 + 𝐻1B vac vac + 𝐻2B . (2.44) This procedure is useful in the context of IMSRG-evolved Hamiltonians (Chapter 3.3), where we may wish to use vacuum coefficients obtained from the unnormal-ordered evolved Hamiltonian as inputs for exact methods (e.g. Full Configuration Interaction calculations in Chapter 2.10). 19 2.8 Density Matrix Normal-ordering When generally correlated reference states are introduced to the normal-ordering, we must mod- ify the normal-ordering procedure to account for occupations which are fractional; in a correlated reference, the occupation factor in Equations 2.27 and 2.28 is no longer an integer, in general. Normal-ordering with respect to a general density matrix is as simple as inserting the full density matrices, where appropriate, into the calculation of the normal-ordered operator in Equations 2.40, 2.41, and 2.41. Summations over hole indices must then be replaced by sums over the full single- particle basis instead, e.g. ∑︁ ∑︁ ⟨𝑝𝑖|𝑣|𝑞𝑖⟩ vac −→ ⟨𝑝𝑟 |𝑣|𝑞𝑠⟩ vac 𝜌 𝑠𝑟 𝑖 𝑟𝑠 ∑︁ ∑︁ ⟨𝑖|𝑡|𝑖⟩ vac −→ ⟨𝑝|𝑡|𝑞⟩ vac 𝜌 𝑞 𝑝 (2.45) 𝑖 𝑟𝑠 1 ∑︁ 1 ∑︁ ⟨𝑖 𝑗 |𝑣|𝑖 𝑗⟩ vac −→ ⟨𝑝𝑞|𝑣|𝑟 𝑠⟩ vac 𝜌𝑟 𝑠𝑝𝑞 . 4 𝑖𝑗 4 𝑝𝑞𝑟 𝑠 2.9 Slater Determinant Occupation Number Representation A common strategy for writing and computing Slater determinants, particularly in code, is the occupation number representation. In this scheme, a Slater determinant with number of single states 𝑁 is written as a string of 1s and 0s of length 𝑁, where 1 indicates an occupied single particle state and 0 an unoccupied single particle state. For example, a Fermi vacuum Slater determinant for 𝑁 = 8 states and 𝐴 = 4 particles states can be written as |Φ0 ⟩ = |11110000⟩ , (2.46) where the ordering of the bit digits is convention, but must be consistent within a given application. In this scheme, the application of creation/annihilation operators replace 0s with 1s and 1s with 0s where appropriate, and the usual rules for antisymmetric exchanges of indices still apply. For example, we can apply an excitation operator 𝐶ˆ45 = 𝑎 †5 𝑎 4 which excites the particle occupied in single particle index 𝑖 = 4 into the single particle index 𝑎 = 5 above the Fermi surface: 𝐶ˆ45 |11110000⟩ = − |11101000⟩ . (2.47) 20 The Slater determinants in this scheme are still orthonormal, meaning that ⟨11110000|11110000⟩ = 1 (2.48) ⟨11110000|𝐶ˆ45 |11110000⟩⟩ = − ⟨11110000|11101000⟩ = 0. 2.10 Full Configuration Interaction Theory to Solve Model Hamiltonians Sometimes we require an exact method in order to validate results from approximate models. The exact method we use often in this dissertation is called Full Configuration Interaction (FCI) theory. The defining feature of an FCI calculation is that the Hamiltonian is constructed in the full basis of Slater determinant configurations of particles [3]. In the occupation number representation, we proceed by constructing the Hamiltonian matrix in the full basis of configurations of occupation strings relative to a “ground state” configuration such as the occupation string in Equation 2.46. After the matrix is built, we can apply a diagonalization method to extract the energy spectrum, or other target observables. This method is prohibitive for practical problems in nuclear theory due to the dimensionality of the model space. For a model with 𝐴 particles in 𝑁 possible single particle states, the total number of Slater determinant configurations is given by [3]   𝑁 𝑁! = , (2.49) 𝐴 (𝑁 − 𝐴)!𝐴! This number quickly explodes with the number of possible single particle states, which may be on the order of 100s or 1000s for typical medium and heavy nuclei. For a toy model, such as the pairing model (see Chapter 4.1), the model space may be so simplified that the dimensionality is not an issue. In the pairing model, we can typically keep our particle number small (on the order of 10 or fewer), and fill half the number of energy levels, in order to extract the physics we require using our theoretical methods. An FCI calculation involving a basis of Slater determinant configurations is performed over several steps: 1. Define the number of single particle states 𝑛 and number of particles 𝑁 21 𝑛 2. Generate 𝑁 number of Slater determinants in the occupation number representation as the basis (see Chapter 2.9) 𝑛 𝑛 3. Compute 𝑁 × 𝑁 square matrix H with elements given by ⟨Φ𝑖 |𝐻|Φ 𝑗 ⟩, where Φ𝑖 is the 𝑖-th Slater determinant in the basis (typically ordered by increasing energy with respect to the ground state, although this ordering may not be known a priori) 4. Diagonalize H to extract the energy spectrum Since we have access to the eigenstates after the diagonalization, we can in principle compute any other observable of interest in this framework. When calculations are feasible, the energy spectrum extracted from FCI is the ultimate benchmark for validating results from other techniques, which will be presented in later sections. 22 BIBLIOGRAPHY [1] H Hergert. “In-medium similarity renormalization group for closed and open-shell nu- clei”. In: Physica Scripta 92.2 (Dec. 2016), p. 023002. issn: 1402-4896. doi: 10.1088/1402-4896/92/2/023002. url: http://dx.doi.org/10.1088/1402-4896/92/2/023002. [2] H. Hergert et al. “The In-Medium Similarity Renormalization Group: A novel ab initio method for nuclei”. In: Physics Reports 621 (Mar. 2016), pp. 165–222. issn: 0370- 1573. doi: 10.1016/j.physrep.2015.12.007. url: http://dx.doi.org/10.1016/j.physrep.2015. 12.007. [3] Morten Hjorth-Jensen, Maria-Paola Lombardo, and Ubirajara van Kolck. An advanced course in computational nuclear physics: bridging the scales from quarks to neutron stars. Springer, 2017. [4] Werner Kutzelnigg and Debashis Mukherjee. “Normal order and extended Wick theorem for a multiconfiguration reference wave function”. In: The Journal of Chemical Physics 107.2 (1997), pp. 432–449. doi: 10.1063/1.474405. eprint: https://doi.org/10.1063/1.474405. url: https://doi.org/10.1063/1.474405. [5] Debashis Mukherjee and Werner Kutzelnigg. “Irreducible Brillouin conditions and con- tracted Schrödinger equations for n-electron systems. I. The equations satisfied by the density cumulants”. In: The Journal of Chemical Physics 114.5 (2001), pp. 2047– 2061. doi: 10.1063/1.1337058. eprint: https://doi.org/10.1063/1.1337058. url: https://doi.org/10.1063/1.1337058. [6] I. Shavitt and R.J. Bartlett. Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory. Cambridge Molecular Science. Cambridge University Press, 2009. isbn: 9780521818322. url: https://books.google.com/books?id=SWw6ac1NHZYC. 23 CHAPTER 3 IN-MEDIUM SIMILARITY RENORMALIZATION GROUP 3.1 Introduction In this chapter, we will present the conceptual and mathematical framework for the in-medium similarity renormalization group (IMSRG), which is the main focus for this body of work. We begin by presenting the idea behind Similarity Renormalization Group methods. 3.2 Similarity Renormalization Group The SRG was conceived by Głazek and Wilson [2] for applications in light-front quantum field theory, and independently developed by Wegner [8] for applications in condensed matter systems. The SRG is constructed such that the Hamiltonian matrix, expressed in an appropriate basis, is driven to a band-diagonal form via a continuous unitary transformation that suppresses off-diagonal elements. In a basis like the harmonic oscillator, these off-diagonal elements encode interactions between states of increasing difference in energy. This process is demonstrated schematically in Figure 3.1. The suppression scales 𝑠0 , 𝑠1 , 𝑠2 represent characteristic relative momentum “bands” which limit relative momentum according to 𝑠 = 𝜆1/4 , where 𝜆 is in units fm−1 . The scale parameter 𝑠 parameterizes the cutoff Λ in the 𝜒EFT theory (see Chapter 4.2). As 𝑠 increases, high momentum modes are decoupled from low momentum modes. The flowing SRG Hamiltonian is expressed as 𝐻 (𝑠) = 𝑈 (𝑠)𝐻 (𝑠 = 0)𝑈 † (𝑠), (3.1) where 𝑈 (𝑠) describes the evolution of the initial, unevolved Hamiltonian to the current value of the continuous flow parameter 𝑠. 𝑑𝐻 (𝑠) 𝑑𝑈 (𝑠) † 𝑑𝑈 † (𝑠) = 𝐻 (0)𝑈 (𝑠) + 𝑈 (𝑠)𝐻 (0) 𝑑𝑠 𝑑𝑠 𝑑𝑠 (3.2) 𝑑𝑈 (𝑠) † 𝑑𝑈 † (𝑠) = 𝑈 (𝑠)𝐻 (𝑠) + 𝐻 (𝑠)𝑈 (𝑠) . 𝑑𝑠 𝑑𝑠 24 q q’ 𝑠2 𝑠1 𝑠0 Figure 3.1 Schematic diagram of 𝑁 𝑁 scattering potential, in a momentum basis, showing how modes of increasing relative momentum (energy) are suppressed faster than low momentum modes. Figure adapted from Bogner, Furnstahl, and Schwenk [1]. Next, we define the 𝑠-dependent anti-Hermitian generator 𝜂(𝑠), which is chosen in a way to “gen- erate” the transformation we want to achieve (i.e. suppression of off-diagonal elements according to chosen decoupling strategy). The generator 𝜂(𝑠) is defined as [4, 3, 5, 1] 𝑑𝑈 (𝑠) † 𝜂(𝑠) = 𝑈 (𝑠) = −𝜂† (𝑠). (3.3) 𝑑𝑠 Inserting Equation 3.3 into Equation 3.2, we naturally arrive at the commutator expression for the operator flow equation: 𝑑𝐻 (𝑠) = [𝜂(𝑠), 𝐻 (𝑠)] . (3.4) 𝑑𝑠 3.3 IMSRG Formalism The computational cost of the SRG, with matrices, in Equation 3.4 scales exponentially with the many-body basis dimension. The IMSRG approach moves this description “in-medium” by shifting the flowing Hamiltonian to a normal-ordered operator basis (Chapter 2), instead of the many-body basis [4]. The flow is implemented at the level of operators, using strings of annhilation and 25 creation operators that are normal-ordered with respect to an appropriate reference state, defining our operator basis. We will see that the IMSRG flow equations are polynomially scaling with single particle basis size. Practically, the purpose of the IMSRG is to decouple a target state from excitations relative to that state, by suppressing the interaction terms in the normal-ordered operator expansion of the Hamiltonian. We start from Equation 3.4, and use the second quantization techniques discussed in Chapter 2.4 to derive the coupled flow equations associated with the IMSRG evolution. We assume that the operator 𝐻 (𝑠) is in normal-ordered form with respect to a chosen reference state (which we might choose to approximate the Hamiltonian ground state); the normal-ordered two-body Hamiltonian is written, ∑︁  1 ∑︁ 𝑓 𝑝𝑞 (𝑠)𝑁 𝑎 †𝑝 𝑎 𝑞 + Γ𝑝𝑞𝑟 𝑠 (𝑠)𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 ,    𝐻IMSRG (𝑠) = 𝐸 (𝑠) + (3.5) 𝑝𝑞 4 𝑝𝑞𝑟 𝑠 where the 𝑓 𝑝𝑞 and Γ𝑝𝑞𝑟 𝑠 operator coefficients are derived from the vacuum representation of the Hamiltonian through the normal ordering procedure, and 𝐸 represents the reference energy, which encodes the expectation value of the vacuum Hamiltonian in the reference state (see Chapter 2.7). We assume that Γ𝑝𝑞𝑟 𝑠 is antisymmetrized, i.e., Γ𝑝𝑞𝑟 𝑠 = −Γ𝑞 𝑝𝑟 𝑠 = −Γ𝑝𝑞𝑠𝑟 = Γ𝑞 𝑝𝑠𝑟 . (3.6) Note that all operator coefficients in the IMSRG Hamiltonian (3.5) depend on the flow parameters 𝑠. In this way, the IMSRG evolves the Hamiltonian within the chosen operator basis toward a desired form which is defined by a chosen decoupling scheme. Evaluation of the commutators in Equation 3.4 induces higher-body terms that are computationally expensive to calculate. In general, we need to truncate the operator basis at a particular level. For example, keeping up to two-body operators throughout the flow defines the IMSRG(2) truncation. In the IMSRG(2) truncation, the anti-Hermitian generator 𝜂(𝑠) has the form, ∑︁  1 ∑︁ 𝜂 𝑝𝑞 (𝑠)𝑁 𝑎 †𝑝 𝑎 𝑞 + 𝜂 𝑝𝑞𝑟 𝑠 (𝑠)𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 .    𝜂(𝑠) = (3.7) 𝑝𝑞 4 𝑝𝑞𝑟 𝑠 26 To proceed, we recall the general form of a product of operators from Equation 2.30. Each product in the commutator of Equation 3.4 generates new operators up to 4-body; however, the symmetry in the commutator eliminates the four-body operators. Again, working in the IMSRG(2) truncation means we keep only up to two-body operators generated by the commutators. h h i h ii As a demonstration, we will derive the term associated with 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 , 𝑁 𝑎 †𝑝2 𝑎 𝑞2 , where we have added the sub-indices 1 and 2 to differentiate between the source of the operator (Hamiltonian versus generator). The logic we use here can be extended to each additional commu- tator term to completely derive the full IMSRG flow equations. We start with the definition of the commutator, 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 , 𝑁 𝑎 †𝑝2 𝑎 𝑞2 =𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 𝑁 𝑎 †𝑝2 𝑎 𝑞2          (3.8) − 𝑁 𝑎 †𝑝2 𝑎 𝑞2 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1     We may apply Wick’s theorem (refer Shavitt and Bartlett [7], Equation 3.194, as well as Chapter 2.5) to a product of normal-ordered operators, which operates similarly to application of Wick’s Í theorem to a vacuum operator. Using the first term as an example (dropping the summations for brevity),  WT  𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 𝑁 𝑎 †𝑝2 𝑎 𝑞2 = 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 𝑎 †𝑝2 𝑎 𝑞2     h i h i † † † + 𝑁 𝑎 𝑝 1 𝑎 𝑞 1 𝑎 𝑠1 𝑎 𝑟 1 𝑎 𝑝 2 𝑎 𝑞 2 + 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 𝑎 †𝑝2 𝑎 𝑞2 h i h i † † † + 𝑁 𝑎 𝑝 1 𝑎 𝑞 1 𝑎 𝑠1 𝑎 𝑟 1 𝑎 𝑝 2 𝑎 𝑞 2 + 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 𝑎 †𝑝2 𝑎 𝑞2 (3.9)     † † † + 𝑁 𝑎 𝑝 1 𝑎 𝑞 1 𝑎 𝑠1 𝑎 𝑟 1 𝑎 𝑝 2 𝑎 𝑞 2 + 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 𝑎 †𝑝2 𝑎 𝑞2     † † † + 𝑁 𝑎 𝑝 1 𝑎 𝑞 1 𝑎 𝑠1 𝑎 𝑟 1 𝑎 𝑝 2 𝑎 𝑞 2 + 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 𝑎 †𝑝2 𝑎 𝑞2 Contractions between the operators inside a particular normal-ordered operator string on the left- hand side of Equation 3.9 are already accounted for by the initial normal ordering, hence only contractions between operators from string 1 and 2 survive. Simplifying Equation 3.9 further using the Wick contractions in Equations 2.22 and 2.25, we 27 find 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 𝑁 𝑎 †𝑝2 𝑎 𝑞2 = 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 𝑎 †𝑝2 𝑎 𝑞2       + 𝜌 𝑝1 𝑞2 𝑁 𝑎 †𝑞1 𝑎 †𝑝2 𝑎 𝑠1 𝑎𝑟1 − 𝜌 𝑞1 𝑞2 𝑁 𝑎 †𝑝1 𝑎 †𝑝2 𝑎 𝑠1 𝑎𝑟1     − (𝛿 𝑠1 𝑝2 − 𝜌 𝑠1 𝑝2 )𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎𝑟1 𝑎 𝑞2 + (𝛿𝑟1 𝑝2 − 𝜌𝑟1 𝑝2 )𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎 𝑞2 (3.10)     − 𝜌 𝑝1 𝑞2 (𝛿 𝑠1 𝑝2 − 𝜌 𝑠1 𝑝2 )𝑁 𝑎 †𝑞1 𝑎𝑟1 + 𝜌 𝑝1 𝑞2 (𝛿𝑟1 𝑝2 − 𝜌𝑟1 𝑝2 )𝑁 𝑎 †𝑞1 𝑎 𝑠1     + 𝜌 𝑞1 𝑞2 (𝛿 𝑠1 𝑝2 − 𝜌 𝑠1 𝑝2 )𝑁 𝑎 †𝑝1 𝑎𝑟1 − 𝜌 𝑞1 𝑞2 (𝛿𝑟1 𝑝2 − 𝜌𝑟1 𝑝2 )𝑁 𝑎 †𝑝1 𝑎 𝑠1 .     h i h i The same procedure applied to the product 𝑁 𝑎 †𝑝2 𝑎 𝑞2 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 yields a similar result, 𝑁 𝑎 †𝑝2 𝑎 𝑞2 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 = 𝑁 𝑎 †𝑝2 𝑎 𝑞2 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1       − 𝜌 𝑝2𝑟1 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎 𝑞2 + 𝜌 𝑝2 𝑠1 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎𝑟1 𝑎 𝑞2     − (𝛿 𝑞2 𝑝1 − 𝜌 𝑞2 𝑝1 )𝑁 𝑎 †𝑞1 𝑎 †𝑝2 𝑎 𝑠1 𝑎𝑟1 + (𝛿 𝑞2 𝑞1 − 𝜌 𝑞2 𝑞1 )𝑁 𝑎 †𝑝1 𝑎 †𝑝2 𝑎 𝑠1 𝑎𝑟1     (3.11) − 𝜌 𝑝2 𝑠1 (𝛿 𝑞2 𝑝1 − 𝜌 𝑞2 𝑝1 )𝑁 𝑎 †𝑞1 𝑎𝑟1 + 𝜌 𝑝2𝑟1 (𝛿 𝑞2 𝑝1 − 𝜌 𝑞2 𝑝1 )𝑁 𝑎 †𝑞1 𝑎 𝑠1     + 𝜌 𝑝2 𝑠1 (𝛿 𝑞2 𝑞1 − 𝜌 𝑞2 𝑞1 )𝑁 𝑎 †𝑝1 𝑎𝑟1 − 𝜌 𝑝2𝑟1 (𝛿 𝑞2 𝑞1 − 𝜌 𝑞2 𝑞1 )𝑁 𝑎 †𝑝1 𝑎 𝑠1 ,     where we have exchanged operators inside the normal-product to match Equation 3.10. At this point, the terms which vanish from the commutator are straightforward to recognize—the three- body normal-product vanishes, as well as any term which only includes the one-body density matrix 𝜌. The commutator becomes 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎𝑟1 , 𝑁 𝑎 †𝑝2 𝑎 𝑞2 =      − 𝛿 𝑠1 𝑝2 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎𝑟1 𝑎 𝑞2 + 𝛿𝑟1 𝑞2 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎 𝑞2     + 𝛿 𝑞2 𝑝1 𝑁 𝑎 †𝑞1 𝑎 †𝑝2 𝑎 𝑠1 𝑎𝑟1 − 𝛿 𝑞2 𝑞1 𝑁 𝑎 †𝑝1 𝑎 †𝑝2 𝑎 𝑠1 𝑎𝑟1     + −𝜌 𝑝1 𝑞2 𝛿 𝑠1 𝑝2 + 𝜌 𝑝2 𝑠1 𝛿 𝑞2 𝑝1 𝑁 𝑎 †𝑞1 𝑎𝑟1    (3.12) + 𝜌 𝑝1 𝑞2 𝛿𝑟1 𝑝2 − 𝜌 𝑝2𝑟1 𝛿 𝑞2 𝑝1 𝑁 𝑎 †𝑞1 𝑎 𝑠1    + 𝜌 𝑞1 𝑞2 𝛿 𝑠1 𝑝2 − 𝜌 𝑝2 𝑠1 𝛿 𝑞2 𝑞1 𝑁 𝑎 †𝑝1 𝑎𝑟1    + −𝜌 𝑞1 𝑞2 𝛿𝑟1 𝑝2 + 𝜌 𝑝2𝑟1 𝛿 𝑞2 𝑞1 𝑁 𝑎 †𝑝1 𝑎 𝑠1 .    Í Now we will re-establish the summation over all single particle states for each operator in the commutator. Since this sum covers all single particle indices for every operator in each normal- 28 product, we can rename and manipulate indices to our advantage. Inside the normal-product, the label of a particular operator is “free” because the indices span the full space. Introducing the appropriate matrix elements associated with each operator reveals that we can combine each term into an associated many-body operator, multiplied by a corresponding matrix element (or expression). Applying this logic, h i 𝐴 (2) , 𝐵 (1) = 1 ∑︁  𝐴 𝑝1 𝑞1𝑟1 𝑠1 𝐵 𝑝2 𝑞2 −𝛿 𝑠1 𝑝2 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎𝑟1 𝑎 𝑞2 + 𝛿𝑟1 𝑝2 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎 𝑞2     4 𝑝 1 𝑞 1 𝑟 1 𝑠1 𝑝 2 𝑞 2 + 𝛿 𝑞2 𝑝1 𝑁 𝑎 †𝑞1 𝑎 †𝑝2 𝑎 𝑠1 𝑎𝑟1 − 𝛿 𝑞2 𝑞1 𝑁 𝑎 †𝑝1 𝑎 †𝑝2 𝑎 𝑠1 𝑎𝑟1     (3.13) + −𝜌 𝑝1 𝑞2 𝛿 𝑠1 𝑝2 + 𝜌 𝑝2 𝑠1 𝛿 𝑞2 𝑝1 𝑁 𝑎 †𝑞1 𝑎𝑟1    + 𝜌 𝑝1 𝑞2 𝛿𝑟1 𝑞2 − 𝜌 𝑝2𝑟1 𝛿 𝑞2 𝑝1 𝑁 𝑎 †𝑞1 𝑎 𝑠1    + 𝜌 𝑞1 𝑞2 𝛿 𝑠1 𝑝2 − 𝜌 𝑝2 𝑠1 𝛿 𝑞2 𝑞1 𝑁 𝑎 †𝑝1 𝑎𝑟1      †  + −𝜌 𝑞1 𝑞2 𝛿𝑟1 𝑝2 + 𝜌 𝑝2𝑟1 𝛿 𝑞2 𝑞1 𝑁 𝑎 𝑝1 𝑎 𝑠1 . For simplicity, we will continue in the eigenbasis of the one-body density matrix such that, 𝜌 𝑝𝑞 = 𝑛 𝑝 𝛿 𝑝𝑞 , (3.14) where the factor 𝑛 𝑝 is 0 or 1 when the single particle state 𝑎 is unoccupied or occupied, respectively, in the reference [4]. Note that working in this basis, we assume that the reference state is a single Slater determinant. We push the associated matrix element 𝐴 𝑝1 𝑞1𝑟1 𝑠1 𝐵 𝑝2 𝑞2 through each term in the sum, carefully 29 removing index dependence that is cancelled by a particular 𝛿, h i 𝐴 (2) , 𝐵 (1) = 1 ∑︁  1 ∑︁ 𝐴 𝑝1 𝑞1𝑟1 𝑢 𝐵𝑢𝑞2 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎𝑟1 𝑎 𝑞2 + 𝐴 𝑝1 𝑞1 𝑢𝑠1 𝐵𝑢𝑞2 𝑁 𝑎 †𝑝1 𝑎 †𝑞1 𝑎 𝑠1 𝑎 𝑞2    − 4 𝑝 𝑞 𝑟 𝑢𝑞 4 𝑝 𝑞 𝑢𝑠 𝑞 1 1 1 2 1 1 1 2 1 ∑︁ 1 ∑︁ 𝐴𝑢𝑞1𝑟1 𝑠1 𝐵 𝑝2 𝑢 𝑁 𝑎 †𝑞1 𝑎 †𝑝2 𝑎 𝑠1 𝑎𝑟1 − 𝐴 𝑝1 𝑢𝑟1 𝑠1 𝐵 𝑝2 𝑢 𝑁 𝑎 †𝑝1 𝑎 †𝑝2 𝑎 𝑠1 𝑎𝑟1     + 4𝑞𝑟 𝑠 𝑝𝑢 4𝑝𝑟𝑠𝑝𝑢 1 1 1 2 1 1 1 2 1 ∑︁ 𝐴𝑣𝑞1𝑟1 𝑢 𝐵𝑢𝑣 (−𝑛𝑣 + 𝑛𝑢 ) 𝑁 𝑎 †𝑞1 𝑎𝑟1   + 4 𝑞 𝑟 𝑢𝑣 (3.15) 1 1 1 ∑︁ 𝐴𝑣𝑞1 𝑢𝑠1 𝐵𝑢𝑣 (𝑛𝑣 − 𝑛𝑢 ) 𝑁 𝑎 †𝑞1 𝑎 𝑠1   + 4 𝑞 𝑠 𝑢𝑣 1 1 1 ∑︁ 𝐴 𝑝1 𝑣𝑟1 𝑢 𝐵𝑢𝑣 (𝑛𝑣 − 𝑛𝑢 ) 𝑁 𝑎 †𝑝1 𝑎𝑟1   + 4 𝑝 𝑟 𝑢𝑣 1 1 1 ∑︁ 𝐴 𝑝1 𝑢𝑣𝑠1 𝐵𝑢𝑣 (−𝑛𝑣 + 𝑛𝑢 ) 𝑁 𝑎 †𝑝1 𝑎 𝑠1 ,   + 4 𝑝 𝑠 𝑢𝑣 1 1 where the new summation indices 𝑢, 𝑣 result from the associated 𝛿 restrictions. The four one-body operator terms in Equation 3.15 are equivalent up to a sign, and the associated indices may be permuted in a way to combine all four terms into a single one-body operator term that reads h i (1) ∑︁ (2) (1) 𝐴𝑣 𝑝𝑢𝑞 𝐵𝑢𝑣 (𝑛𝑣 − 𝑛 𝑏 ) 𝑁 𝑎 †𝑝 𝑎 𝑞 .   𝐴 ,𝐵 = (3.16) 𝑝𝑞𝑢𝑣 The four two-body terms can be rearranged into a single two-body operator with a coefficient that consists of a sum over permuted matrix elements, i.e. h i (2) 1 ∑︁  † † (2) (1)  𝐴 ,𝐵 = 𝑁 𝑎 𝑝 𝑎 𝑞 𝑎 𝑠 𝑎𝑟 × 4 𝑝𝑞𝑟 𝑠𝑢 𝐴 𝑝𝑞𝑟𝑢 𝐵𝑢𝑠 + 𝐴 𝑝𝑞𝑢𝑠 𝐵𝑢𝑟  −𝐴𝑢𝑞𝑟 𝑠 𝐵 𝑝𝑢 − 𝐴 𝑝𝑢𝑟 𝑠 𝐵𝑞𝑢 (3.17) 1 ∑︁ 𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 ×   = 4 𝑝𝑞𝑟 𝑠𝑢  1 − 𝑃𝑟 𝑠 ) 𝐴 𝑝𝑞𝑢𝑠 𝐵𝑢𝑟 − (1 − 𝑃 𝑝𝑞 ) 𝐴𝑢𝑞𝑟 𝑠 𝐵 𝑝𝑢 Here, the operator 𝑃 𝑝𝑞 exchanges the indices 𝑝, 𝑞. We have arrived at the final result for just one commutator that appears in the IMSRG(2) flow equations, with seemingly Herculean effort. 30 Computing additional commutators requires similar, or greater, effort. Luckily, the fundamental commutators which can be used to compute the remaining terms in the flow equation commutator are catalogued in Appendix A of Hergert et al. [4], and tools are now available that can automate these derivations. The zero-, one-, and two-body terms that are produced by the commutators can now be collected into separate, but coupled, flow equations for the reference energy 𝐸, the one-body operator coefficients 𝑓 , and the two-body operator coefficients Γ. The resulting system of IMSRG(2) flow equations reads[4] 𝑑𝐸 ∑︁ = (𝑛 𝑝 − 𝑛𝑞 )𝜂 𝑝𝑞 𝑓𝑞 𝑝 𝑑𝑠 𝑝𝑞 1 ∑︁ + 𝜂 𝑝𝑞𝑟 𝑠 Γ𝑟 𝑠𝑝𝑞 𝑛 𝑝 𝑛𝑞 𝑛¯ 𝑟 𝑛¯ 𝑠 2 𝑝𝑞𝑟 𝑠 𝑑𝑓 𝑝𝑞 ∑︁ = (1 + 𝑃 𝑝𝑞 )𝜂 𝑝𝑟 𝑓𝑟𝑞 𝑑𝑠 𝑟 ∑︁ + (𝑛𝑟 − 𝑛 𝑠 )(𝜂𝑟 𝑠 Γ𝑠𝑝𝑟𝑞 − 𝑓𝑟 𝑠 𝜂 𝑠𝑝𝑟𝑞 ) 𝑟𝑠 1 ∑︁ (3.18) + (𝑛𝑟 𝑛 𝑠 𝑛¯ 𝑡 + 𝑛¯ 𝑟 𝑛¯ 𝑠 𝑛𝑡 )(1 + 𝑃 𝑝𝑞 )𝜂𝑡 𝑝𝑟 𝑠 Γ𝑟 𝑠𝑡𝑞 2 𝑟 𝑠𝑡 𝑑Γ𝑝𝑞𝑟 𝑠 ∑︁ = {(1 − 𝑃 𝑝𝑞 )(𝜂 𝑝𝑡 Γ𝑡𝑞𝑟 𝑠 − 𝑓 𝑝𝑡 𝜂𝑡𝑞𝑟 𝑠 𝑑𝑠 𝑡 − (1 − 𝑃𝑟 𝑠 )(Γ𝑝𝑞𝑡𝑠 𝜂𝑡𝑟 − 𝜂 𝑝𝑞𝑡𝑠 𝑓𝑡𝑟 ))} 1 ∑︁ + (1 − 𝑛𝑡 − 𝑛𝑢 )(𝜂 𝑝𝑞𝑡𝑢 Γ𝑡𝑢𝑟 𝑠 − Γ𝑝𝑞𝑡𝑢 𝜂𝑡𝑢𝑟 𝑠 ) 2 𝑡𝑢 ∑︁ − (𝑛𝑡 − 𝑛𝑢 )(1 − 𝑃 𝑝𝑞 )(1 − 𝑃𝑟 𝑠 )𝜂𝑢𝑞𝑡𝑠 Γ𝑡 𝑝𝑢𝑟 , 𝑡𝑢 where indices run over the full single-particle basis. Note that the occupation numbers 𝑛 𝑝 do not depend on 𝑠 because the reference state is held fixed during the evolution. Since we must loop over all operator coefficients for each flow, the computational effort for implementing the IMSRG(2) evolution naively scales like 𝑂 (𝑁 6 ) where 𝑁 is the number of single particle states. For Slater determinant references and certain types of IMSRG generators, we can leverage the partitioning of the single-particle basis into particle and hole states to achieve an 𝑂 (𝑁 ℎ2 𝑁 𝑝4 ) scaling akin to the Coupled Cluster method with Singles and Doubles excitations 31 a ab abc a ab abc | | i | ij | ijk | | i | ij | ijk U(s)HU† (s) | | a i a i | | ab ab ij ij | | abc s abc ijk ijk | | H H( ) Figure 3.2 Schematic illustration of the evolving structure of the IMSRG(2) Hamiltonian. The 𝑝 𝑝 𝑝′ reference state |Φ⟩ is completely decoupled from 1p-1h and 2p-2h excitations, |Φℎ ⟩ and |Φℎℎ′ ⟩ as a result of the flow. Figure adapted from Hergert et al. [4]. (CCSD), which aims to describe a comparable amount of many-body correlations [4]. 3.4 Generators for the IMSRG The IMSRG generator 𝜂 is an anti-Hermitian operator that we tailor in to implement a particular decoupling scheme and implicitly drive the Hamiltonian matrix toward its desired form. To this end, we need to identify the “off-diagonal” parts of the Hamiltonian that we would like to suppress through the flow. 3.4.1 Choice of Decoupling Scheme Here, we adopt the minimal decoupling scheme, in which we want to decouple a single eigenstate — usually the ground state — from the rest of the many-body Hamiltonian matrix [4]. Figure 3.2 schematically displays the structure of initial and evolved Hamiltonian in the Ap-Ah basis of excitations relative to a reference state |Φ⟩. In this scheme, the off-diagonal matrix elements of the evolving Hamiltonian that are targeted for suppression read [4] ⟨Φ|𝐻 (𝑠)|Φ𝑖𝑎 ⟩ = 𝑓𝑎𝑖 (𝑠) (3.19) ⟨Φ|𝐻 (𝑠)|Φ𝑖𝑎𝑏𝑗 ⟩ = Γ𝑎𝑏𝑖 𝑗 (𝑠), 32 and the normal-ordered off-diagonal Hamiltonian is constructed from these coefficients in the usual way [4], ∑︁  1 ∑︁ h i 𝑓𝑎𝑖 (𝑠)𝑁 𝑎 †𝑎 𝑎𝑖 + Γ𝑎𝑏𝑖 𝑗 (𝑠)𝑁 𝑎 †𝑎 𝑎 †𝑏 𝑎 𝑗 𝑎𝑖 + 𝐻.𝑐..  𝐻 od (𝑠) = (3.20) 𝑎𝑖 4 𝑎𝑏𝑖 𝑗 This definition must be generalized slightly in if we work with a general correlated reference state in the multi-reference IMSRG (MR-IMSRG, see Appendix B), or with a reference state ensemble (Chapter 6). 3.4.2 White Generator Using the decoupling scheme and definition of the off-diagonal Hamiltonian from the previous section, we can construct a class of generators that are inspired by the work of Steven White [9]. They are defined as [4], ∑︁ 𝑓𝑎𝑖 (𝑠)  †  1 ∑︁ Γ𝑎𝑏𝑖 𝑗 (𝑠) h † † i 𝜂W (𝑠) = 𝑁 𝑎 𝑎 𝑎𝑖 + 𝑁 𝑎 𝑎 𝑎 𝑏 𝑎 𝑗 𝑎𝑖 − 𝐻.𝑐.. (3.21) 𝑎𝑖 Δ𝑎𝑖 (𝑠) 4 𝑎𝑏𝑖 𝑗 Δ𝑎𝑏𝑖 𝑗 (𝑠) The anti-Hermiticity enters from the energy denominators Δ, which are inspired by the Epstein- Nesbet (EN) and Møller-Plesset (MP) versions of many-body perturbation theory. The EN denom- inators read (EN) Δ𝑎𝑖 = 𝑓𝑎 − 𝑓𝑖 − Γ𝑎𝑖𝑎𝑖 (EN) Δ𝑖𝑎 = 𝑓𝑖 − 𝑓𝑎 + Γ𝑎𝑖𝑎𝑖 (EN) Δ𝑎𝑏𝑖 𝑗 = 𝑓𝑎 + 𝑓 𝑏 − 𝑓𝑖 − 𝑓 𝑗 (3.22) + Γ𝑎𝑏𝑎𝑏 + Γ𝑖 𝑗𝑖 𝑗 − Γ𝑎𝑖𝑎𝑖 − Γ𝑏 𝑗 𝑏 𝑗 − Γ𝑎 𝑗 𝑎 𝑗 − Γ𝑏𝑖𝑏𝑖 Δ𝑖(EN) 𝑗 𝑎𝑏 = 𝑓𝑖 + 𝑓 𝑗 − 𝑓𝑎 − 𝑓 𝑏 − Γ𝑎𝑏𝑎𝑏 − Γ𝑖 𝑗𝑖 𝑗 + Γ𝑎𝑖𝑎𝑖 + Γ𝑏 𝑗 𝑏 𝑗 + Γ𝑎 𝑗 𝑎 𝑗 + Γ𝑏𝑖𝑏𝑖 , and the MP denominators are given by (MP) Δ𝑎𝑖 = 𝑓𝑎 − 𝑓𝑖 (MP) Δ𝑖𝑎 = 𝑓𝑖 − 𝑓𝑎 (3.23) (MP) Δ𝑎𝑏𝑖 𝑗 = 𝑓𝑎 + 𝑓𝑏 − 𝑓𝑖 − 𝑓 𝑗 Δ𝑖(MP) 𝑗 𝑎𝑏 = 𝑓𝑖 + 𝑓 𝑗 − 𝑓𝑎 − 𝑓 𝑏 . 33 For additional details, we refer to Hergert et al. [4]. An alternative formulation to the White generator regularizes the matrix elements via the   arctan function, so that “rotations” of the Hamiltonian are restricted to the interval − 𝜋4 , 𝜋4 [5]. This formulation proves useful in the reference ensemble (Chapter 6), where this additional constraint also constrains the optimization manifold of correlation parameters. The White-Atan generator reads ∑︁ 1 2 𝑓𝑎𝑖 (𝑠)  †  𝜂WA (𝑠) = arctan 𝑁 𝑎 𝑎 𝑎𝑖 𝑝ℎ 2 Δ 𝑎𝑖 (𝑠) (3.24) 1 ∑︁ 1 2Γ𝑎𝑏𝑖 𝑗 (𝑠) h † † i + arctan 𝑁 𝑎 𝑎 𝑎 𝑏 𝑎 𝑗 𝑎𝑖 − 𝐻.𝑐.. 4 𝑎𝑏𝑖 𝑗 2 Δ𝑎𝑏𝑖 𝑗 (𝑠) Finally, we note that the White generator can be implemented for a general, fractional reference state (e.g. the reference ensemble) by multiplying the relevant particle and hole indices with the corresponding particle and hole occupation numbers. This is the same strategy we use for density matrix normal-ordering in Chapter 2.8. For example, the matrix elements relevant to the one-body denominator Δ𝑎𝑖 in the White generator may be transformed to this generalized single-particle index scheme via, 𝑓𝑖 −→ 𝑓 𝑝 𝑛 𝑝 𝑓𝑎 −→ 𝑓 𝑝 𝑛 𝑝 (3.25) Γ𝑎𝑖𝑎𝑖 −→ Γ𝑝𝑞 𝑝𝑞 𝑛 𝑝 𝑛𝑞 𝑛 𝑝 𝑛𝑞 . By transforming the matrix elements in this way, we ensure the White generator is compatible with the MR-IMSRG flow equations. 3.4.3 Brillouin Generator In MR-IMSRG applications, the so-called Brillouin generator has proven to be a very useful and robust choice. It is defined by considering the change of the energy expectation value in the (potentially correlated) reference state under unitary variations. The evolved Hamiltonian is required to satisfy the so-called irreducible Brillouin conditions (IBCs) [6, 3], which are the 34 stationarity conditions for the unitary variation: ⟨Φ| 𝐻 (∞), 𝑁 𝑎 †𝑝 𝑎 𝑞 |Φ⟩ = 0    ⟨Φ| 𝐻 (∞), 𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 |Φ⟩ = 0    (3.26) . . . (including conjugates). Note that we cannot distinguish particle and hole indices here, since we allow correlated reference states and ensembles. In the limit where |Φ⟩ is a single Slater determinant, the IBCs reduce to the decoupling conditions in Equation 3.19. The Brillouin generator is written in second-quantized form as usual, ∑︁ 𝜂 𝑝𝑞 (𝑠)𝑁 𝑎 †𝑝 𝑎 𝑞   𝜂(𝑠) = 𝑝𝑞 (3.27) 1 ∑︁ 𝜂 𝑝𝑞𝑟 𝑠 (𝑠)𝑁 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 ,   + 4 𝑝𝑞𝑟 𝑠 where the matrix elements 𝜂 𝑝𝑞 (𝑠), 𝜂 𝑝𝑞𝑟 𝑠 (𝑠) are written in terms of the IBCs in Equation 3.26: 𝑎 †𝑝 𝑎 𝑞    𝜂 𝑝𝑞 (𝑠) ≡ ⟨Φ| 𝐻 (𝑠), 𝑁 |Φ⟩   † †  (3.28) 𝜂 𝑝𝑞𝑟 𝑠 (𝑠) ≡ ⟨Φ| 𝐻 (𝑠), 𝑁 𝑎 𝑝 𝑎 𝑞 𝑎 𝑠 𝑎𝑟 |Φ⟩ . Since these matrix elements vanish when the energy becomes stationary, we can think of the Brillouin generator as a kind of gradient of the energy under unitary variations. The commutators in Equation 3.28 are evaluated according to the generalized Wick contraction discussed in Chapter 2.6. The procedure for evaluating the commutators is similar to the derivation outlined in Chapter 3.3, but additional contractions appear that depend on the irreducible two- and three-body density matrices (see Chapter 2.6). Following Hergert [3], the Brillouin matrix elements are written in 35 Equations 3.29 and 3.30, 1 ∑︁  𝜂 𝑝𝑞 =(𝑛𝑞 − 𝑛 𝑝 ) 𝑓𝑞 𝑝 − Γ𝑞𝑟 𝑠𝑡 𝜆 𝑝𝑟 𝑠𝑡 − Γ𝑟 𝑠𝑝𝑡 𝜆𝑟 𝑠𝑞𝑡 (3.29) 2 𝑟 𝑠𝑡 𝜂 𝑝𝑞𝑟 𝑠 =Γ𝑟 𝑠𝑝𝑞 (𝑛 𝑝 𝑛𝑞 𝑛𝑟 𝑛 𝑠 − 𝑛 𝑝 𝑛𝑞 𝑛𝑟 𝑛 𝑠 ) ∑︁  + (1 − 𝑃 𝑝𝑞 ) 𝑓𝑡 𝑝 𝜆𝑡𝑞𝑟 𝑠 − (1 − 𝑃𝑟 𝑠 ) 𝑓𝑟𝑡 𝜆 𝑝𝑞𝑡𝑠 𝑡 1  + (𝜆Γ) 𝑟 𝑠𝑝𝑞 (1 − 𝑛 𝑝 − 𝑛𝑞 ) − (Γ𝜆) 𝑟 𝑠𝑝𝑞 (1 − 𝑛𝑟 − 𝑛 𝑠 ) (3.30) 2 ∑︁ + (1 − 𝑃 𝑝𝑞 − 𝑃𝑟 𝑠 ) (𝑛𝑞 − 𝑛𝑟 )Γ𝑡𝑟𝑢𝑞 𝜆𝑡 𝑝𝑢𝑠 𝑡𝑢 1 ∑︁  + (1 − 𝑃𝑟 𝑠 )Γ𝑟𝑡𝑢𝑣 𝜆𝑡 𝑝𝑞𝑢𝑣𝑠 − (1 − 𝑃 𝑝𝑠 )Γ𝑡𝑢 𝑝𝑣 𝜆𝑡𝑢𝑞𝑣𝑟 𝑠 . 2 𝑡𝑢𝑣 For additional details, we refer to Hergert [3]. For completely general reference states, the use of the Brillouin generator naively increases our storage requirements to 𝑂 (𝑁 6 ) because of its dependence on 𝜆 𝑝𝑞𝑟 𝑠𝑡𝑢 , and the computational cost increases to 𝑂 (𝑁 7 ) because of the summation in the last line of Equation 3.30. Fortunately, we can exploit that the irreducible three-body density matrix is either compact or highly structured for many relevant types of correlated reference states, so that the 𝑂 (𝑁 4 ) storage and 𝑂 (𝑁 6 ) computational cost of the standard IMSRG(2) approach can be preserved in MR-IMSRG(2) applications. 36 BIBLIOGRAPHY [1] S.K. Bogner, R.J. Furnstahl, and A. Schwenk. “From low-momentum interactions to nuclear structure”. In: Progress in Particle and Nuclear Physics 65.1 (2010), pp. 94– 147. issn: 0146-6410. doi: https://doi.org/10.1016/j.ppnp.2010.03.001. url: https: //www.sciencedirect.com/science/article/pii/S0146641010000347. [2] Stanisław D. Głazek and Kenneth G. Wilson. “Renormalization of Hamiltonians”. In: Phys. Rev. D 48 (12 Dec. 1993), pp. 5863–5872. doi: 10.1103/PhysRevD.48.5863. url: https://link.aps.org/doi/10.1103/PhysRevD.48.5863. [3] H Hergert. “In-medium similarity renormalization group for closed and open-shell nu- clei”. In: Physica Scripta 92.2 (Dec. 2016), p. 023002. issn: 1402-4896. doi: 10.1088/1402-4896/92/2/023002. url: http://dx.doi.org/10.1088/1402-4896/92/2/023002. [4] H. Hergert et al. “The In-Medium Similarity Renormalization Group: A novel ab initio method for nuclei”. In: Physics Reports 621 (Mar. 2016), pp. 165–222. issn: 0370- 1573. doi: 10.1016/j.physrep.2015.12.007. url: http://dx.doi.org/10.1016/j.physrep.2015. 12.007. [5] Morten Hjorth-Jensen, Maria-Paola Lombardo, and Ubirajara van Kolck. An advanced course in computational nuclear physics: bridging the scales from quarks to neutron stars. Springer, 2017. [6] Debashis Mukherjee and Werner Kutzelnigg. “Irreducible Brillouin conditions and con- tracted Schrödinger equations for n-electron systems. I. The equations satisfied by the density cumulants”. In: The Journal of Chemical Physics 114.5 (2001), pp. 2047– 2061. doi: 10.1063/1.1337058. eprint: https://doi.org/10.1063/1.1337058. url: https://doi.org/10.1063/1.1337058. [7] I. Shavitt and R.J. Bartlett. Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory. Cambridge Molecular Science. Cambridge University Press, 2009. isbn: 9780521818322. url: https://books.google.com/books?id=SWw6ac1NHZYC. [8] Franz J. Wegner. “Flow equations for Hamiltonians”. In: Nuclear Physics B - Proceedings Supplements 90 (2000). NON-PERTURBATIVE QCD AND HADRON PHENOMENOL- OGY, pp. 141–146. issn: 0920-5632. doi: https://doi.org/10.1016/S0920-5632(00) 00911-7. url: https://www.sciencedirect.com/science/article/pii/S0920563200009117. [9] Steven R. White. “Numerical canonical transformation approach to quantum many-body problems”. In: The Journal of Chemical Physics 117.16 (Oct. 2002), pp. 7472–7482. issn: 0021-9606. doi: 10.1063/1.1508370. eprint: https://pubs.aip.org/aip/jcp/article-pdf/117/ 16/7472/10843158/7472\_1\_online.pdf. url: https://doi.org/10.1063/1.1508370. 37 CHAPTER 4 MODEL HAMILTONIANS 4.1 The Pairing Model Plus Particle-Hole Model In this section, we introduce the model Hamiltonian which is used frequently in the body of work—namely, the pairing model. Information in this Chapter is compiled largely from Hjorth-Jensen et al. [5] and Hjorth-Jensen, Lombardo, and Kolck [6], although this model ap- pears elsewhere outside the context of second-quantized many-body physics. The pairing model is a convenient tool for testing computational and theoretical methods in nuclear physics. For one reason, the model is exactly solvable, because the many-body basis dimension is small enough that the exponential scaling is not an issue.1. We may probe the exact solution to the pairing model to validate the results for a particular method. For another reason, the pairing model captures qualitative features of more realistic nuclear interactions, despite its simplicity [5]. The pairing model consists of equidistant energy levels which contain up to two particles per level with opposite spin values, which we label + (up) and − (down) for convenience. The single- particle states are implicitly assumed to be the eigenstates of the one-particle operator, so that the one-body interaction term is diagonal. The two-body interaction term acts on pairs of particles with the same energy but opposite spin. Thus, the pairing model Hamiltonian can be written as ∑︁ 𝑔 ∑︁ † † 𝐻 𝑝𝑎𝑖𝑟𝑖𝑛𝑔 = 𝛿 ( 𝑝 − 1)𝑎 †𝑝𝜎 𝑎 𝑝𝜎 − 𝑎 𝑎 𝑎 𝑞− 𝑎 𝑞+ , (4.1) 𝑝𝜎 2 𝑝𝑞 𝑝+ 𝑝− where 𝜎 runs over all spin values, 𝛿 controls the energy level spacing, and 𝑔 controls the pairing strength. Note that the explicit − sign in the pairing term implies that we have an attractive pairing interaction for 𝑔 > 0, i.e., the formation of pairs is energetically favorable. Moreover, the interaction term cannot break pairs by moving the particles to different energy levels, and it does not change the spins of the pair it acts upon. This implies that both the number of pairs in the system and the 1 We will demonstrate how to extract the exact solution via Full Configuration Interaction in a later section. 38 p=4 p=3 𝜖 ↑ ↓ p=2 ↑ ↓ p=1 Figure 4.1 Pairing model of four particles at half filling, with two opposite spins per level. The vertical axis indicates increasing energy. The dashed line indicates the Fermi surface. total spin projection 𝐴 ℏ ∑︁ 𝑆𝑧 = 𝜎𝑧 (𝑖) (4.2) 2 𝑖=1 are conserved. To add complexity to our model, we can introduce an additional interaction term that can create one-particle-one-hole excitations and break up pairs. The Hamiltonian of this pairing-plus-particle- hole (P3H) model can be written as 𝑏 ∑︁ † † 𝐻 𝑃3𝐻 = 𝐻 𝑝𝑎𝑖𝑟𝑖𝑛𝑔 − (𝑎 𝑎 𝑎 𝑝 ′ − 𝑎 𝑞+ + 𝑎 †𝑞+ 𝑎 †𝑝 ′ − 𝑎 𝑝− 𝑎 𝑝+ ), (4.3) 2 𝑝 𝑝 ′ 𝑞 𝑝+ 𝑝− where 𝑏 controls the particle-hole / pair-breaking interaction strength. Note that the spin structure of this interaction was chosen such that 𝑆 𝑧 remains conserved. Figure 4.1 shows a ground-state configuration for an attractive pairing model with four particles distributed over four energy levels (for a total of eight single particle states). In the P3H model, on the other hand, we expect richer behavior because of the competition between the pairing and 39 b=0 b 0 Figure 4.2 The left panel features a P3H Hamiltonian with 𝑏 = 0, i.e. no pair-breaking. Black maps to zero. The right panel features a P3H Hamiltonian with 𝑏 ≠ 0. Grey maps to zero. The diagonal elements have been subtracted out to emphasize the off-diagonal structure of each Hamiltonian. pair-breaking terms, and the configuration shown in Figure 4.1 will most likely not be a good approximation for the ground state. In Figure 4.2, we show the qualitative structure of off-diagonal elements in the P3H Hamiltonian. We observe that the nonzero pair-breaking term probes new sectors of the Hamiltonian that were not accessible in the pure pairing interaction. In Figures 4.3 and 4.4, we visualize the eigenstate structure of several P3H Hamiltonians, with no pair-breaking and no pairing, respectively. The qualitative features of these eigenstates are similar, in nature, to wavefunctions in the nuclear shell model, hence the heavy focus on the P3H model in many of this results of this work. Each P3H Hamiltonian has been prepared in a basis of Slater determinants (SD) in an occupation number representation; therefore, each coefficient in the eigenvectors plotted in Figures 4.3 and 4.4 correspond to a particular SD configuration. For example, the first six coefficients indexed by v0 through v5 correspond to fully-paired configurations. Coefficients beyond v5 correspond to unpaired configurations, in the spin 𝑆 = 0 sub-block. By examining the wavefunction structure in the SD basis, we learn how important each SD 40 A 0.8 strong repulsive pairing 0.6 (g, b) = (-2.0000, 0.0000) 0.4 0.2 0.0 1.0 B 0.8 weak attractive pairing 0.6 (g, b) = (0.4828, 0.0000) 0.4 0.2 0.0 0.6 C strong attractive pairing 0.4 (g, b) = (2.0000, 0.0000) 0.2 0.0 v1 v3 v5 v7 v9 v11 v13 v15 v17 v19 v21 v23 v25 v27 v29 v31 v33 v35 Figure 4.3 Three pairing-only examples of wavefunction structure in the P3H model. The horizontal axis is the index of the eigenvector coefficient. The vertical axis is the corresponding squared eigenvector coefficient. configuration is in describing a particular eigenstate. For example, all three panels in Figure 4.3 immediately suggest that unpaired configurations contribute nothing to the eigenstates of the pairing-only P3H model. However, we see that the pairing strength controls how much each paired configuration contributes to the structure of the eigenstate. All three panels in Figure 4.4 show that most, or all, SD configurations are important to the eigenstate structure of the pair-breaking-only P3H model. Contributions for paired and unpaired configurations are especially important in the repulsive pair-breaking case. Understanding the structure of the P3H eigenstates in the SD configuration basis is useful for 41 D 0.3 strong attractive pair-breaking 0.2 (g, b) = (0.0000, -2.0000) 0.1 0.0 E 0.6 weak repulsive pair-breaking (g, b) = (0.0000, 0.4828) 0.4 0.2 0.0 F 0.08 strong repulsive pair-breaking 0.06 (g, b) = (0.0000, 2.0000) 0.04 0.02 0.00 v1 v3 v5 v7 v9 v11 v13 v15 v17 v19 v21 v23 v25 v27 v29 v31 v33 v35 Figure 4.4 Three pair-breaking-only examples of wavefunction structure in the P3H model. The horizontal axis is the index of the eigenvector coefficient. The vertical axis is corresponding squared eigenvector coefficient. the reference-state ensemble discussions in Chapter 6, because it offers insight on how to best to design an ensemble of SD configurations to capture the most important behavior in a particular eigenstate, or the full eigenspectrum. The eigenvectors in Figures 4.3 and 4.4 were obtained via FCI calculations of the P3H Hamil- tonians, which is the subject of Section 2. 4.2 Chiral Effective Field Theory for Nuclear Interactions In this section, we provide a brief overview of chiral effective field theory (𝜒EFT), focusing on breadth rather than depth, in order to build context for the results in Chapters 7 and 8. Machleidt 42 and Entem [7] present a detailed discussion on the derivation of nuclear interactions from 𝜒EFT, and Piarulli and Tews [8] provides a comprehensive review of local chiral nuclear interactions, in particular. We will skip any rigorous mathematical description of 𝜒EFT, and instead summarize where the low-energy constants (LECs), relevant to Chapters 7 and 8, appear in the theory. 𝜒EFT is a low-energy, effective formulation of quantum chromodynamics (QCD)—“effective”, in the sense that the most effective degrees of freedom are strategically chosen to systems of interacting nucleons at the energy scales relevant to nuclear science experiments. QCD is the fundamental theory of the strong interaction in the Standard Model of Particle Physics, and provides tools to describe interactions at the level of quarks and gluons, which are the constituents of nucleons and other strongly interacting composite particles. Since the description of even a single nucleon at this level is extremely challenging, the description of all but the lightest nuclei is all but infeasible. However, it is not necessary to do so because the quark substructure of individual nucleons is irrelevant for the description of nuclear structure and dynamics. Thus, 𝜒EFT pivots the description to the nucleons themselves, and uses pions as the mediators of the strong interaction at low energies. Essentially, we are lowering the resolution at which we describe the system, “blurring” out the substructure of nucleons and pions so that they can be treated as a point-like particles. The conceptual framework of any EFT posits that we need to consider all interactions that are consistent with the symmetries of the underlying theory. In addition to rotational, translational and Lorentz invariance that are required of fundamental interactions, it is necessary to account for the chiral symmetry of QCD [7, 8]. In principle, we can conceive of infinitely many interactions that respect those symmetries, but a properly constructed EFT provides a systematic framework for classifying their importance and truncating interactions that give negligible contributions to the observables of interest — in the case of 𝜒EFT, this is the power counting explained below. In the 𝜒EFT framework, interactions now fall into two broad categories — long-range interactions that are modeled by the exchange of pions (𝜋) between nucleons (𝑁), and short-range contact interactions between nucleons. The parameters associated with corrections to 𝜋𝑁 interaction vertices as well as the contact interactions parameterize the QCD physics that is not directly resolved by the EFT 43 in the most general, unbiased way, so as not to introduce modeling artifacts. Collectively, these parameters are referred to as the low-energy constants (LECs) of the theory. In principle, they could be determined from QCD via direct computations or matching results for observables that are accessible in both theories, but this is currently not feasible. Instead, the LECs are fit to low-energy QCD data, such as NN scattering phase shifts, or 3N observables. In addition to a 𝜒EFT with 𝑁 and 𝜋 degrees of freedom, we will also consider the so-called Δ-full 𝜒EFT, a variant that additionally includes the Δ isobar as a degree of freedom. This particle can be thought of as an excited 𝑁 𝜋 resonant state. The energy and momenta of this excitation are sufficiently low to overlap with the range we seek to describe with 𝜒EFT, so testing its impact on the description of nuclear observables is of high interest. We do note that the Δs are only included as virtual intermediate states — roughly speaking, this means that while Δs will not be observed as constituents of a nucleus on any experimental time scale, they can be formed for brief instants of time during an interaction process that might leave an imprint on the description of nuclear observables. In the Δ-less 𝜒EFT, the effects of the intermediate Δ excitations are implicitly included in LECs of the theory. The interactions relevant to this dissertation have been derived from Weinberg power-counting [9] of the ratio 𝑄/Λ𝑏 , where 𝑄 is a typical momentum scale and Λ𝑏 is the breakdown energy scale, or cutoff, beyond which 𝜒EFT is no longer a valid theory for the chosen degrees of freedom. Λ𝑏 is typically chosen in the range 500 − 1000 MeV, either based on the masses of particles that are not explicit degrees of freedom in the 𝜒EFT scheme, or combinations of parameters that are related to physics beyond the 𝜒EFT description. With typical momentum scales 𝑄 being of the order of the Fermi momentum 𝑘 𝐹 ≈ 1.4 fm−1 (corresponding to approximately 280MeV) or the pion mass 𝑚 𝜋 = 0.7 fm−1 (approx. 140 MeV), the ratio 𝑄/Λ𝑏 will then roughly be on the order of 1/5 to 1/3, so that higher-order processes are suppressed at a reasonable rate. The pion exchange and contact interactions can then be organized according to powers 𝑄/𝜆 𝜒 , leading to the organization scheme shown in Figure 4.5. In the common terminology, the (𝑄/Λ𝑏 ) 𝜈=0 = 1 interactions are referred to as leading order (LO), 𝜈 = 2 as next-to-leading or- 44 Figure 4.5 Figure taken from Piarulli and Tews [8]. On the left column are the NN force diagrams. On the right column are the 3N force diagrams. Solid lines are nucleons, while dashed lines are pions. Double solid lines are Δ-isobars. Each diagram may be identified with an LEC, which controls the strength of the contact force described by the diagram. der (NLO), 𝜈 = 3 as next-to-next-to-leading order (NNLO) and so on. Note the absence of 𝜈 = 1, because parity invariance stipulates that nuclear potentials cannot be linear in momentum [7, 8]. Diagrams with dashed lines involve the exchange of pions, while diagrams with only solid single lines are nuclear contact interactions. Double solid lines indicate Δ isobars in the Δ-full 𝜒EFT. The arrows indicate how Δ-full diagrams emerge from the short-range parameterizations in the Δ-less theory — note that these diagrams are usually promoted in order, which suggests that there should be differences in how rapidly observables converge in the two EFT variants, and which order of the expansion is required to achieve a desired precision. The relative momentum described by 𝑁 𝑁 scattering is decomposed into “partial wave” con- 45 tributions of differing angular momentum, where the S-wave partial waves corresponds to orbital angular momentum equal to 0, P-wave partial waves corresponds to orbital angular momentum equal to 1, D-wave partial waves corresponds to orbital angular momentum 2, and so on. The orbital angular momentum couples to spin to produce the partial waves, characterized by differing total angular momentum. At LO, only one-pion exchange and S-wave contact forces contribute. At NLO, we have two-pion exchanges plus subleading contact forces in the S wave and P waves. In the delta-full 𝜒EFT, a leading 3N force also appears at this order (also called the Fujita-Miywazawa force, named for the work of [4] in its derivation). At NNLO, there is no new contact term in the NN interaction, but higher-order corrections to the 𝜋𝑁 vertex appear: The LECs which control these contributions are 𝑐 1 , 𝑐 2 , 𝑐 3 , 𝑐 4 , where 𝑐 2 is only relevant in the Δ-full theory. The LECs 𝑐 1 , 𝑐 3 and 𝑐 4 also appear consistently in the long-range part of the NNLO 3N interaction, and in addition, a mixed long- and short-range term (proportional to 𝑐 𝐷 ) and a genuine contact (proportional to 𝑐 𝐸 ) contribute. At N3 LO contact forces up to D-waves appear in the NN force, adding to the number of parameters, in addition to more complicated but parameter-free 3N forces as well as a leading 4N force. In Tables 4.1 and 4.2 we list the “baseline” LEC sets for two different NNLO interactions (up to 3B forces) and one N3LO interaction (up to 2B forces), respectively, that will be explored in later chapters of this work. 46 ΔNNLOGO (394) NNLOsat (450) ( 𝑝 𝑝) 𝐶˜1𝑆0 -0.338142 -0.158149 (𝑛𝑝) 𝐶˜1𝑆0 -0.339250 -0.159822 (𝑛𝑛) 𝐶˜1𝑆 0 -0.338746 -0.159150 ( 𝑝 𝑝) 𝐶˜3𝑆1 -0.259839 -0.177674 (𝑛𝑝) 𝐶˜3𝑆1 -0.259839 -0.177674 (𝑛𝑛) 𝐶˜3𝑆 1 -0.259839 -0.177674 𝐶1𝑆0 +2.505389 +2.539368 𝐶3𝑃0 +0.700499 +1.398366 𝐶1𝑃1 -0.387960 +0.555959 𝐶3𝑃1 -0.964856 -1.136095 𝐶3𝑆1 +1.002189 +1.002893 𝐶3𝑆1 −3𝐷 1 +0.452523 +0.600716 𝐶3𝑃2 -0.883122 -0.802300 𝑐1 -0.74 -1.121521 𝑐2 -0.49 – 𝑐3 -0.65 -3.925006 𝑐4 +0.96 +3.765687 𝑐𝐷 +0.081 +0.081 𝑐𝐸 -0.002 -0.002 Table 4.1 ΔNNLOGO (394) [1] and NNLOsat (450) [2] interaction LEC values (including 𝑐 𝐸 and 𝑐 𝐷 relevant to 3B force). We refer to these as “baseline” LEC sets. We thank Andreas Ekström for providing 𝑁 𝑁 interaction routines for the relevant LEC values. 47 NN3 LO(500) ( 𝑝 𝑝) 𝐶˜1𝑆0 -0.145286 (𝑛𝑝) 𝐶˜ 1𝑆0 -0.147167 ˜ (𝑛𝑛) 𝐶1𝑆 0 -0.146285 ( 𝑝 𝑝) 𝐶˜3𝑆1 -0.118972 (𝑛𝑝) 𝐶˜3𝑆1 -0.118972 (𝑛𝑛) 𝐶˜3𝑆 1 -0.118972 𝐶1𝑆0 +2.38 𝐶3𝑃0 +1.487 𝐶1𝑃1 +0.656 𝐶3𝑃1 -0.63 𝐶3𝑆1 +0.76 𝐶3𝑆1 −3𝐷 1 +0.826 𝐶3𝑃2 -0.538 𝐷ˆ 1𝑆0 -2.545 𝐷 1𝑆0 -16.0 𝐷 3𝑃0 +0.245 𝐷 1𝑃1 +5.25 𝐷 3𝑃1 +2.35 𝐷ˆ 3𝑆1 +7.0 𝐷 3𝑆1 +6.55 𝐷 3𝐷 1 -2.8 ˆ 𝐷 3𝑆1 −3𝐷 1 +2.25 𝐷 3𝑆1 −3𝐷 1 +6.61 𝐷 1𝐷 2 -1.77 𝐷 3𝐷 2 -1.46 𝐷 3𝑃2 +2.295 𝐷 3𝑃2 −3𝐹2 -0.465 𝐷 3𝐷 3 +5.66 𝑐1 -0.81 𝑐2 +2.8 𝑐3 -3.2 𝑐4 +5.4 𝑑1 + 𝑑2 +3.06 𝑑3 -3.27 𝑑5 +0.45 𝑑14 − 𝑑15 -5.65 Table 4.2 “Baseline” LEC set for the N3 LO interaction with cutoff Λ = 500 MeV by Entem and Machleidt [3, 7]. Does not include LECs relevant to 3B force. 48 BIBLIOGRAPHY [1] A. Ekström et al. “Δ isobars and nuclear saturation”. In: Phys. Rev. C 97 (2 Feb. 2018), p. 024332. doi: 10.1103/PhysRevC.97.024332. url: https://link.aps.org/doi/10.1103/ PhysRevC.97.024332. [2] A. Ekström et al. “Accurate nuclear radii and binding energies from a chiral interaction”. In: Phys. Rev. C 91 (5 May 2015), p. 051301. doi: 10.1103/PhysRevC.91.051301. url: https://link.aps.org/doi/10.1103/PhysRevC.91.051301. [3] D. R. Entem and R. Machleidt. “Accurate charge-dependent nucleon-nucleon potential at fourth order of chiral perturbation theory”. In: Physical Review C 68.4 (Oct. 2003). doi: 10.1103/physrevc.68.041001. url: https://doi.org/10.1103%5C%2Fphysrevc.68.041001. [4] Jun-ichi Fujita and Hironari Miyazawa. “Pion Theory of Three-Body Forces”. In: Progress of Theoretical Physics 17.3 (Mar. 1957), pp. 360–365. issn: 0033-068X. doi: 10.1143/PTP. 17.360. eprint: https://academic.oup.com/ptp/article-pdf/17/3/360/5252121/17-3-360.pdf. url: https://doi.org/10.1143/PTP.17.360. [5] M Hjorth-Jensen et al. “Many-body interactions and nuclear structure”. In: Jour- nal of Physics G: Nuclear and Particle Physics 37.6 (Apr. 2010), p. 064035. doi: 10.1088/0954-3899/37/6/064035. url: https://doi.org/10.1088%2F0954-3899%2F37% 2F6%2F064035. [6] Morten Hjorth-Jensen, Maria-Paola Lombardo, and Ubirajara van Kolck. An advanced course in computational nuclear physics: bridging the scales from quarks to neutron stars. Springer, 2017. [7] R. Machleidt and D.R. Entem. “Chiral effective field theory and nuclear forces”. In: Physics Reports 503.1 (June 2011), pp. 1–75. doi: 10.1016/j.physrep.2011.02.001. url: https://doi.org/10.1016%5C%2Fj.physrep.2011.02.001. [8] Maria Piarulli and Ingo Tews. “Local Nucleon-Nucleon and Three-Nucleon Interactions Within Chiral Effective Field Theory”. In: Frontiers in Physics 7 (2020). issn: 2296- 424X. doi: 10.3389/fphy.2019.00245. url: https://www.frontiersin.org/articles/10.3389/ fphy.2019.00245. [9] Steven Weinberg. “High-Energy Behavior in Quantum Field Theory”. In: Phys. Rev. 118 (3 May 1960), pp. 838–849. doi: 10.1103/PhysRev.118.838. url: https://link.aps.org/doi/ 10.1103/PhysRev.118.838. 49 CHAPTER 5 BUILDING AN IMSRG MANY-BODY CODE USING TENSOR NETWORKS 5.1 Introduction Tensor network theory, in the context of this dissertation, describes the abstraction of multidi- mensional arrays and mathematical operations on them into a diagrammatic “formalism”, which simplifies the calculations which involve these often complex data structures. The diagrammatic method of tensor networks is somewhat less rigorous, but analogous, to the Goldstone and Hugen- holtz diagrammatic notation used in second-quantized many-body theory (see Chapter 4 of Shavitt and Bartlett [11]) to render many-body calculations less cumbersome. As a result, translating operators and working equations in methods like the IMSRG into tensor networks feels natural. In addition, tensor networks and tensor decompositions have been successfully applied to many-body methods in areas adjacent to this work, including Many-Body Perturbation Theory [14, 13], vari- ational Monte-Carlo [1], Tree Tensor Network States [6], Density Matrix Renormalization Group [12], and Tensor Network Renormalization [2]. In this body of work, we use tensor networks to provide an interface for abstraction of the computation of the IMSRG flow equations (e.g. Chapter 3, Equation 3.18). We want to offload the effort of computing the the IMSRG flow equations to a generalized tensor network library, in order to leverage the efficiency and optimization techniques implemented therein. With the computational load abstracted, we may freely exchange tensor network libraries, as well as sync our many-body codes with future improvements in those libraries. Ultimately, we strive for a separation of the workload between many-body physicists and computer scientists, so that the scientific progress benefits from ongoing optimizations of algorithms, without the friction associated with ad-hoc improvements to the many-body codebase. This strategy offers more long-term sustainability for the code relevant to the science, because of changing compiler and hardware landscapes. In addition, we present two IMSRG codes which were written or contributed to aid the studies that encompass this dissertation. These codes employ key ideas from tensor network theory, and 50 are discussed in Chapter 5.5 and Chapter 5.6. The purpose of these sections is to explain how these codes operate and the logic in their design. Appendix C features a catalog of these two codes, as well as their dependencies. Both IMSRG codes were designed with an object-oriented approach, so that the pieces which comprise the IMSRG “engine” are modular in nature. That is, users may e.g. write their own Hamiltonian class, inheriting from the appropriate interfaces, and plug directly into the rest of the IMSRG machine with relative simplicity. Additionally, one may write their own backends for performing the tensor contractions with solve the IMSRG flow equation terms. In principle, the administrative overhead in executing an object-oriented code is negligible compared to a sequential version of the same code; however, there are certainly areas within the codes we introduce that could benefit from a closer look at optimal memory management and code organization. For now, our goal is to provide simple access to validating results via the P3H model, as well as a template for improving the IMSRG production code with tensor libraries. Chapter 5.4 provides a detailed overview of how tensor networks are integrated into the com- putation of the IMSRG flow equations. These tensor networks are translated into code by tensor network calculation libraries such as TensorNetwork1 [9], built in Python, and the Tensor Algebra Compiler2 [5], built in C++. 5.2 Tensor Network Diagrams We may introduce a tensor, graphically as in Figure 5.1, as a node with a number of dangling edges corresponding to indices of a multidimensional array [9, 7]. Abstractly, the order of the edges does not matter (although they might matter in certain physics applications). As a general rule, dangling edges represent fixed indices, and edges between vertices (or loops) are summed over in an operation called a contraction. Contractions between tensors are performed graphically via connecting two edges. For example, a trace over one index in Figure 5.2 encodes the operation written mathematically, 1 Source: https://github.com/google/TensorNetwork 2 Source: https://github.com/tensor-compiler/taco 51 i j 𝐴 k l Figure 5.1 A tensor 𝐴 with some number of dangling edges. Evaluation of indices 𝑖 . . . 𝑗 . . . 𝑘 . . . 𝑙 gives the value of 𝐴 at that location. 𝑖 𝐴 𝐴 𝑖 𝑗 Figure 5.2 The trace operation over two dangling edges of the tensor 𝐴. ∑︁ tr( 𝐴) = 𝐴𝑖𝑖 . (5.1) 𝑖 Additional examples of common linear algebra operations, translated to a tensor diagram, are displayed in Figure 5.3. Notice that the resultant tensor is immediately identified by the dangling edges of the network, which remain fixed through the operation. A tensor network, however arbitrarily complex, will always involve dangling edges and/or contraction edges. Based on how the nodes and edges are connected in the diagram, we can always translate the diagram to mathematical form. 5.3 Factorization with Tensor Train Decomposition In general, the problems involving tensors require computing resources that scale exponentially in the dimension of the problem. Naive tensor storage (that is, a contiguous multidimensional array) is impractical or impossible to implement in large dimensional problems. As a result, 52 Í 𝑅= 𝑖 𝐴𝑖 𝐵𝑖 𝐴 𝑖 𝐵 Í 𝑖 𝑗 𝑅𝑖 = 𝑗 𝐴𝑖 𝑗 𝐵 𝑗 𝐴 𝐵 Í 𝑖 𝑘 𝑗 𝑅𝑖 𝑗 = 𝑘 𝐴𝑖𝑘 𝐵 𝑘 𝑗 𝐴 𝐵 Figure 5.3 Depicted are the vector-vector multiplication, matrix vector multiplication, and matrix- matrix multiplication in tensor diagrams. we seek efficient and effective representations of tensors that, in general, do not require storage of the complete tensor—instead, the tensor is decomposed into characteristic factors which can be contracted in order to reconstruct the full tensor. The tensor train decomposition (or TT- decomposition) is one such factorization technique introduced by Oseledets [8] as an alternative to the NP-hard problem of the canonical decomposition [8], 𝑟 ∑︁ 𝐴(𝑖1 , 𝑖2 , . . . , 𝑖 𝑑 ) = 𝑈1 (𝑖1 , 𝛼)𝑈2 (𝑖 2 , 𝛼) . . . 𝑈𝑑 (𝑖 𝑑 , 𝛼), (5.2) 𝛼=1 where 𝐴(𝑖1 , 𝑖2 , . . . , 𝑖 𝑑 ) represents the element of 𝐴 associated with the indices 𝑖1 , 𝑖2 , . . . , 𝑖 𝑑 . Note that 𝑖 𝑘 spans an unspecified range 𝑛 𝑘 , which may or may not be equal to another index range 𝑛𝑙 . For most methods involving the canonical decomposition, the decomposition rank 𝑟 must be known and the algorithms for computing the low-rank tensor components are not stable [8]. Figure 5.4 illustrates a schematic diagram of the canonical decomposition. There are 𝑑 canonical factors which correspond to the tensor dimension 𝑑. Each decomposition rank features an outer product over all canonical factors in that rank, and the sum over all ranks generate the full tensor 𝐴. The decomposition rank 𝑟 is taken as the minimum number of terms which reproduce 𝐴 satisfying an error threshold 𝜖. The TT-decomposition seeks a tensor 𝐴 which approximates a target tensor 𝐵 such that 𝐴 ≈ 𝐵. 53 𝑖1 𝑖𝑑 𝑖1 . . . 𝑖 𝑑 Í𝑟 = ( 𝛼) ( 𝛼) Ë Ë 𝐴 𝑢1 ... 𝑢𝑑 𝛼=1 Figure 5.4 Schematic tensor network diagram of the canonical decomposition. The rank 1 tensors (vectors) are referred to as canonical factors. The elements of 𝐴, which is characterized by 𝑑 indices, are computed via 𝐴(𝑖 1 , 𝑖2 , . . . , 𝑖 𝑑 ) = 𝐺 (𝑖1 )𝐺 (𝑖2 ) . . . 𝐺 (𝑖 𝑑 ), (5.3) where 𝐺 (𝑖 𝑘 ) is a matrix with dimension 𝑟 𝑘−1 × 𝑟 𝑘 . The elements of the tensor 𝐴 are expressed as a matrix product. Thus, instead of storing the number of elements corresponding to 𝐴, which might Í be as large as 𝑛 𝑑 , we store a number of elements equal to 𝑑𝑘=1 𝑟 𝑘−1𝑟 𝑘 which is considerably more efficient in most cases. Figure 5.5 draws a TT-decomposition diagram. The sum over decomposition ranks does not appear in the TT-decomposition, in contrast to the canonical decomposition in Figure 5.4. The auxiliary indices 𝛼 which contract between TT cores are fixed. The cores which reproduce the elements of 𝐴, which we denote with 𝐺, are structurally different than the cores which reproduce the full tensor A (in bold-face), which we denote in bold-face G. However, we emphasize that G is a full realization of the parameters that characterize 𝐺. The procedure for generating a TT-decomposition is not unique in general, due to algorithms which minimize a norm deviation, such as the Frobenius error. In these cases minimization may not yield unique results. In this study, we have not used an algorithmic procedure to generate the relevant TT-decomposition, instead relying on ad hoc inspection, because the decomposition are simple enough to derive by hand. The inspection method is best explained by example. Particularly helpful explanations can be found in [10]. Consider a rank 2 tensor 𝐴 with elements given by 𝐴(𝑎, 𝑏) = 𝑛𝑎 − 𝑛 𝑏 , where the right-hand side defines a particular element, corresponding to (𝑎, 𝑏). 54 𝑖1 𝑖2 𝑖𝑑 𝑖1 . . . 𝑖 𝑑 A = G1 𝛼1 G2 𝛼2 ... 𝛼𝑑−1 G𝑑 Figure 5.5 Schematic tensor network diagram of the tensor-train decomposition. When drawing the reconstruction of the full tensor, note that the boundary cores appear as rank 2 tensors, while inner cores appear as rank 3 tensors. The elements 𝑛𝑎 and 𝑛 𝑏 are numbers which come from two distinct (or equivalent) sets. Let the size of these sets be equal to 𝑁. Then, the total number of elements stored by the tensor 𝐴(𝑎, 𝑏) is 𝑁 2. We can generate a TT-decomposition by inspection from the analytical representation of the tensor elements 𝐴(𝑎, 𝑏). We need a number of smaller tensors 𝐺 equal to 2, by Equation 5.3 (𝑑 = 2). For example,   h i 1  𝐴(𝑎, 𝑏) = 𝐺 (𝑛𝑎 )𝐺 (𝑛 𝑏 ) = 𝑛𝑎  = 𝑛𝑎 − 𝑛 𝑏 .   1  (5.4) −𝑛   𝑏   In order to evaluate the tensor element at a particular index location, we need only evaluate the multiplication of the smaller tensor “cores” corresponding to that index location. Note that under the TT-decomposition scheme, we need only store the decomposed tensors, which in this case amounts to 2𝑁 + 2𝑁 = 4𝑁 total elements instead of 𝑁 2 elements to store the full tensor—a reduction in memory scaling from quadratic to linear (refer to Figure 5.6 for an example of how memory scales with dimension). The disadvantage, of course, is the computational cost associated with evaluating each tensor element dynamically. Later we will show that for a sufficient computing environment, this cost is negligible for the TT-decomposition relevant to this body of work. Although, we emphasize that the cost is dependent on the problem for which the TT-decomposition has been applied. 55 Figure 5.6 Memory scaling for the rank 2 tensor in Equation 5.4, comparing storage of the full matrix to the storage for the tensor train components. Notice that the TT-decomposition in Equation 5.4 can be evaluated to compute the full tensor A by building two matrices G(𝑛𝑎 ),G(𝑛 𝑏 ) from all possible 𝑛𝑎 ,𝑛 𝑏 in tensor cores 𝐺 (𝑛𝑎 ),𝐺 (𝑛 𝑏 ). Let 𝑁 = 4, and let the possible 𝑛𝑎 = 𝑛 𝑏 = [1, 1, 0, 0]] for 𝑎 = 1, 2, 3, 4, respectively. Then,   1 1       1 1   1 1 1 1 G(𝑛𝑎 )G(𝑛 𝑏 ) =      0 1 −1 −1 0 0        0 1    (5.5)   0 0 1 1    0 0 1 1 =  = A.  −1 −1 0 0      −1 −1 0 0   In this way, the tensor cores can be fully contracted to recover the full tensor 𝐴. An important feature of this exact decomposition is that we can used the decomposed tensor cores to perform calculations without ever explicitly constructing the full tensor 𝐴. Consider the matrix-vector 56 operation Av, where vT = [1, 2, 3, 4],       0 0 1 1 1  7             0 0 1 1 2  7  Av =   =        (5.6) −1 −1 0 0 3 −3                  −1 −1 0 0 4 −3           1 1 1          1 1  1 1 1 1 2 G(𝑛𝑎 )G(𝑛 𝑏 )v =       0 1 −1 −1 0 0 3           0 1 4      (5.7)     1 1  7          1 1  10  7 = =  .        0 1 −3 −3               0 1  −3     This short example is not meant to formally prove this fact, but to provide intuition for why this method works. The outer axes of the decomposition contain information about the elements of the full tensor. These axes directly interact with the multiplying factor, and the auxiliary axes contract only among themselves. Therefore, as long as the full tensor cores are stored for a particular (exact) TT-decomposition, calculations involving the full tensor can be performed by substituting the TT- decomposition. In this idealized example of exact decomposition, the operation in Equation 5.6 scales 𝑂 (𝑁 2 ), while the operation in Equation 5.7 scales 𝑂 (2 ∗ 𝑁 + 𝑁 ∗ 2) = 𝑂 (4𝑁). In this special case, the complexity is actually reduced (not accounting for the associated cost of performing the decomposition, which in this special case was performed ad hoc by inspection). An additional benefit of the TT-decomposition is efficiency in large-scale calculations. Our hope is that the cost of evaluating the TT-decomposition, in tandem with e.g. a matrix calculation, does not affect performance when compared to standard numerical methods; in Figure 5.7 we show that factorization via TT-decomposition of the two input matrices to a random matrix-matrix calculation achieves comparable time scaling to standard, optimized techniques. Note that the 57 Figure 5.7 Time scaling for performing matrix-matrix multiplication, while comparing several methods. Each curve computes multiplication of two random, dense, square matrices of varying dimension. The calculations were performed in Python using TensorNetwork [9] and numpy. TT-decomposed calculations seem to scale slightly better than the naive 𝑂 (𝑁 3 ), which is likely due to a similar scaling reduction to the example in Equations 5.6 and 5.7. 5.4 Tensor Networks in the IMSRG As discussed in the Introduction (Chapter 5.1), one of the reasons for viewing the IMSRG through a tensor network lens is to sequester the effort to compute the tensor contractions in Equation 3.18, the IMSRG(2) flow equations, into a separate library, or modular environment, where the optimization can be controlled independently of the many-body code. Additionally, we leverage the TT-decomposition discussed in Chapter 5.3 to bundle the decomposed occupation factors in Equation 3.18 with the IMSRG coefficient tensor contraction operations. In the future, we want to leverage tensor factorization beyond the occupation factors, and apply these factorization ideas to the IMSRG flow terms (Equation 3.18) directly. Zhu, Wirth, and Hergert [15] has shown success in factorizing the 𝜒EFT interactions (see Chapter 4.2) used as inputs for the IMSRG, which motivates further effort in a factorized IMSRG flow. Based on the discussion in Chapter 5.2, we observe that any term in Equation 3.18 which 58 𝑎 𝑏 𝜂 ′ (𝑠) 𝑐 Γ(𝑠) 𝑑 Figure 5.8 Tensor diagram corresponding to the two-body contribution to the reference energy in the IMSRG(2) flow equations in Equation 3.18. contains a multiplication of two tensors, as well as a sum over a shared index between those two tensors, can be translated to a tensor network diagram. For example, consider the two-body contribution to the reference energy (ignoring numerical factors):   (2) 𝑑𝐸 (𝑠) ∑︁ = 𝜂𝑎𝑏𝑐𝑑 (𝑠)Γ𝑐𝑑𝑎𝑏 (𝑠)𝑛𝑎 𝑛 𝑏 𝑛𝑐 𝑛 𝑑 . (5.8) 𝑑𝑠 𝑎𝑏𝑐𝑑 While this term actually contains three tensors to contract, we can leverage the fact that the occupation factor in 𝑛𝑎 𝑛 𝑏 𝑛𝑐 𝑛 𝑑 is just a number independent of 𝑠; we may freely commute the occupation factor among the IMSRG(2) coefficients, and combine with the 𝜂 factor via elementwise multiplication to obtain   (2) 𝑑𝐸 (𝑠) ∑︁ = 𝜂′𝑎𝑏𝑐𝑑 (𝑠)Γ𝑐𝑑𝑎𝑏 (𝑠), (5.9) 𝑑𝑠 𝑎𝑏𝑐𝑑 where 𝜂′𝑎𝑏𝑐𝑑 (𝑠) is now 𝜂′𝑎𝑏𝑐𝑑 (𝑠) = 𝑛𝑎 𝑛 𝑏 𝑛𝑐 𝑛 𝑑 · 𝜂𝑎𝑏𝑐𝑑 (𝑠). (5.10) The tensor network diagram corresponding to this term might look like Figure 5.8. Let us express the elements of the left hand tensor 𝜂′ (𝑠) in a network diagram, assuming that the occupation factor 𝑛𝑎 𝑛 𝑏 𝑛𝑐 𝑛 𝑑 has been TT-decomposed (see Appendix A.1, Equation A.5). This particular occupation factor is especially easy to decompose, because the analytical form of the tensor elements is a simple product. Then, the decomposition network would be evaluated for each element of the full tensor 𝜂′ (𝑠). Figure 5.9 illustrates where the TT-decomposition of the occupation factor integrates with the network for this term, which further increases the efficiency of evaluating the full term in a tensor network architecture. 59 𝐺 𝑎 (𝑛1 ) 𝐺 𝑏 (𝑛1 ) 𝐺 𝑐 (𝑛1 ) 𝐺 𝑑 (𝑛1 ) 𝑎 ... 𝑏 𝐺 𝑎 (𝑛 𝑁 ) 𝐺 𝑏 (𝑛 𝑁 ) 𝐺 𝑐 (𝑛 𝑁 ) 𝐺 𝑑 (𝑛 𝑁 ) Γ(𝑠) ⊙ 𝑐 𝜂 𝑑 Figure 5.9 Diagram illustrating where the tensor train decomposition of the occupation factor 𝑛𝑎 𝑛 𝑏 𝑛𝑐 𝑛 𝑑 integrates with the tensor network to compute Equation 5.9. The elements of the occupation factor are generated via repeated evaluations of the corresponding TT-decomposition. Then, each element is multiplied elementwise with the tensor 𝜂(𝑠) to compute 𝜂′ (𝑠). Finally, 𝜂′ (𝑠) is fully contracted with Γ(𝑠) to complete the diagram. As another example, let us consider the not-permuted piece of the two-body contribution to the one-body flowing operator 𝑓 (𝑠), written (ignoring numerical factors)   (2) 𝑑𝑓 𝑝𝑞 (𝑠) ∑︁ = (𝑛𝑎 𝑛 𝑏 𝑛𝑐 + 𝑛𝑎 𝑛 𝑏 𝑛𝑐 )𝜂𝑐 𝑝𝑎𝑏 (𝑠)Γ𝑎𝑏𝑐𝑞 (𝑠) 𝑑𝑠 𝑎𝑏𝑐 ∑︁ (5.11) = 𝜂′𝑐 𝑝𝑎𝑏 (𝑠)Γ𝑎𝑏𝑐𝑞 (𝑠), 𝑎𝑏𝑐 where again we have absorbed the 𝑠-independent occupation factor into the factor corresponding to 𝜂′ using elementwise multiplication. In this case, we wish to elementwise-multiply two tensors of different rank. This process amounts to multiplying the elements of the occupation factor 𝐶 (1) (𝑎, 𝑏, 𝑐) (refer to Appendix A) in the relevant index range of 𝜂𝑐 𝑝𝑎𝑏 , being careful to ensure that the information stored along axis 𝑝 remains unchanged. We can broadcast 𝐶 (1) (𝑎, 𝑏, 𝑐) to a rank 4 tensor, with only identity information stored along the last axis corresponding to 𝑝, via a tensor product with the identity operator, 𝐶 (1) (𝑎, 𝑏, 𝑐, 𝑝) = 𝐶 (1) (𝑎, 𝑏, 𝑐) ⊗ 1l. (5.12) 60 𝑎 𝑝 𝑞 𝜂 ′ (𝑠) 𝑏 Γ(𝑠) 𝑐 Figure 5.10 Tensor diagram corresponding to the two-body contribution to the one-body flowing operator in the IMSRG(2) flow equations in Equation 3.18. Finally, we permute the indices of 𝐶 (1) (𝑎, 𝑏, 𝑐, 𝑝) to align with 𝜂𝑐 𝑝𝑎𝑏 and then the elementwise multiplication can be performed to arrive at 𝜂′𝑐 𝑝𝑎𝑏 . The tensor network corresponding to Equation 5.11 is captured in Figure 5.10. The dangling edges immediately show the result is a rank 2 tensor which is consistent with the required shape of the one-body operator flow equation. 5.5 Tensor-factors IMSRG: Python Implementation for the P3H Model The code called Tensor-factors IMSRG, or tfimsrg, is a Python implementation of the IMSRG(2) (with developing support for IMSRG(3)) for the P3H model (see Chapter 4.1 and also Appendix C.1). The code was originally developed to demonstrate the feasibility of using Google’s TensorFlow tensor contraction API to express the IMSRG flow equations, interpreting the complicated sums as contractions between tensors. However, the project now uses the ncon API of the tensor contraction library TensorNetwork [9], which evaluates the contractions using the numpy backend. TensorNetwork provides several additional backends including JAX, TensorFlow, and PyTorch. The codebase was designed as an importable Python package rather than a standalone exe- cutable. Its overall structure (excluding extraneous files and directories) is: -- tfimsrg/ |-- main.py -- oop_imsrg/ |-- flow.py |-- generator.py |-- hamiltonian.py |-- occupation_tensors.py The primary pieces of the IMSRG flow machine are organized into objects within the 61 oop_imsrg/ directory. Flow.py computes the IMSRG flow equations in Equation 3.18 (as well as an implementation of the multi-reference IMSRG, see Appendix B). Generator.py computes the generators from Chapter 3.4. Hamiltonian.py computes the model Hamiltonian, the P3H model, from Chapter 4.1. Occupation_tensors.py pre-computes the TT-decomposed occupa- tion factors from Appendix A. These components are bundled and passed to the scipy solver, which integrates the solution of the flow equations, in main.py. 5.5.1 Performance Study The performance of the tfimsrg code, evaluating the P3H model, was measured relative to the naive implementation of IMSRG(2) in Python using matrix calculations introduced in Chapter 10 of Hjorth-Jensen, Lombardo, and Kolck [4]. Figure 5.11 shows the total wall time for several implementations of the IMSRG(2) code that (approximately) solve the P3H model with zero pair-breaking, including the matrix implementation, occupation factors constructed with loops, occupation factors constructed with contractions of TT-decomposed tensors (see Chapter 5.3 and Appendix A), and a TT-decomposition implementation exported to the GPU. Exporting the calculations to the GPU offers the most advantageous scaling in the context of these Python codes. All implementations that outsource tensor contractions in the IMSRG(2) flow equation to the TensorNetwork library perform faster than the naive implementations using matrix calculations. In Figure 5.12 we show that the tensor library implementation of the IMSRG(2) code captures the weak scaling feature of the single reference IMSRG, when particle and hole states are distinct. Note that the total number of single particle states 𝑁 = 𝑁 ℎ + 𝑁 𝑝 , where 𝑁 ℎ is the number of hole states and 𝑁 𝑝 is the number of particle states. For large bases in realistic calculations, we have 𝑁 𝑝 ≈ 𝑁. The performance results for this IMSRG(2) implementation are motivation to extend the tensor contraction strategy to a compiled language like C++, as well as modify the tensor contraction approach to take advantage of weak scaling. We want to also leverage this technology to the IMSRG production codes (Appendix C.1) which are written in C/C++. 62 Figure 5.11 Total wall time versus the number of single particle states in the P3H model, excluding the pair-breaking term. Although not depicted here, the relative scaling for a P3H model with nonzero pair-breaking is not appreciably different. In general, the IMSRG(2) scales 𝑂 (𝑁 6 ) with the number of single particle states 𝑁; however using the TensorNetwork library allows us to soften the prefactor on that scaling to reach larger basis sizes. IMSRG(2) execution time versus basis size (fixed num hole states) 8 Nh 2 wall-time microseconds (ln) 7 6 4 5 6 8 4 10 3 14 2 16 1 0 1.5 2.0 2.5 3.0 3.5 num S.P. states (ln) Figure 5.12 Total wall time versus the number of single particle states, for a fixed number of hole states. The P3H model parameters for input Hamiltonian are 𝛿 = 1.0, 𝑔 = 0.5, 𝑏 = 0.0. Although not depicted here, the relative scaling for a P3H model with nonzero pair-breaking is not appreciably different. 63 5.6 Tensorized C++ IMSRG: C++ Implementation for the P3H Model The code called Tensorized C++ IMSRG, or tcimsrg, is a C++ implementation of the IMSRG(2) for solving the P3H model (cf. Appendix C.1). This code is, more or less, a translation of tfimsrg from Python to C++. This code also includes more features for customizing the IMSRG(2) flow at execution time, including an option to use Dynamic Mode Decomposition (DMD) emulation of the IMSRG flow (see Chapter 7.4). This code was designed to run as an executable, with command line options similar to the IMSRG production code [3] (see also Appendix C.1). Again, we have taken an object-oriented approach to keep the code modular for future development and feature addition. The breakdown of objects is similar to tfimsrg, but we have tried to structure it according to best practices for C++: -- tcimsrg/ |-- CMakeLists.txt -- build/ |-- solve_imsrg -- include/ |-- BACKEND.hpp |-- DMD.h |-- flow_imsrg2.hpp |-- generator.hpp |-- pairinghamiltonian.hpp |-- occupation_factors.hpp |-- state_type.hpp |-- system.hpp |-- imsrg_utils.hpp -- src/ |-- BACKEND_taco.cpp |-- BACKEND_ublas.cpp |-- DMD.c |-- flow_imsrg2.cpp |-- white.cpp |-- white_atan.cpp |-- brillouin.cpp |-- pairinghamiltonian.cpp |-- occupation_factors.cpp |-- system.cpp |-- main.cpp The object headers are organized within the include/ directory, and the object source files 64 which implement the instructions in their respective header files are organized under the src/ directory. Running cmake within the build/ directory will build the executable solve_imsrg. BACKEND_taco.cpp and BACKEND_ublas.cpp set up the IMSRG(2) flow equation evaluation us- ing TACO and uBLAS loops, respectively. DMD.c sets up the DMD emulation. Flow_imsrg2.cpp sets up calculation the IMSRG(2) flow equation terms, based on the user’s specification of the back- end. White.cpp, white_atan.cpp, and brillouin.cpp set up the generators specified in Chap- ter 3.4. Pairinghamiltonian.cpp sets up the P3H model Hamiltonian according to Chapter 4.1. Occupation_tensors.cpp pre-computes the TT-decomposed occupation factors from Appendix A. System.cpp organizes the ODE system which is passed to the Boost library’s odeint tool. Finally, the file main.cpp combines all relevant objects for computing a full IMSRG(2) calculation into the executable solve_imsrg. 5.6.1 Performance Study Here we discuss the performance comparison between two different implemented calculation backends, naive loops using uBLAS vectors and contractions of tensors using the TACO library. The proceeding results include TT-decomposed occupation factors contracted directly in the network evaluation of the IMSRG(2) flow terms. Figure 5.13 shows the speedup of the TACO implementa- tion versus the naive loops implementation. Again, the tensor library contraction combined with TT-decomposed occupation factors softens the prefactor on the IMSRG(2) computational scaling, allowing us to reach larger basis sizes. The results in Figure 5.13 reflect implementations that are not multithreaded; however, we expect the results would be similar provided both implementations are optimally multithreaded. In Figure 5.14, we show the wall-time associated with increasing single particle basis size, varying the fixed number of hole states (details as in Figure 5.12). Here, the results are consistent even across the number of fixed hole states, in contrast to the Python implementation. This is likely due to differing tensor library implementations. Again, the weak scaling is captured. 65 Tloops/TTACO 10.0 speedup 7.5 5.0 2.5 0.0 4 6 8 10 12 14 16 18 numStates N Figure 5.13 Associated speedup of TACO tensor contractions to the naive loop calculations, versus the number of single particle states in the P3H model. The red dotted line marks a ratio equal to 1, where the performance of the two implementations is equal. IMSRG(2) execution time versus basis size (fixed num hole states) Nh 20 2 19 4 6 Wall-time (ln) 18 8 10 17 12 16 14 15 1.50 1.75 2.00 2.25 2.50 2.75 3.00 3.25 3.50 num S.P. states (ln) Figure 5.14 Details as in Figure 5.12, except as measured by the TACO C++ implementation. Note that wall-time is longer than the Python implementation due to differing integration strategies between the codes. 66 BIBLIOGRAPHY [1] Stephen R Clark. “Unifying neural-network quantum states and correlator product states via tensor networks”. In: Journal of Physics A: Mathematical and Theoretical 51.13 (Feb. 2018), p. 135301. issn: 1751-8121. doi: 10.1088/1751-8121/aaaaf2. url: http: //dx.doi.org/10.1088/1751-8121/aaaaf2. [2] G. Evenbly and G. Vidal. “Tensor Network Renormalization”. In: Physical Review Letters 115.18 (Oct. 2015). issn: 1079-7114. doi: 10.1103/physrevlett.115.180405. url: http: //dx.doi.org/10.1103/PhysRevLett.115.180405. [3] H. Hergert et al. “The In-Medium Similarity Renormalization Group: A novel ab initio method for nuclei”. In: Physics Reports 621 (Mar. 2016), pp. 165–222. issn: 0370- 1573. doi: 10.1016/j.physrep.2015.12.007. url: http://dx.doi.org/10.1016/j.physrep.2015. 12.007. [4] Morten Hjorth-Jensen, Maria-Paola Lombardo, and Ubirajara van Kolck. An advanced course in computational nuclear physics: bridging the scales from quarks to neutron stars. Springer, 2017. [5] Fredrik Kjolstad et al. “The Tensor Algebra Compiler”. In: Proc. ACM Program. Lang. 1.OOPSLA (Oct. 2017), 77:1–77:29. issn: 2475-1421. doi: 10.1145/3133901. url: http://doi.acm.org/10.1145/3133901. [6] V. Murg et al. “Tree Tensor Network State with Variable Tensor Order: An Efficient Mul- tireference Method for Strongly Correlated Systems”. In: Journal of Chemical Theory and Computation 11.3 (2015). PMID: 25844072, pp. 1027–1036. doi: 10.1021/ct501187j. eprint: https://doi.org/10.1021/ct501187j. url: https://doi.org/10.1021/ct501187j. [7] Román Orús. “A practical introduction to tensor networks: Matrix product states and projected entangled pair states”. In: Annals of Physics 349 (Oct. 2014), pp. 117–158. issn: 0003-4916. doi: 10.1016/j.aop.2014.06.013. url: http://dx.doi.org/10.1016/j.aop.2014.06. 013. [8] I. V. Oseledets. “Tensor-Train Decomposition”. In: SIAM Journal on Scientific Computing 33.5 (2011), pp. 2295–2317. doi: 10.1137/090752286. eprint: https://doi.org/10.1137/ 090752286. url: https://doi.org/10.1137/090752286. [9] Chase Roberts et al. TensorNetwork: A Library for Physics and Machine Learning. 2019. arXiv: 1905.01330 [physics.comp-ph]. [10] Anton Rodomanov. Introduction to the Tensor Train Decomposition and Its Applications in Machine Learning. HSE Seminar on Applied Linear Algebra. Mar. 2016. url: https: //bayesgroup.github.io/team/arodomanov/tt_hse16_slides.pdf. 67 [11] I. Shavitt and R.J. Bartlett. Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory. Cambridge Molecular Science. Cambridge University Press, 2009. isbn: 9780521818322. url: https://books.google.com/books?id=SWw6ac1NHZYC. [12] Szilárd Szalay et al. “Tensor product methods and entanglement optimization for ab initio quantum chemistry”. In: International Journal of Quantum Chemistry 115.19 (2015), pp. 1342–1391. doi: 10.1002/qua.24898. eprint: https://onlinelibrary.wiley.com/doi/pdf/ 10.1002/qua.24898. url: https://onlinelibrary.wiley.com/doi/abs/10.1002/qua.24898. [13] A. Tichai, J. Ripoche, and T. Duguet. “Pre-processing the nuclear many-body problem”. In: The European Physical Journal A 55.6 (June 2019). issn: 1434-601X. doi: 10.1140/epja/ i2019-12758-6. url: http://dx.doi.org/10.1140/epja/i2019-12758-6. [14] A. Tichai et al. “Tensor-decomposition techniques for ab initio nuclear structure calculations: From chiral nuclear potentials to ground-state energies”. In: Physical Review C 99.3 (Mar. 2019). issn: 2469-9993. doi: 10.1103/physrevc.99.034320. url: http://dx.doi.org/10. 1103/PhysRevC.99.034320. [15] B. Zhu, R. Wirth, and H. Hergert. “Singular value decomposition and similarity renormal- ization group evolution of nuclear interactions”. In: Physical Review C 104.4 (Oct. 2021). doi: 10.1103/physrevc.104.044002. url: https://doi.org/10.1103%5C%2Fphysrevc.104. 044002. 68 CHAPTER 6 OPTIMAL REFERENCE STATE ENSEMBLES FOR THE IN-MEDIUM SIMILARITY RENORMALIZATION GROUP 6.1 Introduction In this chapter, we focus on the main methodological improvement to the IMSRG framework that is developed in this thesis. In general, the implementation of any RG method begins with the choice of an operator basis. If we can work in a complete basis, we can capture the RG flow exactly (for example, SRG for finite-size matrices). In most applications, and in particular the IMSRG(2), we only work with truncated bases due to the complexity of tracking induced operators. For example, in the IMSRG(2) we work with a finite set of normal ordered, second-quantized operators. The RG flow comes with built in diagnostics—it can tell us which operators are relevant or irrelevant, and stabilities or divergences can be a sign of truncations that are too severe, or truncations that violate symmetries we would like or need to preserve. The response of the RG flow to the choice of the operator basis impacts the accuracy with which we can describe the observables of interest. Here, we will explore the use of ensembles of reference states to enhance the robustness of the IMSRG operator basis, and to control truncation errors. We show that the use of an ensemble serves to “inform” the IMSRG flow about the excited states of the system, so that the RG improvement of the Hamiltonian is not merely targeted to the description of the ground state alone. Similar ideas have been successfully explored in the context of quantum chemistry [4, 7]. 6.2 Truncation Error in the IMSRG In the IMSRG, the unitarity of the transformation provides a key diagonstic, which is that observables from exact many-body calculations should be invariant. We have shown in Equation 2.30 that the product of two normal-ordered 𝑀-body and 𝑁-body operators produces a sum of up to 𝑀 + 𝑁-body operators. The commutator in the IMSRG flow equation, Equation 3.4, produces a 69 flowing operator of total rank rank(𝜂) + rank(𝐻) − 1. In the IMSRG(2) truncation, we disregard operators for three or more particles which are induced by the commutator. The error in this truncation scheme manifests as a violation of unitarity in the IMSRG transformation, which is evident from the non-preservation of the flowing eigenvalue spectrum that we would obtain in an exact diagonalization. In cases where this loss of unitarity is especially egregious, the assumption is that the terms we truncated from the operator basis are important for maintaining a stable IMSRG flow for that system configuration. Increasing the truncation level to keep three and higher 𝐴-body operators is not feasible in realistic applications of the IMSRG, because the memory requirements (𝑂 (𝑁 4 ) in the IMSRG(2) versus 𝑂 (𝑁 6 ) in the IMSRG(3)) and computational costs (𝑂 (𝑁 6 ) in the IMSRG(2) versus 𝑂 (𝑁 9 ) in the IMSRG(3)) become prohibitive. The solution we propose in this Chapter is to recover physical information lost due to truncation by normal-ordering the IMSRG Hamiltonian with respect to an ensemble of Slater determinant reference states, thereby tuning the chosen operator basis in which we implement the IMSRG(2) flow. We will show that the impact on the cost of the IMSRG(2) flow is minimal. There is a non-negligible cost to the “mixing” optimization procedure, but we will show that improvements to the IMSRG result can still be obtained using less-than-optimal mixing parameters. 6.3 Formulation of the Reference State Ensemble In the reference state ensemble scheme, the many-body Hamiltonian is normal-ordered with respect to an ensemble of orthonormal Slater determinant configurations, {Φ1 , . . . , Φ𝑃 }. We choose a reference state 𝜓, which is a mixture of the Slater determinant configurations such that ∑︁𝑃 |𝜓⟩ = 𝑐𝑖 |Φ𝑖 ⟩ . (6.1) 𝑖=1 The density matrix 𝜌ˆ (dropping explicit indices for brevity) is ∑︁ ∑︁ 𝜌ˆ = |𝜓⟩ ⟨𝜓| = 𝑐𝑖∗ 𝑐𝑖 |Φ𝑖 ⟩ ⟨Φ𝑖 | = 𝑝𝑖 |Φ𝑖 ⟩ ⟨Φ𝑖 | . (6.2) 𝑖 𝑖 70 The trace tr 𝜌ˆ is simply the sum of the diagonal terms, i.e. ∑︁ ∑︁ tr 𝜌ˆ = 𝑝𝑖 ⟨Φ𝑖 |Φ𝑖 ⟩ = 𝑝𝑖 (6.3) 𝑖 𝑖 Expectation values of operators 𝑂 can be computed as traces in the density matrix, tr 𝜌𝑂. ˆ Consider the expectation value of operator pair 𝑎 †𝑝 𝑎 𝑞 , which generates the one-body density matrix, in 𝜓: ! ∑︁ ⟨𝜓|𝑎 †𝑝 𝑎 𝑞 |𝜓⟩ = tr 𝜌𝑎ˆ †𝑝 𝑎 𝑞 = tr 𝑝𝑖 |Φ𝑖 ⟩ ⟨Φ𝑖 | 𝑎 †𝑝 𝑎 𝑞 𝑖 ∑︁ = 𝑝𝑖 ⟨Φ𝑖 |𝑎 †𝑝 𝑎 𝑞 |Φ𝑖 ⟩ (6.4) 𝑖 ∑︁ = 𝑝 𝑖 𝜌𝑖 , 𝑖 where the factor 𝜌𝑖 is the one-body density matrix corresponding to reference state configuration |Φ𝑖 ⟩. Thus, we have shown that the reference state ensemble is a simple sum of one-body density matrices, in a basis of Slater determinants states, with terms weighted by the probability in which they might accurately represent the state targeted for decoupling by the IMSRG (e.g., the ground state). Note that for a “pure”, correlated state |Ψ⟩, the density operator 𝜌ˆ = |Ψ⟩ ⟨Ψ| would generate additional terms in the expectation value that are missing in the “mixed” state 𝜓. In this formulation, the occupation of a particular state in the reference ensemble becomes fractional (according to 𝑝𝑖 ) instead of an integer in the single Slater determinant case as in Equation 2.27. This implies that the two-body and three-body density matrices will not factor completely into one-body density matrices, which generates nonzero IDMs in Equations 2.34 and 2.35. In our reference ensemble scheme, the irreducible terms do not encode pure correlations in a quantum mechanical sense, because they do not result from the “quantum interference” terms that are present within a pure wave function density operator. Since the reference ensemble consists of a simple mixture of Slater determinants, we may interpret the irreducible terms as providing “statistical” correlation information, in a effort to simulate a truly correlated reference state. Concretely, the one-body density matrix corresponding to the reference state ensemble reads 𝑃 ∑︁ 𝝆 (E) = 𝑝 𝑖 𝝆𝑖 , (6.5) 𝑖 71 where 𝑃 is the total number of probabilities in which we expand the reference, and each 𝝆𝑖 corre- sponds to the irreducible one-body density matrix for a particular Slater determinant configuration, which is fixed by the number of hole and particle states in the single particle basis. The total number of Slater determinant configurations available to the ensemble is given by the binomial coefficient 𝑃 = 𝐶 (𝑁, 𝑁 ℎ ), where 𝑁 is the number of single particle states and 𝑁 ℎ is the number of hole states. States which are “excluded” from the ensemble correspond to weights 𝑝𝑖 = 0. In order for particle number to be conserved, we assert that the probabilities are normalized such that 𝑃 ∑︁ 𝑝𝑖 = 1. (6.6) 𝑖 The reference state ensemble adds two additional degrees of freedom to the IMSRG flow, namely the size of the ensemble and the distribution of probabilities across the ensemble. Thus, we may optimize the ensemble size and weighting scheme in order to minimize the truncation error of the IMSRG(2). Truncation error is difficult to measure without an exact method to compare with the IMSRG(2) result. Fortunately, the exact spectrum is accessible for the P3H Hamiltonian (see Chapter 4.1), which we employ to study the efficacy of the reference state ensemble. 6.4 Ensemble Optimization Scheme As a proof-of-concept, we demonstrate one technique which can be used to optimize the prob- ability weights that parameterize the reference state ensemble when an exact method is tractable for comparison. In Algorithm 6.1, we define a residual as the difference between the exact eigen- spectrum and the IMSRG(2) eigenspectrum, computed from a particular reference state ensemble characterized by the probability weight vector p. We perform a least-squares minimization of this residual under variations of p. Of course, this optimization scheme is not practical in a realistic application, where the IMSRG(2) may be expensive to complete (IMSRG emulator techniques may assist in reducing the expense, see Chapter 7.3) or the exact eigenspectrum may not be accessible due to the exponential cost of determining it. 72 Algorithm 6.1 Pseudocode implementation which demonstrates the minimization process, for which the residual involves a complete IMSRG(2) flow as well as a full configuration interaction (FCI) computation to obtain the eigenvalues of the vacuum Hamiltonian. function residual(p, y, M, 𝑚) Í p ← p/sum(p) ⊲ enforce constraint 𝑖 𝑝𝑖 = 1 z = Mp ⊲ compute reference state x = R (z, 𝛿, 𝑔, 𝑏) ⊲ compute eigenvalues from 𝐻 (∞) return x[1 : 𝑚] − y[1 : 𝑚] end function initialize reference ensemble matrix M initialize p y = F (𝛿, 𝑔, 𝑏) ⊲ compute eigenvalues from FCI repeat p ← lsq (residual(p, y, 𝑚)) until residual(p, y, M, 𝑚) < tolerance 6.5 Regularization of Generator Divergences Let us consider the IMSRG(2) flow of a P3H Hamiltonian (see Chapter 4.1) with a strong pairing strength, 𝑔 = 2.0. Let 𝑏 = 0.0. Figure 6.1 illustrates an IMSRG(2) flow of this Hamiltonian. The diagonal matrix elements, which should converge to the true eigenspectrum as the IMSRG(2) flows, diverge from this destination and cross each other unpredictably. This is an issue of the White generator for states which are close to the Fermi surface; the generator, which depends on denominators of energy differences, blows up uncontrollably when the energy differences are small (see Chapter 3.4.2). We note that although this behavior is most obvious with the White generator because of the explicit energy denominators, this issue may persist in other generator forms due to the truncation of induced terms. In this example, the pairing strength squeezes the energy states close together such that the IMSRG(2) flow cannot sufficiently decouple the target ground state from excited states, and the flow fails. In the proceeding results, we attempt to correct this failure by inserting information about the excitations of the system into the IMSRG(2) through the use of an ensemble of Slater determinants. While the use of fractional occupation factors implies that two-body and higher irreducible density matrices (IDMs) introduced in Chapter 2.6 are nonzero, as discussed above, we truncate terms containing the IDMs for now, and refer to this ensemble evaluation scheme as fractional occupation only, or FOO. 73 Diagonal spectrum flow 1-state ensemble # S.P. states = 8 | g, b = 2.00, 0.00 8 6 4 diagonal ME 2 0 2 4 6 0 1 2 3 4 5 6 7 s_vals Figure 6.1 Diagonal matrix elements of a flowing P3H Hamiltonian with 8 single particle states. The dashed horizontal lines to the right represent the true eigenspectrum for this Hamiltonian, obtained through exact diagonalization (see Chapter 2.10). The longest dashed line represents the true ground state. Figure 6.2 demonstrates the effect of the FOO ensemble scheme: We are able to cure the divergent behavior in the IMSRG(2) flow for strong pairing 𝑔 = 2.0, and the ground state energy converges to the desired result; in addition, many of the unpredictable crossings in excited states have been mitigated. One disadvantage we can immediately note from this process is that the dynamical range has doubled compared to Figure 6.1, meaning that more computational effort is required to achieve the same result (without mentioning the computational effort required for optimization). The spectrum of the flowing Hamiltonian in the FOO scheme is not necessarily invariant under changes of 𝑠 (and unitarity is still violated). However, the extent of the unitarity violation is improved compared to the single reference case. Figure 6.3 compares the spectrum evolution according to the single reference flow versus the FOO scheme flow. 6.6 Improving the Eigenspectrum In this section, we demonstrate the complexity in improving all eigenvalues and eigenstates of the IMSRG(2) evolved Hamiltonian simultaneously, using the FOO ensemble scheme from 74 Diagonal spectrum flow 5-state ensemble # S.P. states = 8 | g, b = 2.00, 0.00 8 6 diagonal ME 4 2 0 2 0 2 4 6 8 10 12 14 s_vals Figure 6.2 Diagonal matrix elements of a flowing P3H Hamiltonian with 8 single particle states. In contrast to Figure 6.1, the Hamiltonian has been normal-ordered with respect to a 5-state optimized reference ensemble. Eigenspectrum flow; single reference (left) vs 5-state ensemble (right) nholes,particles = 4,4; g, b = 2.00, 0.00 8 5 6 evolved eigenvalue evolved eigenvalue 0 4 5 2 0 10 2 0 2 4 6 0 2 4 6 8 10 12 14 s_vals s_vals Figure 6.3 Comparison of the flowing spectrum for a single reference flow to a 5-state FOO scheme flow. The red dashed lines represent the true spectrum for this particular P3H Hamiltonian. Blue lines correspond to paired eigenstates. The spectrum is not invariant under IMSRG(2) evolution in either case, but the FOO scheme flow preserves more of the original spectrum. 75 Eigenvalue error (relative to FCI) d=1.0, g=2.0, b=0.0 # S.P. states = 8 0 0 error_eig0 error_eig1 100 100 101 101 102 102 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 x1 x1 0 0 error_eig2 error_eig3 100 100 101 101 102 102 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 x1 x1 0 0 error_eig4 error_eig5 100 100 101 101 102 102 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 x1 x1 Figure 6.4 The vertical axis is the difference between the IMSRG(2) energy and the FCI energy for each eigenstate. The horizontal axis 𝑥 1 is the mixing strength, which controls the 2-state ensemble dependence on the first excited fully paired SD configuration. the previous section. “Improvement” in the spectrum is interpreted as the preservation of the eigenvalues throughout the flow, which assesses the unitarity of the IMSRG transformation. We have shown in Figure 6.2 that the reference state ensemble scheme “regularizes” the flow of matrix elements; compared to the extreme case in Figure 6.1, this accounts for an improvement in preservation of the eigenspectrum (especially, the ground state) which we observed in Figure 6.3. Let us systematically study how a 2-state FOO ensemble changes the eigenspectrum preservation for 6 fully paired1 eigenstates of the IMSRG(2) evolved Hamiltonian. The 2-state FOO ensemble, in this case, reads 𝝆 (E−2) = (1 − 𝑥1 ) 𝝆0 + 𝑥 1 𝝆1 , (6.7) 1 For the P3H Hamiltonian in question, the pair-breaking contribution is zero. Thus, we worry only about the accuracy of the eigenstates accessed by the pairing contribution. 76 where 𝑥 1 is a fraction which controls the proportion of the ground state SD configuration compared to the first excited state SD configuration. Note that when 𝑥1 = 0, we recover the usual single SD reference state. In the FOO ensemble scheme, the weights (1 − 𝑥 1 ), 𝑥1 directly affect the fractional values in the occupation factors, and nothing else. Figure 6.4 shows the error curves for each fully paired eigenstate, for the same Hamiltonian as in Figures 6.1 and 6.2 using the 2-state reference ensemble. The solid horizontal line in each plot denotes the zero crossing, where the error in the energy is zero. These plots demonstrate that, for a single mixing parameter, there is no mixing for which the total error in the IMSRG(2) eigenspectrum is zero. This result suggests that we must make concessions based on the goal of the reference state ensemble scheme; tuning the ensemble to improve a single eigenstate, e.g. the ground state, does not necessarily guarantee all eigenstates will be improved directly proportionally. In the 2-state ensemble, there is maybe no good reason to expect improvements for high-lying excited states if we only use ground-state and first-excited-state Slater determinant configurations. Seeking improvement in high-lying excitations motivates the use of more basis configurations in the ensemble, at the cost of a more complicated multi-parameter optimization problem. Figure 6.7 also displays a sweeping range of control over error in every IMSRG(2) eigenstate for only a single mixing parameter in a 2-state ensemble. The 2-state ensemble begins to see improvement over the single reference beyond around 𝑥1 = 0.1, before which we observe a steep increase in error around 𝑥 1 = 0.08. In some way, this singularity reveals the topology of the optimization problem. The steep increase shows how much we must weight the addition of the excited state in order to decouple the ensemble result from the single reference result, such that the spectrum no longer diverges. These results suggest that we may be able to alleviate particular crossings with a simple 2-state ensemble. 6.7 Ensemble Consistency with Non-Targeted Observables In the previous sections, we have shown that reference ensembles are capable of improving the IMSRG(2) flow of the Hamiltonian—in other words, the reference ensemble can reduce the error 77 in energy eigenvalues, or even improve the entire spectrum robustly. However, most of the time, we wish to also compute the expectation values of observables other than the energies, hence we must ensure that the use of the reference state ensemble does not spoil their description. We may evolve any additional observable 𝑂 alongside the Hamiltonian according to Equation 3.4 [2], 𝑑𝑂 (𝑠) = [𝜂(𝑠), 𝑂 (𝑠)] , (6.8) 𝑑𝑠 where 𝜂(𝑠) is the IMSRG generator constructed from the Hamiltonian. Since the space for storing the full flow history of 𝜂 is prohibitive in realistic applications, we must concatenate the operator flow equations in the same normal-ordered operator basis as the Hamiltonian to the Hamiltonian ODE system. Concatenation of each additional operator implies the addition of another full- size version of the ODE system, which restricts the simultaneous evolution to only a couple of observables in practice. One solution for this issue is the Magnus [5, 1] formulation of the IMSRG demonstrated by Morris, Parzuchowski, and Bogner [6] (also see Hergert et al. [2] and Chapter 10 of Hjorth-Jensen, Lombardo, and Kolck [3]). The Magnus method casts the unitary transformation to a true exponential, 𝑈 (𝑠) = 𝑒 Ω(𝑠) , and then reformulates the flow equations to evolve the Magnus operator Ω(𝑠) instead of the Hamiltonian itself. Future work will implement in the reference state ensemble in the Magnus-IMSRG, but has not been implemented for the results in this study. One disadvantage of the rather abstract P3H toy model is the lack of compelling observables with which to gather relevant data. The most interesting candidate beside 𝐻 itself is the total spin- squared operator 𝑆ˆ2 . We will explore how the use of a reference state ensemble that is optimized for 𝐻 impacts the expectation value ⟨𝑆 2 ⟩; in this way, we show that the reference ensemble flow is “safe” in that important physical information in simultaneous observables is not destroyed. The second-quantized spin-squared operator 𝑆ˆ2 reads ∑︁ 𝑆ˆ2 = ⟨𝑝| 𝑠ˆ2 |𝑞⟩ 𝑎 †𝑝 𝑎 𝑞 𝑝𝑞 (6.9) 1 ∑︁ + ⟨𝑝𝑞| 𝑠ˆ (1) · 𝑠ˆ (2) |𝑟 𝑠⟩ A 𝑎 †𝑝 𝑎 †𝑞 𝑎 𝑠 𝑎𝑟 . 4 𝑝𝑞𝑟 𝑠 The 𝑆ˆ2 matrix elements are constructed in the single particle basis states of the P3H model, with quantum numbers 𝑝, which denotes the energy level, and 𝜎𝑝 , which denotes the spin direction of 78 the particle at level 𝑝. The one-body matrix elements of 𝑆ˆ2 , in units of ℏ2 , read 3 ⟨( 𝑝, 𝜎𝑝 )| 𝑠ˆ2 |(𝑞, 𝜎𝑞 )⟩ = 𝛿 𝑝𝑞 𝛿𝜎𝑝 𝜎𝑞 . (6.10) 4 Using ladder operators 𝑠ˆ− and 𝑠ˆ+ , the two-body part of the operator can be expressed as 1  (1) (2)  (1) (2) 𝑠ˆ · 𝑠ˆ = 𝑠ˆ+ 𝑠ˆ− + 𝑠ˆ− 𝑠ˆ+ + 𝑠ˆ𝑧(1) 𝑠ˆ𝑧(2) . (1) (2) (6.11) 2 The ladder operator terms in Equation 6.11 correspond to a mutual “spin-flipping” interaction, whereas the spin projection 𝑠ˆ𝑧(1) 𝑠ˆ𝑧(2) term is structurally similar to the one-body matrix element. Introducing the two-body antisymmetrized product state 1 |( 𝑝, 𝜎𝑝 )(𝑞, 𝜎𝑞 )⟩ A = √ (|( 𝑝, 𝜎𝑝 )(𝑞, 𝜎𝑞 )⟩ − |(𝑞, 𝜎𝑞 )( 𝑝, 𝜎𝑝 )⟩), (6.12) 2 we can compute the two-body matrix element corresponding to the 𝑠ˆ (1) · 𝑠ˆ (2) term in Equation 6.9. In Figure 6.5 we present the evolution of the expectation value ⟨𝑆 2 ⟩ in the ground eigenstate of the Hamiltonian, compared to evolution of the ground state energy. The vertical axis measures the observables’ difference from the starting point, so that energy and ⟨𝑆 2 ⟩ can be plotted to the same scale. The horizontal axis represents the 𝑠 evolution scale. In panels 𝐴 and 𝐵, we have normal-ordered both operators with respect to the single, standard SD ground state configuration. Note that in these cases, the 𝑆ˆ2 matrix maps to zero because the total spin is zero for the SD ground state configuration. In this way, ⟨𝑆 2 ⟩ is trivially conserved. The energy behaves as expected, in that the strongly interacting Hamiltonian in panel 𝐵 results in a divergent flow. In panel 𝐶, we have normal-ordered both operators with respect to an optimized ensemble. Now that there are 𝑆 ≠ 0 SD configurations mixed into the ensemble, we see that ⟨𝑆 2 ⟩ is no longer conserved. However, we see that the evolution of ⟨𝑆 2 ⟩ closely resembles that of the energy, meaning that the character of the 𝑆ˆ2 has been preserved under the reference ensemble flow. This study suggests that the reference state ensemble is safely preserves simultaneous observables in the IMSRG flow. 6.8 Using Exact Eigenstates as References We conclude the discussion in this chapter by considering the impact of using an exact eigenstate of 𝐻 for setting up the normal ordered operator basis for the IMSRG. In this case, we expect the 79 Single standard reference Single standard reference d=1.0, g=0.5, b=0.0 d=1.0, g=2.0, b=0.1 A # S.P. states = 8 B # S.P. states = 8 diff = 0 Difference from s=0 Difference from s=0 8 E_gs(s=0)-E_gs(s) 0.003 SS(s=0)-SS(s) 6 diff = 0 0.002 E_gs(s=0)-E_gs(s) SS(s=0)-SS(s) 4 0.001 2 0.000 0 0 2 4 6 8 10 12 14 16 0.0 0.2 0.4 0.6 0.8 1.0 1.2 s s Optimized reference ensemble d=1.0, g=2.0, b=0.1 C # S.P. states = 8 diff = 0 Difference from s=0 0.2 E_gs(s=0)-E_gs(s) SS(s=0)-SS(s) 0.0 0.2 0.4 0.6 0 2 4 6 8 10 12 14 s Figure 6.5 Evolution of ⟨𝑆 2 ⟩ and the ground state energy 𝐸 (𝑠) for two P3H model Hamiltonians and two reference schemes. The vertical axis measures the change in the zero-body part of the evolving operator from the starting point 𝑠 = 0, while the horizontal axis plots the 𝑠 evolution range. Panel 𝐴 features a weakly interacting P3H Hamiltonian where energy converges, according to a single standard ground state SD configuration reference. Panel 𝐵 features a strongly interacting P3H Hamiltonian where energy diverges, according to a single standard ground state SD configuration reference. Panel 𝐶 features a strongly interacting P3H Hamiltonian where energy converges, according to an optimized ensemble reference. 80 zero-body part of the operator to be identical to the eigenvalue associated with the reference state, and it should be completely insensitive to the IMSRG truncation. Since the eigenstate is correlated, the IDMs will be non-zero, and we will perform the evolution in the MR-IMSRG(2) scheme (see Appendix B). We will use generalized form of the White-ArcTan generator (Chapter 3.4.2), as well as the Brillouin generator (Chapter 3.4.3), to generate the flow. In Figure 6.6, we show the eigenvalue flow corresponding to the lowest four eigenstates of a weak-pairing P3H Hamiltonian and strong-pairing/weak-pair-breaking P3H Hamiltonian. Note that in both cases, the first and second excited states are degenerate in energy. The color coding indicates which exact eigenstate was used to normal order 𝐻. In the weak-pairing case, we find that all eigenstates remain invariant under evolution if we use the associated occupations and IDMs for constructing the normal ordered operator basis, and we find the expected lack of sensitivity to the MR-IMSRG truncation. Interestingly, the (fully paired) ground state remains invariant regardless of the exact eigenstate we use for the normal ordering. The third excited state, which is fully paired as well, also remains invariant if we use any of the excited states to construct the normal ordered operator basis, while a basis built from the ground state is not able to capture the evolution of this state completely. For the first and second excited states, evolving with the fully paired references only leads to minor deviations from the exact eigenvalue. In the strong-pairing/weak-pair-breaking case, we find that the flow breaks down for all cases except in the case where Hamiltonian is normal-ordered with respect to the ground state. This suggests that for strongly interacting systems, a single exact eigenstate does not provide enough information to account for these induced terms which have been truncated. Note especially the “closeness” of the second and third energies in the strongly interacting P3H Hamiltonian–the eigenvalues in these state seem to be switching places in a level crossing. In Figure 6.7, we systematically study how truncations of terms depending on the IDMs affect the unitarity of the ground state eigenvalue’s evolution. Here, we examine the flow generated by the Brillouin generator and normal order with respect to the exact ground state. In both cases, the inclusion of the fractional occupation numbers alone yields a converged flow, but the deviation 81 Eigenspectrum (lowest four) vs flow parameter Eigenspectrum (lowest four) vs flow parameter # S.P. states = 8 # S.P. states = 8 g, b = 0.50, 0.00 g, b = 1.50, 0.10 3.5 2 3.0 evolved eigenvalue evolved eigenvalue 1 2.5 N-O ground state 0 N-O ground state N-O 1st excited state (unpaired) N-O 1st excited state (unpaired) N-O 2nd excited state (unpaired) N-O 2nd excited state (unpaired) 2.0 N-O 3rd excited state (paired) N-O 3rd excited state (paired) 1 1.5 2 0 5 10 15 20 0 5 10 15 20 s_vals s_vals Figure 6.6 Evolution of the lowest four eigenvalues in two P3H models, varying the exact eigenstate against which the Hamiltonian has been normal-ordered. The flow was evaluated using a White- ArcTan generator. In the left pane, we evolve a weak pairing Hamiltonian with no pair-breaking. In the right pane, we evolve a strong pairing Hamiltonian with weak pair-breaking. from the exact eigenvalue can be substantial, and in the strongly interacting case it is difficult to tell whether the evolution is completely leveled off. In the weakly-interacting case, inclusion of the two-body IDM in the Brillouin generator greatly reduces the violation of unitarity, and inclusion of the three-body IDM guarantees the invariance of the eigenvalue. In the strongly-interacting case, the inclusion of the two-body IDM terms on top of the fractional occupation numbers actually causes the flow to “diverge”—since the Brillouin generator is effec- tively the gradient of the energy under unitary variations, a likely explanation is that we overshoot the stationary point because the gradient is unbalanced. The inclusion of the three-body IDM fixes this issue, hence the takeaway is that IDM truncations in the generator need to be carefully considered. Since the IDMs encode correlations in the reference state that are coupling to the correlations treated by the MR-IMSRG during the evolution, it is not surprising that the balance of the two- and three-body IDMs is more delicate in the strongly interacting case. 82 Figure 6.7 Flow of the ground state eigenvalue when normal ordering with respect to the ground state for two P3H Hamiltonians, with no-pairing breaking. The left pane is a weak pairing case, where the right pane is the strong pairing case. The truncation level of the IDMs in the generator and the MR-IMSRG(2) flow has been varied. 6.9 Outlook In this chapter, we have presented the reference state ensemble scheme for reducing truncation error in the IMSRG(2) flow. We have shown that by introducing information about potential excitations of the system into the IMSRG(2) flow via the reference state ensemble, we can reduce the violation of the transformation’s unitarity due to truncations. This is evident from the improvement in the preservation of the spectrum of evolved P3H Hamiltonians which have been normal-ordered with respect to reference ensembles. We have also shown that the use of reference ensembles that are optimized for the Hamiltonian does not seem to significantly distort observables which are evolved alongside the Hamiltonian, although this matter will have to be revisited in more realistic applications. In the future, we hope to improve the optimization process in Algorithm 6.1 in several ways. For one, the exact spectrum will be inaccessible in applications to realistic systems, so we will need to to conceive of different characteristics for the loss of unitary in observables, e.g., partial 83 traces. Or, perhaps, we could implement a neural network optimization strategy to “learn” the best ensemble mixing from a large set of training points. Of course, this method is also constrained by the computational effort for performing the full (MR-)IMSRG evolutions, as well as the richness of viable evolutions in the parameter surface spanned by the mixing parameters. We may be able to leverage the emulator techniques discussed in Chapter 7 for generating comprehensive training sets for reference ensemble optimization. 84 BIBLIOGRAPHY [1] S. Blanes et al. “The Magnus expansion and some of its applications”. In: Physics Reports 470.5 (2009), pp. 151–238. issn: 0370-1573. doi: https://doi.org/10.1016/j.physrep.2008. 11.001. url: https://www.sciencedirect.com/science/article/pii/S0370157308004092. [2] H. Hergert et al. “The In-Medium Similarity Renormalization Group: A novel ab initio method for nuclei”. In: Physics Reports 621 (Mar. 2016), pp. 165–222. issn: 0370- 1573. doi: 10.1016/j.physrep.2015.12.007. url: http://dx.doi.org/10.1016/j.physrep.2015. 12.007. [3] Morten Hjorth-Jensen, Maria-Paola Lombardo, and Ubirajara van Kolck. An advanced course in computational nuclear physics: bridging the scales from quarks to neutron stars. Springer, 2017. [4] Chenyang Li and Francesco A. Evangelista. “Driven similarity renormalization group for excited states: A state-averaged perturbation theory”. In: The Journal of Chemical Physics 148.12 (Mar. 2018), p. 124106. issn: 0021-9606. doi: 10.1063/1.5019793. eprint: https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5019793/16710196/124106\ _1\_online.pdf. url: https://doi.org/10.1063/1.5019793. [5] Wilhelm Magnus. “On the exponential solution of differential equations for a linear op- erator”. In: Communications on Pure and Applied Mathematics 7.4 (1954), pp. 649–673. doi: https://doi.org/10.1002/cpa.3160070404. eprint: https://onlinelibrary.wiley.com/doi/ pdf/10.1002/cpa.3160070404. url: https://onlinelibrary.wiley.com/doi/abs/10.1002/cpa. 3160070404. [6] T. D. Morris, N. M. Parzuchowski, and S. K. Bogner. “Magnus expansion and in-medium similarity renormalization group”. In: Phys. Rev. C 92 (3 Sept. 2015), p. 034331. doi: 10. 1103/PhysRevC.92.034331. url: https://link.aps.org/doi/10.1103/PhysRevC.92.034331. [7] Thomas J. Jr. Watson and Garnet Kin-Lic Chan. “Correct Quantum Chemistry in a Minimal Basis from Effective Hamiltonians”. In: Journal of Chemical Theory and Computation 12.2 (2016). PMID: 26756223, pp. 512–522. doi: 10.1021/acs.jctc.5b00138. eprint: https://doi.org/10.1021/acs.jctc.5b00138. url: https://doi.org/10.1021/acs.jctc.5b00138. 85 CHAPTER 7 KOOPMAN OPERATOR FOR THE IN-MEDIUM SIMILARITY RENORMALIZATION GROUP 7.1 Introduction In this chapter, we present the results of several algorithms developed from Koopman operator approximation techniques, namely the Dynamic Mode Decomposition (DMD), to emulate the solution to the IMSRG(2) flow (Equation 3.18). There are two algorithms which work in tandem to generate a complete emulator to the parametric IMSRG problem (i.e. emulate IMSRG solution varying couplings): 1) a standard DMD formula, which we show can emulate the solution to a particular IMSRG(2) flow [4, 3] and 2) an algorithm which interpolates DMD operators across a parametric coupling surface which characterizes the model Hamiltonian we wish to evolve via IMSRG(2) [6]. 7.2 The Koopman Operator The Koopman operator [7] is an infinite-dimensional linear operator which acts on a mea- surement subspace of a nonlinear dynamical system. The system itself may be finite or infinite- dimensional. The spectral decomposition of this operator characterizes the behavior of the nonlinear dynamical system [3, 4]. In principle, the measurement subspace is continuous in time; however, for our applications we concern ourselves with the discrete-time Koopman operator which acts on a discrete measurement subspace. Consider a measurement function 𝑔 which evaluates a system 𝑥, where 𝑥 is interpreted as a vector of state variables, at some point in time 𝑥 𝑘 . Forward time-steps in the system are governed by a flow map F, which propagates the system forward in time such that, F(𝑥 𝑘 ) = 𝑥 𝑘+1 . (7.1) The Koopman operator is defined by the composition of 𝑔 and F [4], K𝑔 = 𝑔 ◦ F, (7.2) 86 which shows that the Koopman operator propagates the measurement function forward in time, so that [4], K𝑔(𝑥 𝑘 ) = 𝑔(𝐹 (𝑥 𝑘 )) = 𝑔(𝑥 𝑘+1 ). (7.3) In other words, the Koopman operator characterizes the flow map F, under which 𝑥 evolves, via measurements of the system by 𝑔. In the eigenbasis of K [4], K𝜙(𝑥 𝑘 ) = 𝜆𝜙(𝑥 𝑘 ) = 𝜙(𝑥 𝑘+1 ), (7.4) which shows that evaluation of the Koopman eigenfunctions characterize the dynamics of the system. In other words, the dynamics of the nonlinear system are fully captured by an infinite- dimensional linear operator (the Koopman operator) that acts on the space of measurements of the system. Since the Koopman operator is infinite-dimensional, we must find techniques which can ap- proximate the Koopman operator in a finite dimension. The technique which we use in this chapter is called Dynamic Mode Decomposition, or DMD. 7.3 Dynamic Mode Decomposition for IMSRG The DMD is a finite-dimensional, finite-measurement, approximation to the discrete Koopman operator for a nonlinear dynamical system. It was originally developed by Schmid [12] in the context of fluid dynamics. The DMD is formulated as the best-fit linear operator which connects successive, time-ordered measurements of the target system [12, 11]; as a result, Singular Value Decomposition (SVD) is most often the technique used to solve for the DMD operator, due to its computationally efficient scaling and flexibility in low-rank truncated spaces. The DMD begins with 𝑁 successive observations from the evolving dynamical system, orga- nized into offset data matrices 𝑋 and 𝑋 ′, where 𝑋 ′ is one step forward relative to 𝑋. In matrix form,     | |  | |      ′     𝑋 = 𝑥1 · · · 𝑥 𝑁−1  ;  𝑋 = 𝑥 2 · · · 𝑥 𝑁  .  (7.5)     | |  | |          87 Algorithm 7.1 Pseudocode implementation which demonstrates the standard algorithm for DMD. collect 𝑁 observations from evolving dynamical system, taking steps of fixed width 𝑑𝑡 organize offset matrices 𝑋 and 𝑋 ′, with shapes 𝑀 × 𝑁 − 1 solve 𝑋 † (via SVD) 𝐴 = 𝑋′𝑋† compute eigendecomposition 𝐴Φ = ΛΦ compute DMD mode amplitudes 𝑏 = Φ† 𝑥 0 expand dynamical state 𝑥 𝑘 = ΦΛ𝑘−1 𝑏 ⊲ for 𝑘 > 0, For simplicity, we assert that each observation be taken at a regular intervals 𝑑𝑡 in the dynamical variable. The DMD operator 𝐴 shifts 𝑋 one step forward in time (or whatever the dynamical variable may be), i.e. 𝐴𝑋 = 𝑋 ′ . (7.6) We may proceed to solve for the “best-fit” 𝐴 in a least-squares sense, obtaining 𝐴 = 𝑋 ′ 𝑋 †, (7.7) where † denotes the Moore-Penrose inverse (pseudoinverse) of 𝑋. The Moore-Penrose inverse minimizes the L2 loss function ∥ 𝑋 ′ − 𝐴𝑋 ∥ 2 such that 𝑋 † = (𝑋 𝑇 𝑋) −1 𝑋 𝑇 , (7.8) where we must have that (𝑋 𝑇 𝑋) is invertible. Standard Dynamic Mode Decomposition The standard DMD algorithm is written, in pseu- docode, in Algorithm 7.1. In the final line of Algorithm 7.1, the state vector 𝑥 𝑘 represents the state of the dynamical system at step 𝑘, where 𝑘 = 𝑡 𝑘 /𝑑𝑡. The DMD eigenfunctions, or modes, are denoted Φ and the DMD eigenvalues are denoted Λ. The vector 𝑥 0 represents the initial state of the dynamical system; the product Φ† 𝑥 0 , usually denoted 𝑏, is the DMD amplitudes which represent the DMD eigenfunctions evaluated in the initial conditions. Algorithm 7.1 explains how the dynamical system is expanded in DMD modes, according to the DMD expansion [8, 3, 4, 12], written 𝑥 𝑘 = ΦΛ𝑘−1 𝑏, 𝑘 > 0. (7.9) 88 We may translate this expansion to continuous time by converting the diagonal DMD eigenvalue matrix Λ to its continuous analog Ω, with diagonal elements 𝜔𝑖 = ln (𝜆𝑖 )/𝑑𝑡, so that the continuous expansion reads [4], 𝑥(𝑡) = Φ exp (Ω𝑡)𝑏. (7.10) In order to fully characterize the dynamical system, the necessary information that must be extracted from DMD are the modes Φ and the eigenvalues Λ. Adjacent objects, such as the amplitudes 𝑏, are inferred or directly computed from Φ and Λ. Thus, we refer to a DMD decomposition of a dynamical system as the pair (Φ, Λ). Projected Dynamic Mode Decomposition One prohibiting issue with the standard DMD al- gorithm is that the memory and computational cost scales heavily with the number of dynamical variables 𝑀. In our case, these would be the coefficients of the IMSRG(2) Hamiltonian. The DMD operator in Algorithm 7.1 is a dense matrix with shape 𝑀 × 𝑀; for a typical nuclear many-body calculation using IMSRG(2), 𝑀 might be on the order 106 or greater. The viability of the standard DMD algorithm is quickly challenged for practical uses of the IMSRG(2). Instead, we leverage the flexibility of the SVD to solve the projected DMD, which solves for the DMD operator in a reduced measurement subspace. Beginning with an SVD on the data matrix 𝑋, 𝑋 ≈ 𝑈Σ𝑉 ∗ , (7.11) we can truncate the number of singular values used in the reconstruction at 𝑟 while minimizing the reconstruction error, assuming that the singular values are exponentially decaying [8]. The truncated SVD expression then reads 𝑋trunc ≈ 𝑈𝑟 Σ𝑟 𝑉𝑟∗ . (7.12) Note that the shapes of 𝑋 and 𝑋trunc are equivalent, so we assume that 𝑋 ≈ 𝑋trunc up to some negligible truncation error (due to the exponentially decaying singular values), which is controllable by choosing an appropriate 𝑟. Then we may interpret the matrix 𝑈𝑟 as a transformation from the full measurement space to the reduced measurement subspace. 89 The reduced DMD operator is the projection of the full DMD operator onto the reduced SVD modes [8], 𝐴𝑟 = 𝑈𝑟∗ 𝐴𝑈𝑟 , (7.13) where 𝐴𝑟 is the 𝑟 ×𝑟 projected DMD operator onto the reduced space. From here, we may substitute the expression for the full DMD operator 𝐴, Equation 7.7, as well as the SVD on 𝑋, to arrive at 𝐴𝑟 = 𝑈𝑟∗ (𝑋 ′ 𝑋 † )𝑈𝑟 = 𝑈𝑟∗ 𝑋 ′ (𝑈𝑟 Σ𝑟 𝑉𝑟∗ ) †𝑈𝑟 (7.14) = 𝑈𝑟∗ 𝑋 ′𝑉𝑟 Σ𝑟−1 . The reduced DMD operator 𝐴𝑟 and the full DMD operator 𝐴 necessarily share up to 𝑟 eigen- values. Thus, the eigendecomposition of 𝐴𝑟 reads 𝐴𝑟 𝑊 = 𝑊Λ. (7.15) Tu et al. [13] showed that the full DMD eigenvectors Φ can be recovered from the reduced DMD eigenvectors 𝑊 via Φ = 𝑋 ′𝑉𝑟 Σ𝑟−1𝑊 . (7.16) At this point, the full DMD decomposition (Φ, Λ) has been recovered from the reduced measurement subspace. The dynamical system can then be expanded in the DMD modes according to Equations 7.9 or 7.10. Algorithm 7.2 summarizes the procedure for computing a reduced DMD decomposition. Note that the computational effort for computing the eigendecomposition of the DMD operator has been significantly reduced by solving for the reduced operator 𝐴𝑟 ; the projection softens the scaling from 𝑂 (𝑀 3 ) in the full measurement space, to 𝑂 (𝑟 3 ) in the reduced measurement subspace, where 𝑟 ≪ 𝑀. 7.4 Emulating the IMSRG Flow In this section, we justify how the DMD can be used to “learn” the unitary transformation which propagates the IMSRG flow. We have studied only IMSRG up to 2-body truncations, i.e. IMSRG(2), although we emphasize that the emulation procedure is consistent for IMSRG(𝑁). 90 Algorithm 7.2 Pseudocode implementation which demonstrates the algorithm for projected DMD. collect 𝑁 observations from evolving dynamical system, at fixed width 𝑑𝑡 organize offset matrices 𝑋 and 𝑋 ′, with shapes 𝑀 × 𝑁 − 1 compute truncated SVD on 𝑋 ≈ 𝑈𝑟 Σ𝑟 𝑉𝑟∗ compute 𝐴𝑟 = 𝑈𝑟∗ 𝑋 ′𝑉𝑟 Σ𝑟−1 compute eigendecomposition 𝐴𝑟 𝑊 = 𝑊Λ project back to full space Φ = 𝑋 ′𝑉𝑟 Σ𝑟−1𝑊 or Φ = 𝑈𝑟 𝑊 compute DMD mode amplitudes 𝑏 = Φ† 𝑥 0 expand dynamical state 𝑥 𝑘 = ΦΛ𝑘−1 𝑏 ⊲ for 𝑘 > 0, In Chapter 3.2, Equation 3.1, we expressed the continuous, unitary, 𝑠-dependent transformation 𝑈 (𝑠) as the transformation which drives the target Hamiltonian toward a desired form (e.g., to the minimal decoupling of an eigenvalue), and in the standard IMSRG approach, we “solve” for this transformation implicitly by propagating the system of flow equations for the Hamiltonian. We can also think of the transformation in the form 𝐻 (𝑠) = U (𝑠)𝐻 (0), (7.17) where U (𝑠) is an unknown “super operator” that acts on the algebra of second-quantized operators, and encodes the desired continuous unitary transformation. The expression for the DMD expansion in Equation 7.10 provides a description of exactly this form. Replacing 𝑥 → 𝐻 so that the coefficients of 𝐻 (𝑠) becomes the state variables tracked in 𝑥, we can express U (𝑠) in terms of the Koopman modes and frequencies as 𝐻 (𝑠) = U (𝑠)𝐻 (0) = {Φ exp (Ω𝑡)Φ† }𝐻 (0) . (7.18) Note that the product {Φ exp (Ω𝑡)Φ† } is unitary, provided that the columns of Φ are sufficiently linearly independent. We may take for granted that this transformation provides the desired evolution behavior, simply because the transformation is informed by measurement data of the IMSRG flow. Of course, this assumption is predicated on good measurement data, meaning the IMSRG flow being measured is well-behaved to begin with. In this way, we consider the DMD a data-driven, non-intrusive emulator for the IMSRG—“data-driven” in the sense that the DMD is built from, or informed by, data measured from the evolving system, and “non-intrusive” in the sense that DMD does not modify the IMSRG method itself. 91 DMD eigenvalue dynamics 1.0 0.8 P3H model (d,g,b) = (0.25,0.30,0.2) 100 DMD observations, SVD rank = 6 0.6 0.4 0.2 0.0 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 s Figure 7.1 DMD mode dynamics in 𝑠, for a rank-6 DMD expansion of a strongly interacting P3H model. Each curve correspond to exponentials 𝑒 𝜔1 𝑠 through 𝑒 𝜔6 𝑠 in Equation 7.19, ordered from slowest to fastest frequency. Note that one mode in particular remains approximately constant in 𝑠. Equation 7.18 reveals an interesting feature of DMD emulation in the IMSRG context. Prac- tical implementation of DMD emulation will treat 𝐻 (𝑠) as a one-dimensional vector which contains all three pieces of the normal-ordered operator in the IMSRG(2) truncation, 𝐻 (𝑠) = [𝐸 (𝑠), 𝑓 (𝑠), Γ(𝑠)] 𝑇 . In this way, the rows of Φ correspond to the elements of the vector 𝐻 (𝑠). For example, we need only the first row of Φ to reconstruct the flow of 𝐸 (𝑠), 𝐸 (𝑠) = 𝜙11 e𝜔1 𝑠 𝑏 1 + 𝜙12 e𝜔2 𝑠 𝑏 2 + · · · + 𝜙1𝑟 e𝜔𝑟 𝑠 𝑏𝑟 . (7.19) For the energy to be convergent, the DMD modes must be convergent. This is readily achieved with decaying mode frequencies, 𝜔 𝑘 —and, in fact, this is what we observe from from DMD decomposition of IMSRG flow data. Figure 7.1 displays the exponentials which drive the IMSRG dynamics in the DMD reconstruction of a strongly interacting P3H model. If we take 𝑠 to infinity, DMD expansion will take the Hamiltonian to zero. In some cases, it may be beneficial or necessary to define a termination constraint in the expansion. Figure 7.1 is typical in IMSRG data in the sense that the slowest decaying mode, where the 𝜔’s are ordered from slowest to fastest, tends to “look constant” in the range of relevant dynamics. In this case, we may approximate 𝜔1 ≈ 0 so 92 that 𝑒 𝜔1 ≈ 1. In this way, the DMD expansion will converge to a nonzero result as 𝑠 increases. Formally, we might view this process as multiplying through the expansion a common factor 𝑒 −𝜔1 𝑠 , to arrive at 𝐸 (𝑠) = 𝜙11 e0 𝑏 1 + 𝜙12 e (𝜔2 −𝜔1 )𝑠 𝑏 2 + · · · + 𝜙1𝑟 e (𝜔𝑟 −𝜔1 )𝑠 𝑏𝑟 . (7.20) Note that if 𝜔1 is truly small, then for 𝜔 𝑘 where 𝑘 > 0, 𝜔 𝑘 − 𝜔1 ≈ 𝜔 𝑘 . In this special case that we wish to calculate the theoretically converged result as 𝑠 → ∞, note that 𝐸 (𝑠 → ∞) ≈ 𝜙11 𝑏 1 , (7.21) where the converged energy is simply the first element of Φ multiplied with the first DMD mode amplitude. We refer to this as the no-decay approximation or ND approximation. The approximation may be extended to the full Hamiltonian, evolved to 𝑠 → ∞, by multiplying with the entire DMD eigenvector corresponding to 𝜔1 . Now we may characterize the efficacy of forecasting the IMSRG flow using DMD. There are two primary parameters that can be varied to tune emulation performance, namely the number of observations of the evolving system, 𝑁, and the SVD truncation rank 𝑟. Increasing 𝑁 provides the DMD emulation with more context, in the sense that the dynamics are more constrained to the target problem. Increasing truncation rank 𝑟 improves the fidelity of the reduced DMD operator in use, in the sense that more SVD directions are available to the dynamical expansion. However, increasing 𝑟 too much may introduce noise to the expansion in the form of singular values which are close to zero, causing instability in the dynamical expansion. One must tune these parameters for the optimal performance in the target problem. In the IMSRG context, we want to minimize the number of observations (i.e. polynomial-cost ODE solver steps) while also minimizing the error of the forecasted evolved energy. The DMD forecast is most useful if we can stop observing the flow before the convergence elbow in the energy (see e.g. the top row of Figure 7.4, which plots the IMSRG-evolving energy). Beyond this elbow, the IMSRG is essentially solved and forecasting is no longer useful. The 𝑠 value where the convergence elbow begins and ends depends on the chosen flow generator; in general, however, we 93 10 Nh 4 6 8 8 10 TIMSRG/TDMD 14 6 16 4 2 10 15 20 25 30 35 Number of DMD observations Figure 7.2 The ratio of TACO IMSRG(2) wall time to DMD emulation time. The total wall time for DMD includes the time cost for calculating 𝑥 number of DMD observations. The hue represents the single particle basis size, with darker color corresponding to larger bases. The red dotted line is where the ratio equals 1. find that taking 20-30 observations, at a step-width characteristic to the generator, provides enough observations to safely converge to a reasonably accurate energy. This number of observations also tends to fall close to the onset of convergence. As a first example, we return to the P3H model in Chapter 4.1, calculated using the tcimsrg code from Chapter 5.6. In Figure 7.2, we show the speedup of DMD emulation over the TACO implementation of the IMSRG(2). The hue represents the size of the single particle basis. The horizontal axis reflects the number of observations used to the build the DMD operator; note that a larger measurement matrix means a larger SVD system to compute, which is why speedups decrease for increasing observations. These results show sizeable performance gains using DMD emulation. Figure 7.3 shows that these performance gains come at minimal cost to accuracy. Motivated by the cost-savings from DMD emulation in the P3H model, we continue with an example involving an IMSRG(2) calculation of 16 O, computed with a SRG-softened NN3 LO(500) interaction, and study the accuracy of the DMD forecast on the number of observations. Our results are shown in Figure 7.4. The observation steps were taken in steps of 𝑑𝑠 = 0.05, such that the 94 0.00 (EDMD EIMSRG)/EIMSRG 0.02 Nh 0.04 4 6 8 10 0.06 14 16 10 15 20 25 30 35 Number of DMD observations Figure 7.3 The relative error in the evolved energy versus number of DMD observations. We use the same hue scheme as Figure 7.2. The error is controllable by the number of DMD observations, but generally falls on the order of less than a percent no matter the size of the measurement matrix. max observed s at 1.0 corresponds to 20 observation points. We see from the bottom row of plots that the relative error in the evolved energy improves quickly as the maximum observation point approaches the elbow of convergence; beyond this region, the error improvement saturates. In this way, observations beyond the elbow are not technically relevant for obtaining a good DMD decomposition. Figure 7.5 provides more detail on how error in the predicted energy interacts with the number of observations and reduced DMD truncation rank. For fewer observations, the SVD is more sensitive to noise—in other words, directions which correspond to near-zero singular values appear at lower truncation levels, which causes instability in the DMD result. This is evident from the way relative error increases sharply beyond 𝑟 = 10 in the row corresponding to 20 observations, for example. Figure 7.5 suggests that the most robust, and effective, way to reduce DMD error is to increase the number of observations provided to the operator. In order for DMD emulation of the IMSRG flow to be effective, we must balance the amount of observations with the error that can be tolerated for the target problem. Figure 7.6 plots the Frobenius error relative to the IMSRG result for DMD emulation of the 95 A1 A2 135 DMD reconstruction 135 DMD reconstruction 140 IMSRG evolution 140 IMSRG evolution ND asymptotic result ND asymptotic result 145 Max observed s = 1.00 145 Max observed s = 1.00 E(s) 150 O16, eMax = 6, = 20 MeV 150 O16, eMax = 6, = 20 MeV 155 Max extrapolated s 155 Max extrapolated s 160 160 0 5 10 15 20 0 5 10 15 20 s s 0.25 B1 0.25 B2 Full expansion Full expansion 0.20 ND asymptotic result 0.20 ND asymptotic result Relative error 0.15 0.15 Max extrapolated s = 10.00 Max extrapolated s = 20.00 0.10 0.10 s = 1.00 s = 1.00 0.05 0.05 0.00 0.00 0.6 0.8 1.0 1.2 1.4 0.6 0.8 1.0 1.2 1.4 Max observed s Max observed s Figure 7.4 DMD forecast results for 16 O, in an eMax = 6 harmonic oscillator basis with oscillator width ℏ𝜔 = 20 Mev, with SRG-softened NN3 LO(500) interaction. The top row plots evolution of the ground state energy, according to the IMSRG as well as DMD reconstruction of the IMSRG. The vertical dashed line indicates where direct IMSRG observations stop. The bottom row plots the relative error in the energy at 𝑠 = 10 (panel B1) and 𝑠 = 20 (panel B2). The vertical dashed line indicates the relative error result for the flow in the panel above. 100 10 14.65% 14.02% 38.50% 12.05% 16.40% 28.73% 74.64% 91.52% 91.52% 91.52% 91.52% 91.52% 91.52% 91.52% 20 8.01% 8.11% 21.10% 3.03% 10.17% 4.22% 5.89% 9.13% 16.27%100.00%100.00%100.00%100.00%100.00% 1 10 relative error % 30 4.57% 5.36% 10.49% 5.82% 4.02% 3.47% 9.24% 4.55% 6.38% 7.08% 9.77% 96.05% 72.59% 63.96% NOBS N N40 2.60% 3.78% 5.89% 5.29% 1.57% 3.32% 0.12% 1.31% 0.91% 2.09% 1.62% 4.72% 0.53% 0.38% 2 10 50 1.39% 2.77% 3.59% 4.07% 0.40% 2.39% 0.66% 1.02% 1.42% 1.88% 1.55% 0.39% 5.75% 4.03% 60 0.63% 2.07% 2.34% 2.94% 0.27% 1.59% 0.64% 0.61% 1.44% 1.00% 2.46% 0.28% 1.90% 2.01% 3 70 0.13% 1.58% 1.62% 2.11% 0.63% 0.94% 0.61% 0.53% 0.46% 0.12% 0.05% 1.03% 0.70% 0.97% 10 2 3 4 5 6 7 8 9 10 11 12 13 14 15 r Figure 7.5 DMD emulation of 16 𝑂, details as in Figure 7.4. Relative error in the DMD-predicted evolved energy as a function of number of observations and truncation rank 𝑟. 96 0.08 0.025 HDMD HIMSRG / HIMSRG 0.020 0.06 0.015 0.04 0.010 0.02 0.005 O16 Chiral NN N3LO Ne20 Chiral NN N3LO 0.000 O16 AV18 0.00 Ne20 AV18 0 100 200 300 400 0 200 400 600 800 1000 Number of iterations Number of iterations Figure 7.6 DMD emulation of the IMSRG(2) flow for 16 O 20 Ne in two different interactions AV18 and a chiral NN interaction. The right pane is eMax = 6 basis size, while the right pane is eMax = 10. The vertical axis plots the Frobenius error relative to the IMSRG(2) result. flow in different interactions and different nuclei, 16 O and 20 Ne. Note that 20 Ne open-shell, which involves stiffer IMSRG(2) flow equations because the reference state is harder to decouple. Also note that AV18 [14] is a very “hard” phenomenological interaction, compared to the chiral interactions in Chapter 4.2. We show the DMD is able to reproduce the full flowing Hamiltonian with high accuracy in every case, meaning the DMD is able to fully capture the dynamics of the problem regardless of the interaction. The DMD emulates the results within the finite subset of the operator algebra that characterizes the IMSRG flow. 7.5 Emulating the Parametric IMSRG Flow The “parametric” IMSRG is an extension to the IMSRG problem, which includes the cou- pling constants that parameterize the evolving the Hamiltonian. For the purposes of uncertainty quantification and sensitivity analyses, we need to explore the dependence of the IMSRG results on the parameters of the chiral interactions (Chapter 4.2). The parameter space is very large at 15+ dimensions and we require millions, or even billions, of samples to characterize the space. Collecting this volume of samples is not feasible with direct IMSRG calculations. Thus, we want to leverage DMD reduced-order models of the IMSRG flow to emulate the results. The super operator 97 we wish to emulate becomes 𝐻 (𝑐 1 , . . . , 𝑐 𝑁𝑐 ; 𝑠) = U (𝑐 1 , . . . , 𝑐 𝑁𝑐 ; 𝑠)𝐻 (𝑐 1 , . . . , 𝑐 𝑁𝑐 ; 0) (7.22) where 𝑐𝑖 represents a particular coupling constant (e.g. the parameters of the P3H Hamiltonian, or the LECs of the chiral interactions). The coupling constants which characterize the evolving Hamiltonian also characterize the unitary transformation. In this emulation scheme, we apply the DMD to emulate the unitary transformation as a trajectory 𝑠 along the coupling manifold characterized by 𝑐 1 , . . . , 𝑐 𝑁𝑐 . The primary algorithm we employ to tackle the parametric problem is Reduced Koopman Operator Interpolation (rKOI), introduced by Huhn et al. [6]. In the rKOI scheme, we directly interpolate the DMD operator by sampling a training set of DMD operators computed from varying IMSRG trajectories along the Hamiltonian coupling manifold. Algorithm 7.3 provides an overview of the rKOI algorithm. The key feature of this strategy is that the reduced DMD operator is interpolated directly. The DMD expansion computed from Algorithm 7.3 Pseudocode implementation which demonstrates the algorithm for rKOI DMD. vary coupling constants 𝑐 1 , . . . , 𝑐 𝑁𝑐 𝑁train times and compute IMSRG trajectory 𝐷 𝑖 for each realization for each 𝐷 𝑖 in 𝑁train size set do compute projected DMD in Algorithm 7.2 record 𝐴𝑟(𝑖) , 𝑈𝑟(𝑖) , 𝑏 (𝑖) end for train interpolators to predict I𝐴𝑟 , I𝑈𝑟 , I𝑏 on training data {𝐴𝑟(𝑖) }, {𝑈𝑟(𝑖) }, {𝑏 (𝑖) } compute DMD expansion in Algorithm 7.2 from interpolated operators interpolated operators reads 𝑥 I𝑘 = 𝑈𝑟I 𝑊 I (ΛI ) 𝑘−1 𝑏 I (7.23) where we explicitly include how the interpolated operators enter into the expansion. Note that 𝑊 I , ΛI are computed by the eigendecomposition of 𝐴𝑟I . We refer to Algorithm 7.3, or the process of ingesting training data to compute the interpolated operators, as the interpolation engine (IE). The computed operators from which we may construct the DMD expansion, and therefore predict the evolution of the system given a target set of coupling 98 constants, we refer to as the interpolation engine model, or IEM. The trained IEM fully captures the transformation to reduced space 𝑈𝑟 , the dynamical trajectory through the coupling constant manifold encoded by 𝐴𝑟 , and the space of initial conditions—this means that a prediction using this interpolation scheme is data-free, in the sense that the only input necessary to predict every feature of the dynamical system is the target set of coupling constants. Of course, the accuracy of the prediction depends on the success of the interpolation method, which itself depends on the well-behaved nature of the manifold we wish to interpolate. 7.6 Interpolating DMD using RBFs We have chosen to interpolate the required DMD operators using a radial basis function (RBF) interpolator1, in linear and thin plate spline kernels. The mesh-free nature of RBF interpolation was the deciding factor for its implementation in DMD interpolation—the RBF-IEM can freely provide predictions not constrained to a mesh. Instead of a mesh, we employ Latin hypercube sampling for optimal coverage of the parameter space, and the nature of the RBF interpolation allows us to generalize the model inside and outside of the sampled region. In the following, we demonstrate the success of the RBF-IEM in one and two dimensions, corresponding to the interaction parameters 𝑔, 𝑏 in the P3H model (Chapter 4.1), as well as the success of interpolation in a realistic nuclear system. Interpolating in 1D Toy Model Figure 7.7 illustrates how the RBF-IEM reproduces the evolved energy in the P3H model across a wide range of pairing strengths 𝑔, with pair-breaking strength 𝑏 fixed at zero. In the 1D case, the training points have been sampled uniformly except at 𝑔 = 0 (we exclude 𝑔 = 0 from the training space because the energy is trivially zero). The predicted energy matches the IMSRG energy with reasonable accuracy, at least on the order 10−3 , with a notable loss in accuracy beyond the training point boundaries—this is a disadvantage of the mesh-free interpolation. 1 SciPy documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RBFInterpolator. html#scipy.interpolate.RBFInterpolator. 99 Training points IMSRG(2) evolved 2.5 interpolated DMD reconstructed 2.0 E evolved 1.5 1.0 interpolated NOBS: 50 TRUNC: 3 0.5 EMU: rKOI 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 g Figure 7.7 IMSRG(2) ground-state energy as a function of the pairing strength 𝑔, comparing the interpolated DMD prediction against the IMSRG(2)-evolved result. The red dots indicate IMSRG(2) evolution training points for the interpolation engine. Each training point represents 𝑁 = 50 IMSRG(2) iterations which are used to construct a training DMD operator at SVD truncation level 𝑟 = 3. Interpolating in 2D Toy Model In Figure 7.8, we now vary 𝑔 and 𝑏 in the P3H model. The training points have been collected using a Latin hypercube sampling strategy. The Latin hypercube sampling promotes optimal coverage by dividing the parameter space into equal partitions, and then sampling one point from each partition. Note that the error range is fairly wide, with most predictions received at an error 10−1 to 10−4 . The distribution of training points seems to affect the locality of error; that is, regions which are less dense with training points impute a larger error in those regions. At the same time, there are particular regions, such as the lower left quadrant, where increasing the number of training points does not seem to affect the error significantly. This is likely improvement saturation constrained by the reconstruction error of the individual DMD 100 |(EDMD - EIMSRG)/EIMSRG| (NTRAIN = 50, RMS error = 0.0331) |(EDMD - EIMSRG)/EIMSRG| (NTRAIN = 300, RMS error = 0.0211) -1.0 -1.0 -0.9 -0.9 -0.8 -0.8 10 2 -0.7 10 2 -0.7 -0.6 -0.6 -0.5 -0.5 -0.4 -0.4 -0.3 -0.3 10 3 -0.2 -0.2 -0.1 10 3 -0.1 g 0.0 g 0.0 0.1 0.1 0.2 0.2 0.3 0.3 10 4 0.4 0.4 0.5 10 4 0.5 0.6 0.6 0.7 0.7 0.8 0.8 10 5 0.9 0.9 1.0 1.0 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 b b 0.8 0.9 1.0 0.8 0.9 1.0 Figure 7.8 Absolute relative error in the evolved energy, comparing DMD emulation to IMSRG(2) evolution. Each training point represents a rank 3 DMD operator computed from 50 observations. (LEFT) 50 training points. (RIGHT) 300 training points. training points. In other words, the parametric emulator can only predict as well as the DMD operators on which it has been trained. In Figure 7.9, we plot the absolute relative error of direct DMD forecasts on the same grid points as Figure 7.8. The results of DMD forecasting reveals that the RBF-IEM error is constrained by the accuracy of the DMD itself. We might improve the error of the DMD forecast by increasing the truncation rank. However, a larger DMD operator results in a more complicated parameter space to interpolate. Interpolating a Realistic Nucleus Ultimately, our goal is to construct a practical IMSRG(2) emulator for realistic (nuclear) many-body calculations. There is an array of additional challenges associated with the IE applied in this case, including memory requirements (compare the P3H Hamiltonian of 103 IMSRG(2) coefficients to a realistic nucleus which requires 106 −109 IMSRG(2) coefficients), but we will demonstrate that they can be overcome. As an exampple, we will construct an RBF-IEM to model realistic IMSRG(2) calculations of the closed-shell 16 O nucleus, varying 24 low energy constants (LECs) in the underlying chiral N3 LO NN interaction by Entem and Machleidt 101 |(EDMD - EIMSRG)/EIMSRG| (RMS error = 0.0506) -1.0 0.035 -0.9 -0.8 -0.7 0.030 -0.6 -0.5 -0.4 0.025 -0.3 -0.2 -0.1 0.020 g 0.0 0.1 0.2 0.015 0.3 0.4 0.5 0.010 0.6 0.7 0.8 0.005 0.9 1.0 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 b 0.8 0.9 1.0 Figure 7.9 Absolute relative error in the evolved energy for DMD emulation, measured by IM- SRG(2). Each grid point is the evolved energy predicted by a rank 3 truncated DMD operator directly computed from the flow, for the corresponding grid point. (Chapter 4.2). Figure 7.10 summarizes the interpolation error on 16 O. The RBF-IEM was trained on 300 points, where each training point is a rank 8 truncated DMD operator computed from 50 observations of the evolving 16 O Hamiltonian. We plot two error measurements, relative error in the evolved energy and flow reconstruction error (mean-squared error of discrete flow matrix reproduced by DMD). The relative error in the energy measures how well the RBF-IEM predicts the target evolved energy, whereas the reconstruction error measures how well the RBF-IEM models the full dynamics of the associated IMSRG(2) flow. Note the variability is rather high in both error measurements. We measure the error according to direct IMSRG(2) calculations using parameters which were not in the training set. Validation points with particularly high error likely are “far” from the training 102 |(EDMD - EIMSRG)/EIMSRG| Flow reconstruction error (MSE) 100 100 abs_rel_E_error 10 1 mse 10 1 10 2 0 10 20 30 40 50 0 10 20 30 40 50 Validation point (index) Validation point (index) Figure 7.10 Error measurements from interpolation engine trained on 300 training points. In both panels, the error measurement is distributed across 50 random validation points. The left pane plots absolute relative error in the evolved energy. The right pane plots mean-squared error in the reconstructed flow from the interpolated DMD operators. points, meaning that the model does not seem to generalize well, and that more training points could tighten the error spread. Figure 7.11 provides some evidence for this claim. Note that variations in the reconstruction error seem to correlate, although not perfectly, with variations in the average distance from each validation point to the training set. Additionally, interpolating the parameter space of the operator 𝑈𝑟 , which transforms the DMD operator to reduced space, is difficult since the size of the operator is equal to the number of state variables in the dynamical system. This is in contrast to the parametric space of the DMD operator itself, which is 82 parameters in this case. To tackle this issue, we modify our approach by performing a reduction on the training space of IMSRG(2) flows, taking inspiration from the partitioned approach to parametric DMD introduced by Andreuzzi, Demo, and Rozza [1]. Instead of interpolating the full dimension of the training space, we reduce the feature dimension of the training space to the principal components with the help of an SVD. 103 1.0 0.8 0.6 0.4 0.2 Average distance to training point Reconstruction error 0.0 0 10 20 30 40 50 Validation index Figure 7.11 Illustrated in this plot is the reconstruction error (mean-squared error) in Figure 7.10, as well as the validation point average distance to the training points. The vertical axis for each measurement has been scaled so that the data appear between 0 and 1, so that variations can be compared. Borrowing the notation of Andreuzzi, Demo, and Rozza [1], let    | ... |    h i 𝒄𝑖   𝑋 B 𝑥  𝒄 𝑖 𝒄 𝑖  ∈ R𝑚×𝑁 , X1 B 𝑋 𝒄1 . . . 𝑋 𝒄 𝑝 ∈ R𝑚×𝑁 𝑝 , (7.24) . . . 𝑥 𝑁  1   | ... |      𝒄 where 𝑋 𝑖 represents 𝑁 observations of the IMSRG(2) flow according to parametric set 𝒄𝑖 (e.g., a particular set of LECs). There is evidence of low rank structure in the IMSRG(2) flow (we thank Boyao Zhu and Heiko Hergert for their fruitful discussions in this regard; in addition, see [5]), hence we expect that the training matrix of IMSRG(2) flows given by X1 will also exhibit low rank structure. Indeed, Figure 7.12 shows that its singular value spectrum is decaying exponentially. We proceed with SVD on X1 , truncating at 𝑡 singular values. We emphasize that 𝑡 denotes the reduction of the feature space associated with the training data, while 𝑟 denotes the truncation of individual reduced DMD operators. Note that in the IMSRG context 𝑡 is likely to be much greater than the truncation rank 𝑟 of the reduced DMD operator. The truncated SVD of X1 reads X1 ≈ L𝑡 𝚺𝑡 R𝑡∗ . (7.25) 104 10 1 10 3 10 5 10 7 si/s1 10 9 10 11 10 13 10 15 0 2000 4000 6000 8000 10000 12000 14000 Singular value index i Figure 7.12 Singular value spectrum of the training matrix X1 , scaled by the largest singular value. Finally, the state variable information in X1 is compressed using L𝑡 , X̃1 = L𝑡∗ X1 ∈ R𝑡×𝑁 𝑝 , (7.26) and parametric training may proceed with Algorithm 7.3 using the reduced training matrix. We refer to the model produced by the reduction on the training space as reduced IEM, or rIEM. RBF-rIEM is the RBF-interpolated version of the model. We might notate the truncation rank of the reduced training space directly in the acronym for clarity, e.g. RBF-r(𝑡)IEM In Figure 7.13, we have trained the RBF-r(1000)IEM on 300 rank-2 reduced DMD operators in the reduced training space corresponding to truncation 𝑡 = 1000. Comparing with Figure 7.10, we observe smaller variations in both measurements, with fluctuations staying between the orders 10−2 and 10−1 . We have show that the modified parametric DMD algorithm for training space reduction in- creased prediction accuracy across a wide range of validation points. As a result of the training space reduction, the time to compute the DMD operator training points, as well as the time to train 105 |(EDMD - EIMSRG)/EIMSRG| Flow reconstruction error (MSE) 10 1 10 1 abs_rel_E_error mse 10 2 10 2 0 10 20 30 40 0 10 20 30 40 Validation point (index) Validation point (index) Figure 7.13 Error measurements from interpolation engine trained on 300 reduced training points. In both panels, the error measurement is distributed across 50 random validation points. (LEFT) Absolute relative error in the evolved energy. (RIGHT) Mean-squared error in the reconstructed flow from the interpolated DMD operators. the IEM, was also reduced. 7.7 Optimized Training for Large Bases using Streaming SVD The reduced IE method introduces the new cost of computing the SVD on the training space in Equation 7.25. More specifically, we need the left singular vectors L𝑟 in order to perform the reduction on the training space. Thus, the training algorithm performance hinges on the computation of a very “tall and skinny” matrix X1 , which likely will not completely fit into memory. To address this problem, we employ the streaming SVD, or Sequential Karhunen-Loeve (SKL) algorithm introduced by Levey and Lindenbaum [9], for computing the left singular vectors 𝑈 of the full SVD via looping over data rectangles, stored on disk, for which each iteration R-SVD is performed to extract the major singular vectors from a QR-decomposition of the rectangle. We show the method’s pseudocode in Algorithm 7.4 and provide a schematic diagram in Figure 7.14. The forget factor ff is a fraction between 0 and 1 that controls the contribution of each previous set of singular vectors to the next. Algorithm 7.4 partitions calculation of 𝑈 such that each partition of columns of the full matrix can be introduced into memory sequentially. The complexity of computing the left singular values this way is comparable to the complexity of the full SVD, except 106 Algorithm 7.4 Pseudocode implementation which demonstrates the algorithm for streaming SVD, originally conceived by Levey and Lindenbaum [9]. Algorithm adapted from Maulik and Mengaldo [10]. Initialization (rectangle 𝑘 = 0): Compute QR-decomposition 𝐴0 = 𝑄 0 𝑅0 Compute SVD 𝑅0 = 𝑈0′ 𝐷 (𝑘=0) 𝑉0𝑇 (truncating at desired rank) Compute 𝑈 (𝑘=0) = 𝑄 0𝑈0′ (where 𝑈 𝑘 corresponds to the full SVD modes at update 𝑘) for each rectangle 𝑘 >= 1 and matrix 𝐴 𝑘 do Append new data 𝐴 𝑘 to previous basis (ff · 𝑈 𝑘−1 𝐷 𝑘−1 | 𝐴 𝑘 ) = 𝐴′𝑘 Compute QR-decomposition 𝐴′𝑘 = 𝑄 𝑘 𝑅 𝑘 Compute SVD 𝑅 𝑘 = 𝑈 𝑘′ 𝐷 (𝑘=0) 𝑉𝑘∗ (truncating at desired rank) Update singular vectors 𝑈 (𝑘) = 𝑄 𝑘 𝑈 𝑘′ end for Off-line: Data partitions, U, and D, are k= 0 … k= j-1 k= j … k= P stored on disk On-line: U and D built from 𝐷(𝑗−1) previous k=j-1 (𝑗−1) 𝑅𝑗 rectangles are updated 𝑈 + 𝐴𝑗 = 𝑄𝑗 from k=j rectangle, via local R-SVD Update U 𝑅𝑗 = local R-SVD 𝑈′𝑗 and truncate ∗ 𝑈 (𝑗) = 𝑄𝑗 at rank T 𝑈′𝑗 𝐷𝑗 𝑉𝑗′ Figure 7.14 Schematic representation of the SKL algorithm [9]. In the on-line stage, the SVD modes and singular values are updated are updated according to the local data at 𝑘 = 𝑗. The updated modes and singular values can be written to disk and unloaded from memory until the next update. that the memory constraint has been relaxed. Note that the modes produced by the SKL algorithm are robust to noise introduced by concatenation of new data [9]. In general, one does not worry about significant singular vectors being lost in the concatenation. This fact has been tested and proven by Levey and Lindenbaum [9]. An intuitive interpretation of this fact may be that for a “good” random initialization point of the stream, stochastic introduction of new data tends to 107 300 rectangles streamed 104 250 rectangles streamed 200 rectangles streamed 150 rectangles streamed 102 100 rectangles streamed si/s1 50 rectangles streamed full SVD 100 10 2 0 200 400 600 800 1000 Singular value index i Figure 7.15 The singular value spectrum of 300 training points corresponding to 16 𝑂, in an eMax = 6 harmonic oscillator basis size and oscillator width ℏ𝜔 = 20. suppress modes which are significantly different from previous modes, in an orthogonality sense. Additionally, the rate of convergence of the streaming modes depends on the underlying density of the low rank structure in the streaming data; that is, streaming data which exhibits a quickly decaying singular value spectrum will converge faster under the SKL algorithm than streaming data for which the spectrum decays more slowly. Figure 7.15 displays the singular value spectrum convergence of the streaming SVD as new rectangles are incorporated. Note the dashed line which displays the spectrum of the equivalent full SVD, i.e. performing SVD on the full dataset at once. The largest 25-50 singular values converge the fastest in this truncation scheme. This means that, depending on the desired truncation rank in the streaming SVD, we may not need to stream every available data point to arrive at a satisfactory low rank approximation of the dataset. 7.7.1 A Parallel Implementation of the Streaming SVD While current parallel implementations of streaming SVD focus on parallelization of the QR algorithm itself [2, 10], we present a different approach that focuses on increasing the throughput of the SKL algorithm. We simply call this approach parallel-SKL, or PSKL. We start with the fact that the SKL algorithm makes no assumptions on the ordering of the data partitions stored on disk, 108 P=0 P=1 P=2 P=3 Each task 𝑃=0 𝑃=1 𝑃=2 𝑃=3 completes 𝐷𝑗=𝑘 𝑃=0 𝐷𝑗=𝑘 𝑃=1 𝐷𝑗=𝑘 𝑃=2 𝐷𝑗=𝑘 𝑃=3 𝑃=0 𝑃=1 𝑃=2 𝑃=3 PSKL 𝑈𝑗=𝑘 𝑃=0 𝑈𝑗=𝑘 𝑃=1 𝑈𝑗=𝑘 𝑃=2 𝑈𝑗=𝑘 𝑃=3 algorithm for assigned partition Update final SVD modes 𝐷(𝑃=0) Streaming × 𝐷𝑃=1 Streaming × 𝐷𝑃=2 Streaming × 𝐷𝑃=3 𝑈 (𝑃=0) update with 𝑈𝑃=1 update with 𝑈𝑃=2 update with from each 𝑈𝑃=3 P=1 task P=2 task P=3 task task into each other Figure 7.16 Diagram of a 4-task parallel-SKL algorithm. Each task 𝑃 is assigned a partition of data indices, which together span the total range of the target data. Then, each task 𝑃 computes its own streaming SVD of its data partition. Finally, the modes computed from each task are updated into each other, sequentially. nor does the on-line update of SVD modes in Figure 7.14 account for which data partitions are used, or when they are ingested. However, the SVD modes at update step 𝑘 do contain information about all data partitions before step 𝑘 (to the extent controlled by the forget-factor ff), in the form of iterative projections of old information onto new information. Therefore, the parallel structure we present is to simply partition the ingestion of off-line data over several, parallel tasks, each of which perform their own streaming SVD of the data they are assigned. Once each task is completed, these separate, parallel SVD modes may be sequentially updated into each other, generating a single set of SVD modes which contain information from all tasks. Figure 7.16 displays a diagram which illustrates the PSKL algorithm. Note that the information required for sequential updating after each task is completed is stored in the matrix multiplication of U𝑃 D𝑃 . Figure 7.17 confirms that information from each separate task is preserved during the final combination step. The PSKL singular value spectrum, computed using 2 tasks in this example, is consistent with the spectrum of a full SVD calculation. 109 105 full SVD 2 task PSKL 103 si/s1 101 10 1 0 200 400 600 800 1000 Singular value index i Figure 7.17 Singular value spectrum comparison between the full SVD and a 2-task PSKL run, for a training matrix of 300 points corresponding to 16 𝑂, in an eMax = 6 harmonic oscillator basis size and oscillator width ℏ𝜔 = 20. Assuming each task has access to resources identical to the sequential version, then the PSKL algorithm performance scales at least linear in the number of tasks. Of course, this resource dupli- cation across tasks also implies linear scaling in the total memory consumption of the algorithm. However, the performance gains may be worth the memory costs on systems where large amounts of memory are available. 7.8 Outlook In this chapter, we have presented a method for constructing an IMSRG(2) emulator, using DMD as the driving mechanism for predicting the IMSRG dynamics. We have shown that DMD may be used to forecast a single IMSRG(2) evolution (robust to different interactions and nuclei), reducing the number of polynomial-cost ODE solver iterations required for a converged result. We have also shown that DMD may be used to construct an interpolation engine for the purpose of training a model that emulates the flow for any set of parameters that characterize the Hamiltonian including resultion scales and cutoffs. In this case, we desire a rich set of IMSRG(2) evolutions in order to populate the IEM training set with high-fidelity DMD operators. The larger the training 110 set, the more difficult data ingestion becomes; we have presented a solution to this problem using PCA to reduce the training set, in tandem with the SKL algorithm for streaming SVD. With all these improvements combined, we are able to construct a confident model of the IMSRG(2) evolutions of a target nucleus and the associated Hamiltonian, parameterized by the LECs in chiral EFT. In Chapter 8, we present a practical application of RBF-rIEM for IMSRG(2) applied to nuclear Hamiltonians: namely, global sensitivity analysis of energies to variations of the LECs. This analysis is important to determining the contribution of the variance of particular LECs to the total variance in the evolved energy. In future work, this sensitivity knowledge may let us generate more efficient training bases for higher fidelity IEMs, as well as construct efficient sampling distributions for a proper statistical uncertainty quantification of observables that are computed with IMRSG methods. 111 BIBLIOGRAPHY [1] Francesco Andreuzzi, Nicola Demo, and Gianluigi Rozza. A dynamic mode decomposition extension for the forecasting of parametric dynamical systems. 2021. doi: 10.48550/ ARXIV.2110.09155. url: https://arxiv.org/abs/2110.09155. [2] Austin R. Benson, David F. Gleich, and James Demmel. “Direct QR factorizations for tall- and-skinny matrices in MapReduce architectures”. In: 2013 IEEE International Conference on Big Data. 2013, pp. 264–272. doi: 10.1109/BigData.2013.6691583. [3] Steven L. Brunton and J. Nathan Kutz. Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control. Cambridge University Press, 2019. doi: 10.1017/9781108380690. [4] Steven L. Brunton et al. Modern Koopman Theory for Dynamical Systems. 2021. arXiv: 2102.12086 [math.DS]. [5] J. Hoppe et al. “Importance truncation for the in-medium similarity renormalization group”. In: Phys. Rev. C 105 (3 Mar. 2022), p. 034324. doi: 10.1103/PhysRevC.105.034324. url: https://link.aps.org/doi/10.1103/PhysRevC.105.034324. [6] Quincy A. Huhn et al. Parametric Dynamic Mode Decomposition for Reduced Order Mod- eling. 2022. doi: 10.48550/ARXIV.2204.12006. url: https://arxiv.org/abs/2204.12006. [7] B. O. Koopman. “Hamiltonian Systems and Transformation in Hilbert Space”. In: Proceedings of the National Academy of Sciences 17.5 (1931), pp. 315–318. doi: 10.1073/pnas.17.5.315. eprint: https://www.pnas.org/doi/pdf/10.1073/pnas.17.5.315. url: https://www.pnas.org/doi/abs/10.1073/pnas.17.5.315. [8] J. Nathan Kutz et al. Dynamic Mode Decomposition. Philadelphia, PA: Society for Industrial and Applied Mathematics, 2016. doi: 10.1137/1.9781611974508. eprint: https://epubs. siam.org/doi/pdf/10.1137/1.9781611974508. url: https://epubs.siam.org/doi/abs/10.1137/ 1.9781611974508. [9] A. Levey and M. Lindenbaum. “Sequential Karhunen-Loeve basis extraction and its appli- cation to images”. In: IEEE Transactions on Image Processing 9.8 (2000), pp. 1371–1374. doi: 10.1109/83.855432. [10] Romit Maulik and Gianmarco Mengaldo. PyParSVD: A streaming, distributed and random- ized singular-value-decomposition library. 2021. arXiv: 2108.08845 [cs.MS]. [11] Peter Schmid and Joern Sesterhenn. “Dynamic Mode Decomposition of numerical and experimental data”. In: APS Division of Fluid Dynamics Meeting Abstracts. Vol. 61. APS Meeting Abstracts. Nov. 2008, MR.007, MR.007. 112 [12] Peter J. Schmid. “Dynamic mode decomposition of numerical and experimental data”. In: Journal of Fluid Mechanics 656 (2010), pp. 5–28. doi: 10.1017/S0022112010001217. [13] Jonathan H. Tu et al. “On dynamic mode decomposition: Theory and applications”. In: Journal of Computational Dynamics 1.2 (2014), pp. 391–421. issn: 2158-2491. doi: 10.3934/jcd.2014.1.391. url: /article/id/1dfebc20-876d-4da7-8034-7cd3c7ae1161. [14] R. B. Wiringa, V. G. J. Stoks, and R. Schiavilla. “Accurate nucleon-nucleon potential with charge-independence breaking”. In: Phys. Rev. C 51 (1 Jan. 1995), pp. 38–51. doi: 10.1103/PhysRevC.51.38. url: https://link.aps.org/doi/10.1103/PhysRevC.51.38. 113 CHAPTER 8 GLOBAL SENSITIVITY ANALYSIS USING EMULATION OF THE IMSRG 8.1 Introduction In this chapter, we present a practical application of the parametric IMSRG emulator: A global sensitivity analysis (GSA) of energies to variations in the low-energy constants (LECs) which define the nuclear interactions constructed from Chiral EFT, as explained in Chapter 4.2. We present several GSAs for selected nuclei at interaction chiral orders N2 LO (including NN and 3N interactions) and N3 LO (using NN forces only). Our strategy for GSA mirrors the strategy of Ekström and Hagen [3], so that our results for 16 O with the NNLOsat interaction should be directly comparable. Note that this GSA would not be feasible without the improvements achieved by the reduced-basis interpolation strategy of Chapter 7.6. The sensitivity analysis was performed with aid from the SALib1 Python package [8, 6]. We follow the Sobol’ approach to global sensitivity analysis [10, 9]. The sensitivity measures the variance in the model output captured by a change in the model parameters. To begin, the total variance 𝐷 of the model output 𝑌 , with respect to the model parameters, is decomposed into a sum of partial variances relative to simultaneous changes in an increasing number of parameter 𝑐𝑖 . Denoting the variance for varying a single parameter while keeping all others fixed 𝐷 𝑖 , for two parameters 𝐷 𝑖 𝑗 , etc., the decomposition reads ∑︁ ∑︁ 𝐷= 𝐷𝑖 + 𝐷 𝑖, 𝑗 + · · · . (8.1) 𝑖 𝑖, 𝑗,𝑖≠ 𝑗 The first-order sensitivity to parameter 𝑐𝑖 is given by 𝐷𝑖 𝑆𝑖 = , (8.2) 𝐷 and the total sensitivity to the parameter 𝑐𝑖 is given by a sum over all orders of sensitivity, 𝑆𝑇,𝑖 = 𝑆𝑖 + 𝑆𝑖 𝑗 + · · · 𝑆𝑖...𝑁 . (8.3) 1 Documentation can be found at https://salib.readthedocs.io/en/latest/index.html. 114 The first-order sensitivity integral with respect to parameter 𝑐𝑖 is [9] Var𝑖 [E∼i [𝑌 | 𝑐𝑖 ]] 𝑆𝑖 = , (8.4) 𝐷 where the expectation value 𝐸 index ∼ i represents the set of all parameters excluding 𝑐𝑖 . The total sensitivity integral with respect to the parameter 𝑐𝑖 is [9] Var∼i [E𝑖 [𝑌 | 𝑐 ∼i ]] 𝑆𝑇,𝑖 = 1 − . (8.5) 𝐷 𝑆𝑖 describes the sensitivity of the model output 𝑌 to 𝑐𝑖 while fixing all other parameters, i.e. 𝑐 ∼i . 𝑆𝑇,𝑖 describes the sensitivity of 𝑌 to 𝑐𝑖 including the first-order effect plus all higher orders. Thus, the total sensitivity encodes correlations between parameters in the model output. In the Sobol’ approach to global sensitivity analysis, the integrals in Equations 8.4 and 8.5 are estimated via Sobol’ sequence sampling of the model output, which is used for quasi-MC methods. The Sobol’ sequence successively divides the parameter space in half, for an optimal number of partitions, such that random samples in each partition optimally cover the full parameter space. Note that sampling the RBF-rIEM can be very fast, provided the reduced model is being sampled, i.e. the model before transformation back to the full IMSRG space via L𝑟 (see Equation 7.26). Once all relevant samples have been collected, the reduced samples may be transformed back to the full space in a single batch by applying the learned transformation. 8.2 GSA with a Chiral N3LO Nucleon-Nucleon Interaction In Figure 8.1, we present a Sobol’ sensitivity analysis of the 16 O ground-state energy from IMSRG(2) calculations with chiral N3LO NN interactions with cutoff Λ𝑏 = 500 MeV. Our baseline is the NN3 LO(500) interaction by Entem and Machleidt [4], for nucleus 16 O. We trained a RBF- r(50)IEM of 16 O on 300 DMD operator training points built from 50 observations of corresponding IMSRG(2) evolutions. We used an eMax = 12 harmonic oscillator basis size with oscillator width ℏ𝜔 = 20 Mev. The corresponding LEC training set has been sampled from a Latin hypercube with ±20% relative variations around the standard N3LO LEC set in Table 4.2. The particular contact LECs which were varied appear on the horizontal axis of 8.1. The analysis was performed with 115 60 O16 60000 1 : -55.89±14.81 MeV Energy first-order sensitivity S1 Sensitivity % 40 Energy total sensitivity ST 40000 20 Count 20000 0 0 100 C1 np C33D1 D1 3 D3 P2 D3 75 50 25 0 S C3 0np S C11 C3 S C3 0 P C1 0 P1 P C3 1 S1 D1 2 D1 D3 D3 S D31 D3 1 P S0 S D3 0 P D1 0 P1 P1 S S1 D1 D13D1 D D3 2 D33F2 S1 D3 D2 P2 C3 S D31 D3 Figure 8.1 First-order and total sensitivity of the IMSRG(2) ground-state energy for 16 O energy to variation of the LECs around the standard N3 LO(500) interaction of Entem and Machleidt [4]. The left panel contains the sensitivity information per LEC. The right panel plots the total variance in the energy. The shaded region in the variance plot represents the 1𝜎 and 2𝜎 range. 13,631,488 Sobol’ sequence samples of the RBF-r(50)IEM, totaling about 10 minutes of wall- time on a high performance computer, in Figure 7.13. The same number of complete IMSRG(2) calculations, on the same high performance computer, would require approximately 7700 years of compute time. We sample within a ±10% relative variation of the standard N3 LO LECs. The two most sensitive parameters in 𝐶˜3𝑆1 and 𝐶1𝑆0 are consistent with the N2 LO sensitivities presented by Ekström and Hagen [3]. The first-order sensitivity matches almost exactly to the total sensitivity in all LECs, meaning that variations in single LECs do not affect the others, at least according to the fitted RBF-r(50)IEM. The top four most sensitive parameters, in this model, are 𝐶 e3𝑆1 , 𝐶1𝑆0 , 𝐶3𝑆1 −3𝐷 1 , and 𝐶3𝑃2 . Note e𝑛𝑝 . 𝐶1𝑆0 is the subleading contact that we neglect the isospin breaking in the leading 1 𝑆0 contact 𝐶1𝑆0 in the 1 𝑆0 partial wave, so it probably renormalizes the core. 𝐶 e3𝑆1 and 𝐶3𝑆1 −3𝐷 1 are in the deuteron partial waves — the mixed waves is impacted by the strength of the tensor interaction. The S=1 P-waves (𝐶3𝑃 𝑥 ), and 𝐶3𝑃2 in particular, probably renormalize short-range spin-orbit physics. 8.3 GSA with Chiral NNLO Two- Plus Three-Nucleon Interactions In this section, we present a Sobol’ sensitivity analysis of the Chiral NNLO LECs by Ekström et al. [2, 1], up to the 3B forces. Here we have performed the analysis using interactions constructed from “delta-full” and “delta-less” chiral EFT. In this “delta-full” ΔNNLOGO (394) interaction we explicitly add fixed resonance saturation values to the 𝑐 3 and 𝑐 4 LECs in the Fujita-Miyazawa force, because the low cutoff allows us to absorb into the adjust of 𝑐 3 and 𝑐 4 , without the need to 116 implement the long-range term explicitly. The LECs 𝑐 3 , 𝑐 4 are adjusted by an amount 𝑐Δ3 , 𝑐Δ4 [1], 4ℎ2𝐴 𝑐Δ3 = −2𝑐Δ4 = = −2.972246 GeV−1 , (8.6) 9𝛿 for final values that read, Δ 𝑐3𝑁 3 = 𝑐 3 + 𝑐 3 = −3.61 (8.7) Δ 𝑐3𝑁 4 = 𝑐 4 + 𝑐 4 = +2.44 Additionally, we must include the 𝑐 2 𝜋𝑁 LEC in the variation, for a total of 19 variation parameters (see Table 4.1, Chapter 4.2). The 𝑐 3 and 𝑐 4 LECs in the “delta-less” NNLOsat (450) interaction, which read [2] 𝑐 3 = −3.93 (8.8) 𝑐 4 = 3.77, take on similar values to the resonance-saturation-adjusted 𝑐3𝑁 3𝑁 3 and 𝑐 4 in the ΔNNLOGO (394) interaction after introduction of the Fujita-Miyazawa term. The NNLOsat (450) does not include the 𝑐 2 LEC, for a total of 18 variation parameters. For both interactions, we have trained RBF-r(50)IEMs for each nucleus that was studied. The models were trained on 130 DMD training operators, which were constructed from 50 observations of each corresponding IMSRG(2) evolution. We used an eMax = 12 harmonic oscillator basis size with oscillator width ℏ𝜔 = 20 Mev. Every LEC training set has been sampled from a Latin hypercube with ±20% relative around the corresponding standard LEC set for the chosen interaction. 8.3.1 Delta-full NNLO The trained RBF-r(50)IEM has been Sobol’-sampled 11,010,048 times within a ±10% relative variation around the corresponding standard LEC in Table 4.1. The total wall-time for each Sobol’ sampling routine is on the order of a few minutes, on a high performance computer. A full IMSRG(2) calculation for the same number of Sobol’ samples, on the same high performance computer, would require approximately 6000 years of compute time. The sensitivity analyses in Figure 8.3 show near consistent parameter sensitivity across the studied neutron-rich Calcium isotopes. The most sensitive parameter is the 𝐶 e(𝑛𝑝) . There seems 3𝑆1 117 Ca48 Energy first-order sensitivity S1 1 : -440.68±150.42 MeV 40 Sensitivity % Energy total sensitivity ST 40000 20 Count 20000 0 0 C1 pp S C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11 S1 P2 C33D1 c1 c2 c3 c4 cD cE 1000 500 0 C1 C3 S C3 0 P0 P1 P C3 1 S1 C3 40 Ca50 Energy first-order sensitivity S1 1 : -442.72±143.99 MeV Sensitivity % Energy total sensitivity ST 40000 20 Count 20000 0 0 1000 C1 pp S C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11 S1 P2 C33D1 c1 c2 c3 c4 cD cE 500 0 C1 C3 S C3 0 P0 P1 P C3 1 S1 C3 40 Ca52 Energy first-order sensitivity S1 1 : -453.51±149.21 MeV Sensitivity % Energy total sensitivity ST 40000 20 Count 20000 0 0 C1 pp S C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11 S1 P2 C33D1 c1 c2 c3 c4 cD cE 1000 500 0 C1 C3 S C3 0 P0 P1 P C3 1 S1 C3 40 Ca54 Energy first-order sensitivity S1 1 : -477.24±147.76 MeV Sensitivity % Energy total sensitivity ST 40000 20 Count 20000 0 0 1000 C1 pp S C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11 S1 P2 C33D1 c1 c2 c3 c4 cD cE 500 0 C1 C3 S C3 0 P0 P1 P C3 1 S1 C3 40 Ca56 Energy first-order sensitivity S1 1 : -456.23±158.41 MeV Sensitivity % Energy total sensitivity ST 40000 20 Count 20000 0 0 C1 pp S C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11 S1 P2 C33D1 c1 c2 c3 c4 cD cE 1000 500 0 C1 C3 S C3 0 P0 P1 P C3 1 S1 C3 Figure 8.2 Summary of ΔNNLOGO (394) LEC sensitivity as sampled by the RBF-r(50)IEMs corresponding to each nucleus in neutron-rich Calcium isotopes, based on samples of interactions around the ΔNNLOGO (394) baseline. The left column plots of the sensitivities of each nucleus per LEC. The right column plots the total variance in the energy as measured by the model for that nucleus. The shaded region in the variance plot represents the 1𝜎 and 2𝜎 range. 118 40 E(Ca50 - Ca48) Energy first-order sensitivity S1 1 : -2.04±17.84 MeV Sensitivity % Energy total sensitivity ST 40000 20 Count 20000 0 0 C1 pp S0 C1 np S0 C1 nn S0 C3 pp S1 C3 np S1 C3 nn S1 C1 S1 C3D1 3 P2 c1 c2 c3 c4 cD cE 0 100 C3 C1 C3 C3 S0 P0 P1 P1 S1 C3 E(Ca52 - Ca50) Energy first-order sensitivity S1 1 : -10.78±7.34 MeV Sensitivity % Energy total sensitivity ST 40000 20 Count 20000 0 0 C1 pp S0 C1 np S0 C1 nn S0 C3 pp S1 C3 np S1 C3 nn S1 C1 S1 C3D1 3 P2 c1 c2 c3 c4 cD cE 50 25 0 25 C3 C1 C3 C3 S0 P0 P1 P1 S1 C3 50 E(Ca54 - Ca52) Energy first-order sensitivity S1 50000 1 : -23.74±34.01 MeV Sensitivity % Energy total sensitivity ST 25 Count 25000 0 C1 pp S0 c1 0 200 0 C1 np S0 C1 nn S0 C3 pp S1 C3 np S1 C3 nn S1 C1 S0 S1 C3D1 3 P2 c2 c3 c4 cD cE C3 C1 C3 C3 P0 P1 P1 S1 C3 50 E(Ca56 - Ca54) Energy first-order sensitivity S1 1 : 21.01±38.25 MeV 50000 Sensitivity % Energy total sensitivity ST 25 Count 25000 0 0200 C1 pp S0 C1 np S0 C1 nn S0 C3 pp S1 C3 np S1 C3 nn S1 C1 S1 C3D1 3 P2 c1 c2 c3 c4 cD cE 0 200 C3 C1 C3 C3 S0 P0 P1 P1 S1 C3 Figure 8.3 Summary of two-neutron separation energies (𝑆2𝑛 ) as sampled by the RBF-r(50)IEMs corresponding to each nucleus in neutron-rich Calcium isotopes, based on samples of interactions around the ΔNNLOGO (394) baseline. The left column plots of the sensitivities of each S2N per LEC. The right column plots the total variance in the energy as measured by the model for that S2N. The shaded region in the variance plot represents the 1𝜎 and 2𝜎 range. to be very little dependence on the 𝑐𝑖 ’s. As in the NN3 LO(500) interaction, the LEC sensitivities are fully first-order in the RBF-r(50)IEM model output. In Figure 8.3 we plot the two-neutron separation 𝑆2𝑛 energies corresponding to the Calcium isotopes presented. We calculate the 𝑆2𝑛 energy by simply subtracting between the sampled model outputs of two nuclei (carefully ensuring that each subtraction is between two outputs on the same input LEC set), e.g. 𝐸 ( 50 Ca −48 Ca) = 𝐸 ( 50 Ca) − 𝐸 ( 48 Ca), 𝐸 ( 52 Ca −50 Ca) = 𝐸 ( 52 Ca) − 𝐸 ( 50 Ca), and so on. The sensitivity patterns change here, because small differences in the separate nuclear model outputs are enhanced by division of the total variance in the subsequent sensitivity analysis. A major differing feature is that there appears to be second- and higher-order sensitivity effects in play, since the first-order sensitivity does not fully account for the total sensitivity. This is expected, because the 𝑆2𝑛 energy 119 50 O16 Sensitivity % Energy first-order sensitivity S1 1 : -123.58±34.21 MeV 10000 Count Energy total sensitivity ST 0 0 C1 pp S C3 0pp S C1 1np S C3 0np S C1 n1n S C3 n0n S C11 S1 C33D1 P2 c1 c3 c4 200 100 0 S C3 0 P C1 0 P C3 1 P C3 1 cE cD S1 C3 Figure 8.4 Summary of NNLOsat (450) LEC sensitivity in 16 O as sampled by the RBF-r(50)IEMs around the NNLOsat (450) baseline. The left panel plots of the sensitivities per LEC. The right panel plots the total variance in the energy as measured by the RBF-IEM. The shaded region in the variance plot represents the 1𝜎 and 2𝜎 range. output depends on the outputs of two independent nuclear models. 8.3.2 Delta-less NNLO Each trained RBF-r(50)IEM has been Sobol’-sampled 1,310,720 times within a ±10% relative variation around the corresponding standard NNLOsat (450) LEC set in Table 4.1. The total wall- time for each Sobol’ sampling routine is on the order of a few minutes, on a high performance computer. A full IMSRG(2) calculation for the same number of Sobol’ samples, on the same high performance computer, would require approximately 750 years of compute time. The sensitivity analysis for 16 O in Figure 8.4 was performed using the same interaction as in Ekström and Hagen [3] (with slightly different variation parameters). We note slight differences between our results and Ekström and Hagen [3], especially in the 𝐶1𝑆0 contact, which may be due to differences between our learned models, as well as parametric model implementations. In this RBF-r(50)IEM model, there appears to be some higher-order effects in the 𝑐 3 and 𝑐 4 contacts which are important. There does not appear to be much dependence on the 𝑐 𝐸 , 𝑐 𝐷 contacts relevant to the 3N force. As in the previous section for ΔNNLOGO (394), in Figure 8.5 we show the sensitivity analyses for several neutron-rich Calcium isotopes as sampled by learned RBF-r(50)IEMs working in the NNLOsat (450) interaction. In contrast to the previous section, we find increased importance, as well as higher-order sensitivities, in the 𝑐 3 , 𝑐 4 contacts. These effects carry into the 𝑆2𝑛 energy 120 40 Ca48 Energy first-order sensitivity S1 1 : -462.33±93.11 MeV Sensitivity % Energy total sensitivity ST 10000 20 Count 5000 0 C1 pp S c1 0 750 500 250 C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11S0 S1 C3D1 3 P2 c3 c4 cD cE C3 C1 C3 C3 P0 P1 P1 S1 C3 Ca50 Energy first-order sensitivity S1 1 : -426.47±103.17 MeV 40 Sensitivity % Energy total sensitivity ST 10000 20 Count 5000 0 0 C1 pp S C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11 S1 C3D1 3 P2 c1 c3 c4 cD cE 750 500 250 C3 C1 C3 C3 S0 P0 P1 P1 S1 C3 40 Ca52 Energy first-order sensitivity S1 1 : -441.09±112.09 MeV Sensitivity % Energy total sensitivity ST 10000 20 Count 5000 0 0 C1 pp S C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11 S1 C3D1 3 P2 c1 c3 c4 cD cE 750 500 250 0 C3 C1 C3 C3 S0 P0 P1 P1 S1 C3 40 Ca54 Energy first-order sensitivity S1 1 : -439.97±109.19 MeV Sensitivity % Energy total sensitivity ST 10000 20 Count 5000 0 0 C1 pp S C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11 S1 C3D1 3 P2 c1 c3 c4 cD cE 750 500 250 0 C3 C1 C3 C3 S0 P0 P1 P1 S1 C3 40 Ca56 Energy first-order sensitivity S1 1 : -433.56±107.65 MeV Sensitivity % Energy total sensitivity ST 10000 20 Count 5000 0 0 C1 pp S C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11 S1 C3D1 3 P2 c1 c3 c4 cD cE 750 500 250 C3 C1 C3 C3 S0 P0 P1 P1 S1 C3 Figure 8.5 Summary of NNLOsat (450) LEC sensitivity as sampled by the RBF-r(50)IEMs cor- responding to each nucleus in neutron-rich Calcium isotopes, based on samples of interactions around the NNLOsat (450) baseline. The left column plots of the sensitivities of each nucleus per LEC. The right column plots the total variance in the energy as measured by the model for that nucleus. The shaded region in the variance plot represents the 1𝜎 and 2𝜎 range. sensitivities in Figure 8.5. 121 40 E(Ca50 Energy- first-order Ca48) sensitivity S1 1 : 35.86±57.90 MeV Sensitivity % Energy total sensitivity ST 10000 20 Count 5000 0 0 C1 pp S C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11 S1 C3D1 3 P2 c1 c3 c4 cD cE 0 200 C3 C1 C3 C3 S0 P0 P1 P1 S1 C3 20 E(Ca52 - Ca50) Energy first-order sensitivity S1 1 : -14.62±11.11 MeV Sensitivity % Energy total sensitivity ST 10000 10 Count 5000 0 0 C1 pp S C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11 S1 C3D1 3 P2 c1 c3 c4 cD cE 50 0 C3 C1 C3 C3 S0 P0 P1 P1 S1 C3 E(Ca54 Energy- first-order Ca52) sensitivity S1 1 : 1.12±14.30 MeV 30 Sensitivity % Energy total sensitivity ST 10000 20 Count 5000 10 0 0 C1 pp S C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11 S1 C3D1 3 P2 c1 c3 c4 cD cE 50 0 50 100 C3 C1 C3 C3 S0 P0 P1 P1 S1 C3 30 E(Ca56 Energy- first-order Ca54) sensitivity S1 1 : 6.41±13.18 MeV Sensitivity % Energy total sensitivity ST 10000 20 Count 5000 10 0 0 C1 pp S C1 0np S C1 n0n S C3 0pp S C3 1np S C3 n1n S C11 S1 C3D1 3 P2 c1 c3 c4 cD cE 50 0 50 C3 C1 C3 C3 S0 P0 P1 P1 S1 C3 Figure 8.6 Summary of two-neutron separation energies (𝑆2𝑛 ) as sampled by the RBF-r(50)IEMs corresponding to each nucleus in neutron-rich Calcium isotopes, based on samples of interactions around the NNLOsat (450) baseline. The left column plots the sensitivities of each nucleus per LEC. The right column plots the total variance in the energy as measured by the model for that nucleus. The shaded region in the variance plot represents the 1𝜎 and 2𝜎 range. 122 8.4 Discussion and Outlook The right panel total variances reveal important information about the fitted RBF-rIEM for each nucleus. First of all, since in every case the LEC variation samples at a fixed width around standard LEC sets (Tables 4.1 and 4.2), we would expect the most probable value, e.g. the peak of the distribution, to reflect the IMSRG(2) outcome for the interactions associated with the standard LEC sets. However, there seems to be systematic error in the expected result from each emulator; in particular, the expected energies are consistently overbound relative the IMSRG(2) result (see Table 8.1). The reason seems to be an overbinding bias in the training sets generated to train the RBF-rIEM for each nucleus. We measured over/underbinding bias by checking the proportion of IMSRG(2) NNLO training points which resulted in an evolved energy outside an arbitrary fixed energy width of ±20% relative to the IMSRG(2) evolved energy using the standard LEC set. Then, an overbound data point is 𝐸 EMU < 𝐸 IMSRG(2) − 20%, whereas an underbound data point 𝐸 EMU > 𝐸 IMSRG(2) + 20%. A “safe” data point is inside the 20% tolerance. We found that the training sets for the calcium isotopes in Figures 8.2 and 8.5 typically presented more overbound data points relative to safe and underbound data points. We did not check the bias for the N3LO training points for 16 O explicitly (see Figure 8.1) although the safe assumption is that this training set presented an underbinding bias, because the mode of the distribution is above the baseline result from Hergert et al. [5], who report a converged energy of 𝐸 = −126.01 MeV. The takeaway message is that the emulator is successfully capturing the sensitivity in the chosen model space, as expected. By exploring the robustness of the emulator in the number of training points, varying observables and training space reduction rank, we have shown that these biases are likely not caused by the emulation machinery, but the crude uniform prior from which we sampled the training points. This motivates better choices in the sampling distribution. In the future, we intend to implement a Bayesian history matching approach akin to Vernon et al. [11], which has been used in recent nuclear physics applications [7]. In this approach, a training set is constructed from so-called non- implausible interactions by validating the corresponding randomly-sampled LEC sets according to 123 (A)Ca 𝐸 IMSRG(2) (MeV) 𝐸 exp (Mev) 48 -403.351178 -416.001197 (19) 50 -411.5000324 -427.508 (2) 52 -419.1716484 -438.32786 (68) 54 -421.911536 -445.37 (5) 56 -423.3308106 -449.85 (2) (A)O ΔNNLOGO (394) 𝐸 IMSRG(2) (MeV) 𝐸 NNDC (Mev) 16 -108.8349651 127.6193152 (32) Table 8.1 IMSRG(2) evolved energies according to the ΔNNLOGO (394)interaction. The experi- mental values from the NNDC are tabulated here for reference, but are not necessarily relevant to the emulator results. some implausibility criteria. Ideally, we use data for selected observables and expert knowledge to define the implausibility measure in appropriate fashion. Naively, we could validate a non-implausible LEC set based on the converged energy of IMSRG(2) evolutions for selected nuclei and families of interaction; however, full IMSRG(2) calculations are expensive. Instead, we might choose to determine non-implausibility by evaluating NN scattering phase-shifts which result from each LEC sample. Once the training data is “clean” in the implausibility sense, we may employ importance re-sampling to dynamically update the LEC sampling distributions to more realistic ranges, based on the results from the GSAs in the previous sections. The goal of combining these methods is to progress toward a quantifiable theoretical error on a prediction provided by the emulator. 124 BIBLIOGRAPHY [1] A. Ekström et al. “Δ isobars and nuclear saturation”. In: Phys. Rev. C 97 (2 Feb. 2018), p. 024332. doi: 10.1103/PhysRevC.97.024332. url: https://link.aps.org/doi/10.1103/ PhysRevC.97.024332. [2] A. Ekström et al. “Accurate nuclear radii and binding energies from a chiral interaction”. In: Phys. Rev. C 91 (5 May 2015), p. 051301. doi: 10.1103/PhysRevC.91.051301. url: https://link.aps.org/doi/10.1103/PhysRevC.91.051301. [3] Andreas Ekström and Gaute Hagen. “Global Sensitivity Analysis of Bulk Properties of an Atomic Nucleus”. In: Physical Review Letters 123.25 (Dec. 2019). doi: 10.1103/ physrevlett.123.252501. url: https://doi.org/10.1103%5C%2Fphysrevlett.123.252501. [4] D. R. Entem and R. Machleidt. “Accurate charge-dependent nucleon-nucleon potential at fourth order of chiral perturbation theory”. In: Physical Review C 68.4 (Oct. 2003). doi: 10.1103/physrevc.68.041001. url: https://doi.org/10.1103%5C%2Fphysrevc.68.041001. [5] H. Hergert et al. “The In-Medium Similarity Renormalization Group: A novel ab initio method for nuclei”. In: Physics Reports 621 (Mar. 2016), pp. 165–222. issn: 0370- 1573. doi: 10.1016/j.physrep.2015.12.007. url: http://dx.doi.org/10.1016/j.physrep.2015. 12.007. [6] Jon Herman and Will Usher. “SALib: An open-source Python library for Sensitivity Analysis”. In: The Journal of Open Source Software 2.9 (Jan. 2017). doi: 10.21105/joss. 00097. url: https://doi.org/10.21105/joss.00097. [7] Baishan Hu et al. “Ab initio predictions link the neutron skin of 208Pb to nuclear forces”. In: Nature Physics 18.10 (Aug. 2022), pp. 1196–1200. doi: 10.1038/s41567-022-01715-8. url: https://doi.org/10.1038%5C%2Fs41567-022-01715-8. [8] Takuya Iwanaga, William Usher, and Jonathan Herman. “Toward SALib 2.0: Advancing the accessibility and interpretability of global sensitivity analyses”. In: Socio-Environmental Systems Modelling 4 (May 2022), p. 18155. doi: 10.18174/sesmo.18155. url: https: //sesmo.org/article/view/18155. [9] Andrea Saltelli et al. “Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index”. In: Computer Physics Communications 181.2 (2010), pp. 259–270. issn: 0010-4655. doi: https://doi.org/10.1016/j.cpc.2009.09.018. url: https://www.sciencedirect.com/science/article/pii/S0010465509003087. [10] I.M Sobol´. “Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates”. In: Mathematics and Computers in Simulation 55.1 (2001). The Second IMACS Seminar on Monte Carlo Methods, pp. 271–280. issn: 0378-4754. doi: https: 125 //doi.org/10.1016/S0378-4754(00)00270-6. url: https://www.sciencedirect.com/science/ article/pii/S0378475400002706. [11] Ian Vernon et al. “Bayesian Emulation and History Matching of JUNE”. in: medRxiv (2022). doi: 10.1101/2022.02.21.22271249. eprint: https://www.medrxiv.org/content/early/2022/ 02/22/2022.02.21.22271249.full.pdf. url: https://www.medrxiv.org/content/early/2022/ 02/22/2022.02.21.22271249. 126 CHAPTER 9 CONCLUSIONS AND OUTLOOK In this work, we have made significant progress toward the theoretical and computational im- provement of the IMSRG, as outlined in Chapter 1. Here, we summarize the main scientific contributions presented in this work and outline the future directions that naturally arise from the studies conducted therein. Chapter 5 presented a tensor network architecture scheme for factorizing data structures relevant the to numerical solution of the IMSRG flow. We showed that the tensor-train decomposition can be used to factorize the occupation factors, which are constant in the flow, thereby reducing the memory cost (e.g. 𝑁 2 to 𝑁 in some cases) of storing on disk. For large single particle bases, this can result in several orders of magnitude of memory savings. We recast the IMSRG flow equations as tensor contractions to be performed by a tensor library, which when combined with factorization of the occupation factors has favorable scaling over naive loop implementations in a toy model. Success here motivates the continued exploration of reduced representations of the IMSRG itself. In the future, we wish to leverage the performance improvements provided by these techniques in the IMSRG production codes. Tensor library integration will also guarantee that the IMSRG codes take advantage of improvements developed by the maintainers of the library, so that the IMSRG calculation is always as efficient as possible. Improvements to the speed and efficiency of the IMSRG production codes would carry into all other advances presented in this work, since they rely on full-scale IMSRG calculations, e.g., for the construction of training data. Chapter 6 presented our main methodological improvement to the IMSRG, called the reference state ensemble. We formulated the reference state ensemble as a mixture of independent one-body density matrices for selected Slater determinant configurations, which are weighted according to tunable probabilities. We showed that we can reduce error due to truncation of induced terms by informing the standard IMSRG operator basis with excited state information using the ensemble— the evidence being how the reference ensemble stabilizes the spectrum of an evolving Hamiltonian, 127 where a corresponding single reference IMSRG flow would diverge. Ultimately, we use the reference ensemble to attempt to maintain the unitarity of the IMSRG flow as measured by the invariance of the Hamiltonian’s spectrum. Future work in this area will involve a new formulation of the optimization problem using a proxy observable (i.e. not the spectrum) whose evaluation does not require costly exact diagonalizations of the flowing Hamiltonian. Following this breakthrough, we may extend the reference ensemble to realistic nuclear calculations, as well as integrate it into our emulator techniques presented in Chapter 7. Chapter 7 presented our main computational improvement to the IMSRG, which is a method for emulating the evolution with the help of the so-called Dynamic Mode Decomposition (DMD). We showed that the DMD expansion, informed by direct measurements of the IMSRG-evolved Hamil- tonian, may be interpreted as approximating the unitary transformation generated by integration of the IMSRG flow equations. Using this emulator technique, we presented a method for emulating the parametric IMSRG based on an “interpolation engine” powered by DMD. We improved the interpolation method further by effectively performing a Principal Component Analysis (PCA) to reduce the feature space of the IMSRG training points—and we showed, critically, how to imple- ment the PCA calculation via a streaming Singular Value Decomposition in order to mitigate the memory constraints that come with larger realistic IMSRG model spaces. The IE generalized well enough to the full parameter space to perform sensitivity analyses on the parameters of the chiral interactions, the LECs, which we presented in Chapter 8. We showed that a statistically significant number of calculations for GSA could be performed in just a few minutes using our parametric emulator, whereas the same number of full IMSRG(2) calculations would take thousands of years. Future work will attempt to apply the DMD method to the Magnus formulation of the IMSRG, which has the critical advantage that the unitary transformation driving the IMSRG flow is mea- sured directly, instead of implicitly through measurements of the flowing Hamiltonian (as is the case in the standard IMSRG). Emulating the Magnus operator has the advantage that we may easily emulate the flow of observables besides the Hamiltonian, without ever taking direct measurements of their evolution. Additionally, we wish to improve the parametric emulator by informing the 128 baseline, uniform prior with nonimplausibility criteria, e.g., through history-matching, in order to curtail the training set biases and optimize the useful physical information present in the training space. In this way, we may draw better, physically-motivated conclusions from the sensitivity analyses. Improvements in all these areas will take us closer to proper theoretical error bars on the IMSRG results, which will allow us to use the method to address scientific challenges in nuclear physics, nuclear astrophysics, and the exploration of the fundamental symmetries of the Standard Model of Particle Physics. 129 APPENDIX A TENSOR TRAIN DECOMPOSITION OF THE OCCUPATION FACTORS A.1 Tensor Train Decompositions for Occupation Factors in the IMSRG(2)   h i 1  𝐴(𝑎, 𝑏) = 𝑛𝑎 − 𝑛 𝑏 = 𝑛𝑎 1     (A.1) −𝑛   𝑏     h i 1  𝐵(𝑎, 𝑏) = 1 − 𝑛𝑎 − 𝑛 𝑏 = 1 − 𝑛𝑎 1     (A.2) −𝑛   𝑏   𝐶 (1) (𝑎, 𝑏, 𝑐) = 𝑛𝑎 𝑛 𝑏 𝑛¯ 𝑐 + 𝑛¯ 𝑎 𝑛¯ 𝑏 𝑛𝑐 = 𝑛𝑎 𝑛 𝑏 + 𝑛𝑐 − 𝑛 𝑏 𝑛𝑐 − 𝑛𝑎 𝑛𝑐    𝑛 𝑏  0 0 0   1     h i 0 1 0 0  𝑛𝑐  (A.3) = 𝑛𝑎  1 1 𝑛𝑎    0  0 −𝑛 𝑏 0  𝑛𝑐        0 0 0 −1 𝑛𝑐     𝐶 (2) (𝑎, 𝑏, 𝑐) = 𝑛𝑎 𝑛¯ 𝑏 𝑛¯ 𝑐 + 𝑛¯ 𝑎 𝑛 𝑏 𝑛𝑐 = 𝑛 𝑏 𝑛𝑐 + 𝑛𝑎 − 𝑛𝑎 𝑛𝑐 − 𝑛𝑎 𝑛 𝑏    𝑛 𝑐  0 0 0   1     h i 0 1 0 0  𝑛𝑎  (A.4) = 𝑛𝑏  1 1 𝑛𝑏    0  0 −𝑛𝑐 0  𝑛𝑎        0 0 0 −1 𝑛𝑎         h i 𝑛 𝑏  h i 1 − 𝑛 𝑑  𝐷 (𝑎, 𝑏, 𝑐, 𝑑) = 𝑛𝑎 𝑛 𝑏 𝑛¯ 𝑐 𝑛¯ 𝑑 = 𝑛𝑎     0   1 − 𝑛𝑐 0   (A.5) 0  0          130 𝐷 (2) (𝑎, 𝑏, 𝑐, 𝑑) = 𝑛¯ 𝑎 𝑛¯ 𝑏 𝑛𝑐 𝑛 𝑑 − 𝑛𝑎 𝑛 𝑏 𝑛¯ 𝑐 𝑛¯ 𝑑 =   1 0 0 0 0 0     0 1 0 0 0 0     i 0 0 −𝑛 𝑏  h 0 0 0 = 1 −𝑛𝑎 1 −𝑛𝑎 𝑛𝑎 𝑛𝑎  ×  0 0 0 𝑛𝑏 0 0     0 0 0 0 𝑛𝑏 0       0 0 0 0 0 𝑛𝑏       𝑛 𝑐 0 0 0 0 0 𝑛 𝑑      0 𝑛 0 0 0 0 𝑛 𝑑   𝑐        0 0 𝑛𝑐 0 0 0 𝑛 𝑑        0 0 0 1 0 0  1      0 0 0 0 𝑛𝑐 0  1         0 0 0 0 0 1  𝑛 𝑑     (A.6)       h i 𝑛 𝑏  h i 1 − 𝑛 𝑑  h i 1 − 𝑛 𝑓  𝐸 (𝑎, 𝑏, 𝑐, 𝑑, 𝑒, 𝑓 ) = 𝑛𝑎 𝑛 𝑏 𝑛𝑐 𝑛¯ 𝑑 𝑛¯ 𝑒 𝑛¯ 𝑓 = 𝑛𝑎       0   𝑛𝑐 0   1 − 𝑛𝑒 0   (A.7) 0  0   0              131 APPENDIX B MULTI-REFERENCE IMSRG(2) FLOW EQUATIONS B.1 Multi-Reference IMSRG(2) Flow Equations The multi-reference IMSRG(2) flow equations include irreducible density matrices (IDMs), introduced in Chapter 2.6, which encode correlations in the reference state that couple to the correlations described by the IMSRG flow. Here, we record the MR-IMSRG(2) flow equations as they are introduced by Hergert [1]: 𝑑𝐸 ∑︁ 1 ∑︁ = (𝑛𝑎 − 𝑛 𝑏 )𝜂𝑎𝑏 𝑓𝑏𝑎 + (𝜂𝑎𝑏𝑐𝑑 Γ𝑐𝑑𝑎𝑏 − Γ𝑎𝑏𝑐𝑑 𝜂𝑐𝑑𝑎𝑏 ) 𝑛𝑎 𝑛 𝑏 𝑛𝑐 𝑛 𝑑 𝑑𝑠 𝑎𝑏 4 𝑎𝑏𝑐𝑑   (B.1) 1 ∑︁ 𝑑 1 ∑︁ + Γ𝑎𝑏𝑐𝑑 𝜆 𝑎𝑏𝑐𝑑 + (𝜂𝑎𝑏𝑐𝑑 Γ𝑘𝑙𝑎𝑚 − Γ𝑎𝑏𝑐𝑑 𝜂 𝑘𝑙𝑎𝑚 ) 𝜆 𝑏𝑘𝑙𝑐𝑑𝑚 , 4 𝑎𝑏𝑐𝑑 𝑑𝑠 4 𝑎𝑏𝑐𝑑𝑘𝑙𝑚 𝑑𝑓𝑖 𝑗 ∑︁ ∑︁  = (𝜂𝑖𝑎 𝑓𝑎 𝑗 − 𝑓𝑖𝑎 𝜂𝑎 𝑗 ) + 𝜂𝑎𝑏 Γ𝑏𝑖𝑎 𝑗 − 𝑓𝑎𝑏 𝜂 𝑏𝑖𝑎 𝑗 (𝑛𝑎 − 𝑛 𝑏 ) 𝑑𝑠 𝑎 𝑎𝑏 1 ∑︁  + 𝜂𝑖𝑎𝑏𝑐 Γ𝑏𝑐 𝑗 𝑎 − Γ𝑖𝑎𝑏𝑐 𝜂 𝑏𝑐 𝑗 𝑎 (𝑛𝑎 𝑛 𝑏 𝑛𝑐 + 𝑛𝑎 𝑛 𝑏 𝑛𝑐 ) 2 𝑎𝑏𝑐 (B.2) 1 ∑︁ ∑︁ + (𝜂𝑖𝑎𝑏𝑐 Γ𝑑𝑒 𝑗 𝑎 − Γ𝑖𝑎𝑏𝑐 𝜂 𝑑𝑒 𝑗 𝑎 )𝜆 𝑑𝑒𝑏𝑐 + (𝜂𝑖𝑎𝑏𝑐 Γ𝑏𝑒 𝑗 𝑑 − Γ𝑖𝑎𝑏𝑐 𝜂 𝑏𝑒 𝑗 𝑑 )𝜆 𝑎𝑒𝑐𝑑 4 𝑎𝑏𝑐𝑑𝑒 𝑎𝑏𝑐𝑑𝑒 1 ∑︁ 1 ∑︁ − (𝜂𝑖𝑎 𝑗 𝑏 Γ𝑐𝑑𝑎𝑒 − Γ𝑖𝑎 𝑗 𝑏 𝜂𝑐𝑑𝑎𝑒 )𝜆 𝑐𝑑𝑏𝑒 + (𝜂𝑖𝑎 𝑗 𝑏 Γ𝑏𝑐𝑑𝑒 − Γ𝑖𝑎 𝑗 𝑏 𝜂 𝑏𝑐𝑑𝑒 )𝜆 𝑎𝑐𝑑𝑒 , 2 𝑎𝑏𝑐𝑑𝑒 2 𝑎𝑏𝑐𝑑𝑒 𝑑Γ𝑖 𝑗 𝑘𝑙 ∑︁ = (𝜂𝑖𝑎 Γ𝑎 𝑗 𝑘𝑙 + 𝜂 𝑗 𝑎 Γ𝑖𝑎𝑘𝑙 − 𝜂𝑎𝑘 Γ𝑖 𝑗 𝑎𝑙 − 𝜂𝑎𝑙 Γ𝑖 𝑗 𝑘𝑎 − 𝑓𝑖𝑎 𝜂𝑎 𝑗 𝑘𝑙 − 𝑓 𝑗 𝑎 𝜂𝑖𝑎𝑘𝑙 + 𝑓𝑎𝑘 𝜂𝑖 𝑗 𝑎𝑙 + 𝑓𝑎𝑙 𝜂𝑖 𝑗 𝑘𝑎 ) 𝑑𝑠 𝑎 1 ∑︁ + (𝜂𝑖 𝑗 𝑎𝑏 Γ𝑎𝑏𝑘𝑙 − Γ𝑖 𝑗 𝑎𝑏 𝜂𝑎𝑏𝑘𝑙 )(1 − 𝑛𝑎 − 𝑛 𝑏 ) 2 𝑎𝑏 ∑︁   + (𝑛𝑎 − 𝑛 𝑏 ) 𝜂𝑖𝑎𝑘 𝑏 Γ 𝑗 𝑏𝑙𝑎 − Γ𝑖𝑎𝑘 𝑏 𝜂 𝑗 𝑏𝑙𝑎 − 𝜂 𝑗 𝑎𝑘 𝑏 Γ𝑖𝑏𝑙𝑎 − Γ 𝑗 𝑎𝑘 𝑏 𝜂𝑖𝑏𝑙𝑎 . 𝑎𝑏 (B.3) Note that these equations reduce to the single reference flow equations, in Equation 3.18, when we set the IDMs equal to zero. In contrast to the regular IMSRG(2), the MR-IMSRG(2) naively 132 scales as 𝑂 (𝑁 7 ) when the term proportional to the three-body IDM is included in the zero-body flow equation B.1. For common types of reference states, we can often exploit the specific structure of the IDMs to reduce the scaling of the associated terms in the flow equations, and thereby recover the same 𝑂 (𝑁 6 ) computational cost as in the conventional IMSRG(2). Evolving the MR-IMSRG(2) flow equations usually requires one to two orders of magnitude more integration steps [1]. 133 BIBLIOGRAPHY [1] H Hergert. “In-medium similarity renormalization group for closed and open-shell nu- clei”. In: Physica Scripta 92.2 (Dec. 2016), p. 023002. issn: 1402-4896. doi: 10.1088/1402-4896/92/2/023002. url: http://dx.doi.org/10.1088/1402-4896/92/2/023002. 134 APPENDIX C CONTRIBUTED CODE PROJECTS In this appendix, we record the code projects that were created, developed, or improved upon as contribution to this thesis. These code projects were also used in the studies presented within the main body of this work. C.1 IMSRG Codes The codes in Table C.1 were developed to implement the IMSRG(2) flow of schematic Hamil- tonians like the P3H model. One project was written in Python, the other was written in C++. Although originally developed for the P3H model, both codes have been designed with modularity in mind; in principle, one can apply either code to a different model Hamiltonian or Observable with the addition of a new Hamiltonian or Observable class that implements the required features. In chapter 6, this functionality was used in the Python code to evolve and evaluate the total spin-squared operator 𝑆ˆ2 alongside the P3H Hamiltonian. The Python project Tensor-factor IMSRG (TFIMSRG) outsources tensor contraction at the backend to the library TensorNetwork [2]. The C++ project Tensorized C++ IMSRG (TCIMSRG) employs the Tensor Algebra Compiler (TACO) [1] as a backend for tensor contraction computation. Both codes implement the tensor-train decomposition of occupation factors described in Chapter 5.3 and presented explicitly in Appendix A. The summations for the IMSRG flow equations have been converted to tensor contractions and computed with these libraries. Future work may include the extension of TCIMSRG with a new, distributed version of TACO called Distributed Tensor Algebra Compiler (DISTAL) [3], as well as the deployment of the tensor network framework in a new IMSRG production code. The P3H model IMSRG codes in Table C.1 are available at https://github.com/davis9ja/tfimsrg and https://github.com/davis9ja/tcimsrg. The IMSRG production code, which implements the IMSRG for realistic nuclear Hamiltonians, is available on request to the author, H. Hergert. 135 Project Language Dependencies Description Name tfimsrg Python PyCI, Implements the IMSRG(2) flow for the TensorNetwork1 P3H model. Object oriented approach. Machinery for IMSRG(3) implemented, but not ready for experiment. tcimsrg C++ dmcpp, TACO2 Implements the IMSRG(2) flow for the P3H model. Object oriented approach. IMSRG C MEQ, ME3BX, HFB, Feature-complete production code project GSL3, SUNDIALS4 which was updated with emulation func- tionality as part of this dissertation. Table C.1 Table of IMSRG codes developed for this dissertation. The listed dependencies are non-standard libraries or code projects which must be installed and are listed in Tables C.2 and C.3. C.2 Full Configuration Interaction Code for Schematic Hamiltonians This code implements the Full Configuration Interaction approach, i.e., exact diagonalization, for schematic Hamiltonians, as described in Chapter 2.10. The Hamiltonian is fully constructed and diagonalized in the many-body basis. This is often possible for a schematic model like the P3H model. For realistic interactions, the model space dimensions required for the FCI approach quickly become prohibitive, even if the sparsity of the Hamiltonian is exploited and computations are limited to a few low-lying states, e.g., through the use of Lanczos or Arnoldi methods. In the present work, we used this code for diagnostic purposes in IMSRG calculations with reference ensembles, as discussed in Chapter 2.10. The code was also used to study the structure of exact eigenstate of the P3H Hamiltonian as a function of its parameters (e.g., for the determination of phase diagrams of this model). We developed two FCI implementations for this work: one in Python and one in C++. The Python code, located at https://gitlab.msu.edu/imsrg/pyci/-/tree/jacob-dev/imsrg_ci is forked from 1 Tensor contraction library that wraps TensorFlow, among other backends: https://github.com/google/ TensorNetwork/ [2] 2 Tensor Algebra Compiler (TACO) which performs efficient tensor contraction in C++: https://github.com/ tensor-compiler/taco [1] 3 GNU Scientific Library, https://www.gnu.org/software/gsl/ 4 SUite of Nonlinear and DIfferential/ALgebraic equation Solvers, https://github.com/LLNL/sundials 136 a rudimentary version originally developed by A. Zitzelberger, another student in the working group. The C++ code is located at https://github.com/davis9ja/cci. Project Language Dependencies Description Name PyCI Python None Constructs the P3H Hamiltonian in a Slater determinant basis via FCI. CCI C++ DMCPP Constructs the P3H Hamiltonian in a Slater determinant basis via FCI. Table C.2 Table of FCI codes developed for this dissertation. Refer to Chapter 2.10 for information on FCI. C.3 Density Matrix Codes Both P3H IMSRG codes in Table C.1 implement reference state ensemble functionality, as well as the normal ordering with respect to exact FCI eigenstates. This requires the computation of density matrices and their use in normal ordering schemes. The following codes may be used to generate one-, two-, and three-body density matrices, according to a user-specified Slater determinant basis, for use with reference state ensemble studies. The DMCPP code is located at https://github.com/davis9ja/dmcpp. Project Language Dependencies Description Name PyCI Python None The density_matrix module may be ex- ported to compute density matrices inde- pendent of FCI use. DMCPP C++ None Compute density matrices for specified SD basis. Compiled as a static library and intended for import to IMSRG C++ code. Table C.3 Table of density matrix codes developed for this dissertation. These codes are necessary for generating the density matrices required for generalized normal-ordering. C.4 IMSRG Emulation using Dynamic Mode Decomposition The functionality for emulating an IMSRG flow using Dynamic Mode Decomposition (DMD) is wrapped into the IMSRG production code as well as TCIMSRG. For independent implementations 137 and parametric emulation, the package IMSRG-EMU was developed as part of this dissertation. The package is written in Python and provides a simple command-line interface for running emulations on archived data; however, the classes which perform the emulations can be exported to user Python codes for additional functionality and independent experimentation. The IMSRG-EMU code is located at https://github.com/davis9ja/imsrg_emu. Project Language Dependencies Description Name IMSRG-EMU Python None Provides interface for computing a stan- dard, projected DMD emulation on an in- put dataset (Chapter 7.2) as well as com- puting a parametric emulation with train- ing on an appropriate training set (Chapter 7.3). Table C.4 Table of IMSRG emulation codes, using DMD, developed as contribution to this dissertation work. C.5 Full Table of Contributed Codes For reference, we have compiled all contributed codes outlined in Appendix sections C.1 through C.4 into a single table. Refer to each code’s user documentation for a complete description of their syntax and usage. Project Language Dependencies Location Name TFIMSRG Python PyCI, https://github.com/davis9ja/tfimsrg TensorNetwork TCIMSRG C++ dmcpp, TACO https://github.com/davis9ja/tcimsrg IMSRG C MEQ, ME3BX, HFB, Available upon request to the author, H. GSL, SUNDIALS Hergert PyCI Python None https://gitlab.msu.edu/imsrg/pyci/-/tree/ jacob-dev/imsrg_ci CCI C++ DMCPP https://github.com/davis9ja/cci DMCPP C++ None https://github.com/davis9ja/dmcpp IMSRG-EMU Python None https://github.com/davis9ja/imsrg_emu Table C.5 Complete list of code projects and contributed extensions developed as part of this dissertation. 138 BIBLIOGRAPHY [1] Fredrik Kjolstad et al. “The Tensor Algebra Compiler”. In: Proc. ACM Program. Lang. 1.OOPSLA (Oct. 2017), 77:1–77:29. issn: 2475-1421. doi: 10.1145/3133901. url: http://doi.acm.org/10.1145/3133901. [2] Chase Roberts et al. TensorNetwork: A Library for Physics and Machine Learning. 2019. arXiv: 1905.01330 [physics.comp-ph]. [3] Rohan Yadav, Alex Aiken, and Fredrik Kjolstad. “DISTAL: The Distributed Tensor Algebra Compiler”. In: Proceedings of the 43rd ACM SIGPLAN International Conference on Programming Language Design and Implementation. PLDI 2022. San Diego, CA, USA: Association for Computing Machinery, 2022, pp. 286–300. isbn: 9781450392655. doi: 10.1145/3519939.3523437. url: https://doi.org/10.1145/3519939.3523437. 139