Nuclear matter from chiral effective field theory Kernmaterie basierend auf chiraler effektiver Feldtheorie Zur Erlangung des Grades eines Doktors der Naturwissenschaften (Dr. rer. nat.) genehmigte Dissertation von Christian Drischler aus Rüsselsheim Tag der Einreichung: 17.10.2017, Tag der Prüfung: 15.11.2017 Darmstadt 2017 — D 17 1. Gutachten: Prof. Ph.D. Achim Schwenk 2. Gutachten: Prof. Dr. Jens Braun Fachbereich Physik Institut für Kernphysik Theoriezentrum Nuclear matter from chiral effective field theory Kernmaterie basierend auf chiraler effektiver Feldtheorie Genehmigte Dissertation von Christian Drischler aus Rüsselsheim 1. Gutachten: Prof. Ph.D. Achim Schwenk 2. Gutachten: Prof. Dr. Jens Braun Tag der Einreichung: 17.10.2017 Tag der Prüfung: 15.11.2017 Darmstadt 2017 — D 17 Bitte zitieren Sie dieses Dokument als: URN: urn:nbn:de:tuda-tuprints-69845 URL: http://tuprints.ulb.tu-darmstadt.de/id/eprint/6984 Dieses Dokument wird bereitgestellt von tuprints, E-Publishing-Service der TU Darmstadt http://tuprints.ulb.tu-darmstadt.de tuprints@ulb.tu-darmstadt.de Die Veröffentlichung steht unter folgender Creative Commons Lizenz: Namensnennung – Keine kommerzielle Nutzung – Keine Bearbeitung 4.0 international https://creativecommons.org/licenses/by-nc-nd/4.0/ Abstract Nuclear matter is an ideal theoretical system that provides key insights into the physics of different length scales. While recent ab initio calculations of medium-mass to heavy nuclei have demonstrated that realistic saturation properties in infinite matter are crucial for reproducing experimental binding energies and charge radii, the nuclear-matter equation of state allows tight constraints on key quantities of neutron stars. In the present thesis we take advantage of both aspects. Chiral effective field theory (EFT) with pion and nucleon degrees of freedom has become the modern low-energy approach to nuclear forces based on the symmetries of quantum chromodynamics, the fun- damental theory of strong interactions. The systematic chiral expansion enables improvable calculations associated with theoretical uncertainty estimates. In recent years, chiral many-body forces were derived up to high orders, allowing consistent calculations including all many-body contributions at next-to-next- to-next-to-leading order (N3LO). Many further advances have driven the construction of novel chiral potentials with different regularization schemes. Here, we develop advanced methods for microscopic calculations of the equation of state of homogeneous nuclear matter with arbitrary proton-to-neutron ratio at zero temperature. Specifically, we push the limits of many-body perturbation theory (MBPT) considerations to high orders in the chiral and in the many-body expansion. To address the challenging inclusion of three-body forces, we introduce a new partial-wave method for normal ordering that generalizes the treatment of these contributions. We show improved predictions for the neutron-matter equation of state with consistent N3LO nucleon-nucleon (NN) plus three-nucleon (3N) potentials using MBPT up to third order and self-consistent Green’s function theory. The latter also provides nonperturbative benchmarks for the many-body convergence. In addition, we extend the normal-ordering method to finite temperatures. Calculations of asymmetric matter require in addition reliable fit values for the low-energy couplings that contribute to the 3N forces. This was not the case for N3LO calculations. We present a novel Monte-Carlo framework for perturbative calculations with two-, three-, and four-nucleon interactions, which, including automatic code generation, allows to compute successive orders in MBPT as well as chiral EFT in an efficient way. The performance is such that it can be used for optimizing next-generation chiral potentials with respect to saturation properties. As a first step in this direction, we study nuclear matter based on chiral low-momentum interactions, exhibiting a very good many-body convergence up to fourth order. We then explore new chiral interactions up to N3LO, where simultaneous fits to the triton and to saturation properties can be achieved with natural 3N low-energy couplings. We perform a comprehensive Weinberg eigenvalue analysis of a representative set of modern local, semilo- cal, and nonlocal chiral NN potentials. Our detailed comparison of Weinberg eigenvalues provides various insights into idiosyncrasies of chiral potentials for different orders and partial waves. We demonstrate that a direct comparison of numerical cutoff values of different interactions is in general misleading due to the different analytic form of regulators. This shows that Weinberg eigenvalues also can be used as a helpful monitoring scheme when constructing new interactions. Furthermore, we present solutions of the BCS gap equation in the channels 1S0 and 3P2−3F2 in neutron matter. Our studies are based on nonlocal NN plus 3N interactions up to N3LO as well as the aforemen- tioned local and semilocal chiral NN interactions up to N2LO and N4LO, respectively. In particular, we investigate the impact of N3LO 3N forces on pairing gaps and also derive uncertainty estimates by taking into account results at different orders in the chiral expansion. In addition, different methods for obtaining self-consistent solutions of the gap equation are discussed. Besides the widely-used quasilinear method we demonstrate that the modified Broyden method is well applicable and exhibits a robust convergence behavior. Zusammenfassung Kernmaterie ist ein ideales, theoretisches System, das zentrale Einblicke in die Physik verschiedener Län- genskalen gewährt. Während ab initio Rechnungen von mittelschweren bis schweren Kerne vor kurzem gezeigt haben, dass realistische Saturierungseigenschaften in unendlicher Materie für die experimentalen Bindungsenergien und Kernradien wichtig sind, erlaubt die Zustandsgleichung von Kernmaterie starke Einschränkungen von Schlüsselgrößen für Neutronensternen. In der vorliegenden Dissertation machen wir uns beide Aspekte zu nutze. Die chirale effektive Feldtheorie (EFT), mit Pionen und Nukleonen als Freiheitsgraden, wurde zum moder- nen, niederenergetischen Zugang zu Kernkräften, basierend auf den Symmetrien der Quantenchromody- namik, der fundamentalen Theorie der Starken Wechselwirkung. Die systematische, chirale Entwicklung ermöglicht verbesserbare Rechnungen mit theoretischen Unsicherheitsabschätzungen. In den letzten Jah- ren wurden chirale Vielteilchenkräfte bis in hohe Ordnungen ausgearbeitet. Konsistente Rechnungen mit allen Vielteilchenbeiträgen auf der sogenannten „next-to-next-to-next-to-leading order“ (N3LO) sind daher möglich. Viele weitere Fortschritte führten zu der Entwicklung von neuen chiralen Potentialen innerhalb verschiedener Regularisierungschemen. Wir entwickeln fortgeschrittene Methoden für mikroskopische Berechnungen der Zustandsgleichung von homogener Kernmaterie mit beliebigem Proton-zu-Neutron-Verhältnis am Temperatur-Nullpunkt. Im Spe- ziellen setzen wir die Maßstäbe für Vielteilchen-störungstheoretische Ansätze zu hohen Ordnungen in der chiralen und der Vielteilchen-Entwicklung. Um die schwierige Einbindung von Dreinukleon-Kräften zu ermöglichen, führen wir eine neue Partialwellen-Methode für die Normalordnung ein, welche die Behandlung dieser Beiträge verallgemeinert. Wir treffen verbesserte Vorhersagen für die Zustandsglei- chung von Neutronenmaterie mit konsistenten N3LO Nukleon-Nukleon- (NN) plus Dreinukleon- (3N) Potentialen bis zur dritten Ordnung in der Vielteilchenstörungstheorie und in der Methode der selbstkon- sistenten Greens-Funktionen. Letztere erlaubt nichtperturbative Vergleichsrechnungen zum Bestimmen der Vielteilchenkonvergenz. Wir erweitern außerdem die Methode zur Normalordnung zu endlichen Temperaturen. Berechnungen von asymmetrischer Materie setzen zusätzlich zuverlässige Fitwerte für die Niederenergie- kopplungen, die zu den 3N Kräften beitragen, voraus. Das war bisher nicht der Fall für N3LO Rechnun- gen. Wir zeigen ein neues Monte-Carlo Framework für perturbative Rechnungen mit Zwei-, Drei- und Viernukleonen-Wechselwirkungen, das es auf effiziente Weise erlaubt (zusammen mit dem automatischen Generieren von Quellcode), verschiedene Ordnungen in der Störungstheorie als auch in der chiralen EFT zu berechnen. Die Leistungsfähigkeit der Methode ist so angelegt, dass sie zur Optimierung von neuen chiralen Potentialen im Sinne von Saturierungseigenschaften benutzt werden kann. Als einen ersten Schritt in diese Richtung, studieren wir Kernmaterie basierend auf evolvierten chiralen Wech- selwirkungen, welche eine sehr gute Vielteilchenkonvergenz bis zur vierten Ordnung aufweisen. Wir untersuchen dann neue chirale Wechselwirkungen bis zur N3LO, wobei simultane Fits an das Triton und an Saturierungseigenschaften mit natürlichen 3N-Niederenergiekopplungen erzielt werden können. Wir führen eine umfassende Weinberg-Eigenwert-Analyse mit einer repräsentativen Menge von moder- nen lokalen, semilokalen und nichtlokalen chiralen NN-Wechselwirkungen durch. Unsere detaillierten Vergleiche anhand von Weinberg-Eigenwerten gewähren verschiedenste Einblicke in die Eigenheiten von chiralen Potentialen für verschiedene Ordnungen und Partialwellen. Wie zeigen, dass ein direkter Vergleich von numerischen Cutoff-Werten für verschiedene Wechselwirkungen, aufgrund der unterschied- lichen analytischen Form von Regulatoren, täuschen kann. Das zeigt, Weinberg-Eigenwerte können auch als nützliche Hilfsmittel beim Erstellen von neuen Wechselwirkungen dienen. Wir präsentieren weiterhin Lösungen der BCS-Gleichung der Energielücke in den Kanälen 1S0 und 3P2−3F2 in Neutronenmaterie. Unsere Studien basieren sowohl auf nichtlokalen NN- plus 3N-Wechselwirkungen bis zur N3LO, als auch auf den zuvor genannten lokalen und semilokalen chiralen NN-Wechselwirkungen bis zur N2LO beziehungsweise bis zur N4LO. Im Speziellen untersuchen wir den Einfluss von N3LO 3N-Kräften auf die Energielücke und leiten Unsicherheitsabschätzungen her, in dem wir Ergebnisse auf verschiedenen Ordnungen der chiralen Entwicklung berücksichtigen. Zusätzlich werden verschiedene Methoden zum Lösen der BCS-Gleichung der Energielücke diskutiert. Neben einer oft benutzten quasili- nearen Methode, zeigen wir, dass die modifizierte Broyden-Methode gut anwendbar ist und ein robustes Konvergenzverhalten aufweist. 4 Contents 1 Introduction 7 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2 Chiral effective field theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2.1 Quantum chromodynamics and chiral symmetry . . . . . . . . . . . . . . . . . . . . . 13 1.2.2 Hierarchy of chiral many-body forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.2.3 Theoretical uncertainties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.3 Infinite nuclear matter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.3.1 Many-body perturbation theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.3.2 Normal ordering with 3N forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 1.3.3 Other approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 1.3.4 Similarity renormalization group . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2 Nuclear forces from chiral effective field theory 35 2.1 Nucleon-nucleon interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.1.1 Regularization schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.1.2 Partial-wave decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.2 Three-nucleon interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 2.3 Four-nucleon interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3 Weinberg eigenvalue analysis of chiral NN interactions 53 3.1 Weinberg eigenvalue analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.1.1 Perturbative convergence and eigenvalue equation . . . . . . . . . . . . . . . . . . . . 53 3.1.2 Solving the eigenvalue equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.1.3 Connection to phase shifts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.2 Comprehensive study of local, semilocal, and nonlocal interactions . . . . . . . . . . . . . . 58 3.2.1 Observations from repulsive and attractive eigenvalues . . . . . . . . . . . . . . . . . 58 3.2.2 Interpretation and general conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.3 Evolved potentials and universality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4 Neutron-matter equation of state based on chiral NN and 3N interactions up to N3LO 67 4.1 Improved normal-ordering method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.2 Decomposed energy relations up to third order in MBPT . . . . . . . . . . . . . . . . . . . . . 71 4.2.1 Hartree-Fock level . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.2.2 Second-order contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.2.3 Third-order hole-hole and particle-particle contributions . . . . . . . . . . . . . . . . 75 4.3 Many-body convergence: MBPT vs. SCGF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.4 Full N3LO calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.5 Perspectives for finite-temperature calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5 An efficient Monte-Carlo framework for MBPT calculations 87 5.1 Chiral forces in a single-particle basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.2 Automated computation of energy contributions for infinite-matter calculations . . . . . . . 94 5.3 Multidimensional integration using Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.4 Saturation properties with NN and 3N forces at N3LO . . . . . . . . . . . . . . . . . . . . . . . 100 Contents 5 6 BCS pairing gaps in neutron matter: uncertainties and 3N forces 107 6.1 BCS pairing and the energy gap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.2 Solving the BCS gap equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.2.1 Numerical challenges and direct-iteration method . . . . . . . . . . . . . . . . . . . . 110 6.2.2 Khodel’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6.2.3 Modified Broyden’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 6.2.4 Benchmarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 6.3 Results with local and semilocal NN potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.4 Results with 3N forces up to N3LO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 7 Summary and outlook 123 A Normal-ordering symmetry factors 127 B Partial-wave decomposition of the BCS gap equation 129 C List of abbreviations 131 References 133 Acknowledgments 145 Curriculum Vitae 147 6 Contents 1 Introduction 1.1 Motivation Nuclear matter is an ideal environment for testing nuclear forces with important consequences for finite nuclei as well as nuclear astrophysics in general [1]. The equation of state allows the prediction of key quantities relevant for neutron stars and their formation in core-collapse supernovae [2, 3]. Being able to interconnect the microscopic length scales of nuclei with the macroscopic scales of neutron stars [4] makes studies of the equation of state of nuclear matter very exciting. Pioneered by Steven Weinberg and others in the early 1990’s chiral effective field theory (EFT) nowa- days has become the standard approach to nuclear forces at low-energy scales [5, 6], where momenta are of the order of the pion mass Q ∼ mπ, i.e., lower than the breakdown scale Λb ∼ 500 MeV of the theory. Chiral EFT provides a systematic and improvable expansion in powers of Q/Λb ∼ 1/3 based on the symmetries of the fundamental theory of strong interactions, quantum chromodynamics (QCD). The dominant implementation of chiral EFT considers nucleons and pions as degrees of freedom and a hierarchy of many-body forces specified by Weinberg power counting in which nucleon-nucleon inter- actions provide the most important contributions. Not explicitly resolved high-energy physics (such as ∆-excitations) is captured in low-energy constants (LECs) fit to experimental data. In recent years, significant progress has been made in deriving chiral forces to high orders (see Ref. [7] for a recent review). While consistent nucleon-nucleon (NN), three- (3N), and four-nucleon (4N) forces can in principle be applied, state-of-the-art calculations are still at N2LO or based on inconsistent N3LO NN plus N2LO 3N forces [1]. Apart from computational aspects, this is due to the lack of reliable fit values at N3LO for the two 3N LECs that contribute to the N2LO 3N forces. All other many-body contributions are predicted by pion-pion, pion-nucleon, and nucleon-nucleon physics. Meanwhile, the NN interactions have been worked out order-by-order up to N4LO [8–11] and even partly at N5LO [12]. Chiral interactions require moreover an ultraviolet regularization scheme associated with a cutoff scale, or possibly different scales in separate many-body sectors [9]. The broad freedom in choosing the functional form of the regulator has resulted in various new families of chiral NN potentials with conceptionally different regulators schemes: either in momentum [8, 11, 13–15] or coordinate space [16, 17], but also a mixed approach [9, 10] has been introduced recently. Although the regulator choice should be arbitrary, in practice, the actual form can have important consequences at a given order in the chiral expansion, referred to as regulator artifacts [18]. In addition to the development of new fit procedures, e.g., simultaneously optimizing NN plus 3N forces [13, 14], uncertainty quantification using statistical methods [15] (such as Bayesian analyses [19– 21]) has started to obtain attention in recent years. All of these new potentials combined with further advances in the many-body frameworks, renormalization group techniques, and the continuously increas- ing computational resources eventually led nuclear physics to an era of precision. Microscopic predictions of finite nuclei are the central goals of nuclear-structure studies [1]. The left panel in Fig. 1 compares ground-state energies and charge radii of several nuclei as obtained in ab initio calculations (dark-gray symbols) with experimental data (horizontal lines) [14]. In theory, the gen- eral issue is that many interactions overbind medium-mass to heavy nuclei and underestimate their radii [14, 15, 24]. Therefore, Ref. [14] included also heavier nuclei in the optimization protocol of the potential called NNLOsat. The agreement with experiment is manifested by the red diamonds in Fig. 1 (see also Refs. [4, 25, 26]). Nuclear-matter calculations moreover showed that NNLOsat has at same time real- istic saturation properties, especially, the saturation density is within the empirically-known range [14]. Saturation refers here to the empirical fact that both, the density as well as the binding energy of heavy 1 Introduction 7 4He 16O 22O 24O 36Ca 40Ca 48Ca 52Ca 54Ca 60Ca 48Ni 56Ni 68Ni 78Ni -9 -8 -7 -6 E/ A (M eV ) 2.0/2.0 (PWA) 2.2/2.0 (EM) 2.0/2.0 (EM) 1.8/2.0 (EM) AME 2012 " extrapol. Figure 1: Left: Ground-state energy per particle E/A (upper row) and computed charge radii relative to experiment∆rch (lower row) for selected nuclei based on chiral Hamiltonians (symbols). The experimental results are depicted by horizontal lines. The figure has been taken from Ref. [14]. Right: Ground-state energy per particle E/A up to heavier nuclei calculated in IM-SRG (extrapo- lated for 48,78Ni) based on four chiral Hamiltonians used in this work (colored symbols). Experi- mental data [22] is depicted by black symbols. The figure has been taken from Ref. [23] nuclei are approximately constant. This gives evidence that reproducing binding energies, charge radii, and saturation properties should be considered on an equal footing. Using the ab initio in-medium similarity renormalization group (IM-SRG), Ref. [23] recently studied binding energies (as well as charge radii) also for heavier systems based on four established chiral NN plus 3N Hamiltonians [27] used in this work. The results are shown in the right panel of Fig. 1. It is important to stress that the 3N forces were only fit to few-body data. Even though these interactions reproduce few-body data with same accuracy the situation can be different for medium-mass nuclei. To be specific, only the Hamiltonian “1.8/2.0” (purple symbols) agrees with experimental energies (black symbols), whereas the others follow the same trend but are shifted up to ∼ 1.5 MeV per particle [23]. The charge radii show a similar systematic behavior [23]. Qualitatively, these observations could be related to the Coester-like band of the corresponding saturation points, which we have calculated in Ref. [28]. These findings [23] indeed suggest that saturation properties are important for the accurate description of finite nuclei [14, 23, 29, 30]. Moreover they underline that saturation properties (e.g., from perturbative calculations) can be used as a guidance tool for constraining fits of next-generation chiral interactions. This required, however, a better understanding of the many-body convergence [31, 32] as well as the development of a novel efficient framework to be feasible in practice [32]. Homogeneous nuclear matter represents the theoretical testbed for benchmarking saturation properties. It is an idealized system obtained by extrapolating (heavy) nuclei to infinite mass number A and disregarding Coulomb interactions. Figure 2 depicts schematically the energy per particle E/A(β , n) of nuclear matter as a function of the total nucleon density n= np+nn. Here, np labels the proton and nn the neutron density. To measure the neutron excess of the system the asymmetry parameter β = 1− 2x can be introduced, or likewise the proton fraction x = np/n. Neutron matter (cyan line) is unbound and represents the neutron-rich extreme (β = 1) in Fig. 2. Adding protons (0 < β ¶ 1) lowers the energy per particle for (isospin-)asymmetry matter (green line) until bound (isospin-)symmetric nuclear matter (brown line) with equal nucleon densities (β = 0) is obtained. The energy difference between neutron and symmetric matter, S(n) = E A (β = 1, n)− E A (β = 0, n) (1.1) 8 1.1 Motivation symmetric matte r arbitrary proton fraction neutron matte r ~n2⁄3 15 MeV nn0 E A -16 MeV Figure 2: Illustration of the energy per particle E/A(β , n) of nuclear matter as a function of the asymme- try β and the total nucleon density n: neutron matter (β = 1, cyan line), isospin-asymmetric matter (0 < β ¶ 1, green line), and isospin-symmetric matter (β = 0, brown line). At nuclear saturation density n0 ' 0.16 fm−3, E/A(β = 0, n) has a minimum with respect to density and asymmetry. is called (nuclear) symmetry energy. We empirically know that the energy per particle in symmetric matter has a minimum with respect to n at so-called saturation density n0, which is associated with the saturation energy E/A(β = 0, n0) = −aV . The latter corresponds to the volume term in the Bethe- Weiszäcker mass formula [33] in absence of Coulomb interactions. Since symmetric matter minimizes in addition the energy with respect to β , one usually expands E/A(β , n) in a Taylor series about β = 0 (see, e.g., Refs. [34, 35]), E A (β , n) = E A (β = 0, n) + β2 Esym(n) +O � β4 � , (1.2) where we have neglected higher-order corrections in β4,6,... to good approximation [36] such that Esym(n) ' S(n). However, Refs. [35, 37] showed that these contributions should be investigated fur- ther. Only even powers in β contribute due to (approximate) isospin symmetry. Expanding also the density dependence of the coefficients in Eq. (1.2) about n0, E A (β = 0, n) = −aV + K 18 � n− n0 n0 �2 + . . . , and Esym(n) = Sv + L 3 � n− n0 n0 � + . . . , (1.3) imposes the incompressibility K , the slope parameter L of the symmetry energy, and the symmetry energy at saturation density Sv = Esym(n0). Within these truncations, only a few parameters, i.e., n0 ' 0.16 fm−3 (or ρ0 ' 2.7 · 1014 g cm−3), aV ' 16 MeV [38], K ' 240 MeV [39], Sv ' 32 MeV, and L ' 55 MeV [40], parametrize the nuclear-matter equation of state. Reproducing empirical ranges is a key benchmark for nuclear interactions. Figure 3 summarizes experimental constraints on the Sv − L correlation, which overlap in the rather small white area. Nevertheless, the strikingly wide ranges from different extractions indicate that especially L is not well constrained. In fact, the most precise constraints are due to the neutron-matter calculations by Hebeler et al. [41–43] based on many-body perturbation theory (MBPT) and chiral NN plus 3N interactions, and from quantum Monte Carlo calculations by Gandolfi et al. [44] for Hamiltonians adjusted to a fiducial Sv range. Additional constraints are therefore crucial [40, 45], e.g., 1.1 Motivation 9 Figure 3: Several experimental constraints on the symmetry energy Sv at saturation density and its slope parameter L: from nuclear masses [49], from analyses of isobaric analog states (IAS) [50], from the neutron skin thickness of tin (Sn) isotopes [51], from heavy-ion collisions (HIC) [52], from the dipole polarizability of 208Pb [53, 54], from giant dipole resonances (GDR) [55], and from modeling M − R relations (“astrophysics”) [56]. In addition, the two blue regions depict the- oretical constraints from neutron-matter calculations by Hebeler et al. (H) [41, 42] and Gan- dolfi et al. (G) [44]. The figure has been taken from Ref. [56] (see also Ref. [57]). from measurements of 48Ca nuclei [46, 47] that can also be reached in ab initio calculations [4]. On the theory side, this requires, in particular, an improved treatment of 3N forces in many-body calculations [28], reliable fit values for the 3N LECs at N3LO [32], and studies of asymmetric matter [28, 36, 37, 48]. Neutron stars [2, 58] are fascinating stellar objects. With masses of M ∼ 1.5 M� (solar masses) and radii of R ∼ 12 km, for instance, they are one of the densest forms of matter in the observable uni- verse [59]. Constraining the neutron-star equation of state as well as the M −R relation is thus a frontier in astrophysics [3, 60–64]. Two observations already demonstrated that their maximum mass is at least ∼ 2 M� [58, 65], whereas precise radius measurements have not been achieved yet [3]. Taking advan- tage of the lower maximum mass limit, nuclear-matter calculations based on chiral interactions lead to tight constraints [41, 42]. This is remarkable because typical central densities in neutron stars of sev- eral times saturation density exceed the low-energy region of chiral EFT. Therefore, large extrapolations in density are necessary to obtain the M − R relation as solution of the Tolman-Oppenheimer-Volkoff (TOV) equations. Causality and observational data are key constraints to rule out unphysical equations of state [42]. Interestingly, Ref. [66] found an empirical correlation between R and the pressure p(n) of neutron-star matter, as indicated in Fig. 4 for several equations of state. The typical proton fraction x ∼ 5 % in 10 1.1 Motivation – 44 – Fig. 3.— Empirical relation between pressure, in units of MeV fm−3, and R, in km, for EOSs listed in Table 1. The upper panel shows results for 1 M⊙ (gravitational mass) stars; the lower panel is for 1.4 M⊙ stars. The different symbols show values of RP−1/4 evaluated at three fiducial densities. Figure 4: Empirical relation between pressure p(n) and neutron-star radius R for several equations of state. The upper (lower) panel plots the approximate constant Rp−1/4 given a 1.0 M� (1.4 M�) neutron star and the densities n = n0, 1.5 n0, as well as 2 n0 (symbols). The figure has been taken from Ref. [66]. Note that the updated analysis in Ref. [57] also accounts for the observa- tion of the first 2 M� neutron star by Demorest et al. [65]. β-equilibrated matter is rather small. At saturation density, for a 1.4 M� neutron star, the (updated) empirical correlation reads [57] R(n0, 1.4 M�) = (9.52± 0.49) � p(n0) 1 MeV fm−3 � 1 4 km , (1.4) where p(n0) = n2 0 ∂nE/A(β ∼ 0.9, n)|n=n0 is mainly given by L. In consequence, radius predictions (for a given mass) are sensitive to uncertainties in this quantity, consistent to the narrow ranges in Fig. 3 from neutron-matter calculations and the corresponding constraints in Refs. [41, 42]. Using the empirical correlation (1.4) combined with tightly-constrained Sv and L, we also obtained a narrow range for the radius of 1.4 M� neutron star, 11.1 km¶ R1.4 M� ¶ 12.7 km [4], which underlines that pinning down the symmetry energy and its density dependence is highly-relevant for astrophysical applications. The present thesis addresses all of these exciting physics challenges. We initiate detailed asymmetric- matter studies using novel methods in order to extent previous state-of-the-art MBPT calculations to high orders in the chiral as well as the many-body expansion. To this end, we also focus on the improved treatment of 3N forces, particularly at N3LO, while making them accessible in general nuclear-matter considerations. We present several applications after assessing the many-body convergence. This naturally leads us to the construction of chiral Hamiltonians involving all many-body forces at N3LO with realistic saturation properties. For tight astrophysical predictions, we perform for the first time neutron- and symmetric-matter calculations up to fourth order in MBPT, enabled by the development of an efficient 1.1 Motivation 11 Monte-Carlo framework. Further, we provide quantitative insights into recent potentials within different regularization schemes. This paves the way for guiding optimizations of next-generation chiral potentials in terms of empirical saturation properties and triggers ab initio studies of finite nuclei based on these interactions. We organize the thesis as follows. In Secs. 1.2 and 1.3, we briefly introduce chiral EFT and set the stage for our improved nuclear-matter calculations. Section 2 reviews the operatorial definitions of the chiral many- body forces which are relevant for all of our calculations. In this section, we also discuss the details of the recent local, semilocal, and nonlocal chiral NN potentials, developed by Gezerlis, Tews et al. (2013 [16, 17]), Epelbaum, Krebs, and Meißner (2015) [9, 10], Entem, Machleidt, and Nosyk (2017) [11] as well as Carlsson, Ekström et al. (2016) [15], respectively. In Sec. 3, we perform a comprehensive Weinberg eigenvalue analysis based on this representative set of NN potentials. Our detailed comparison in terms of Weinberg eigenvalues provides crucial insights into their idiosyncrasies for different orders and partial waves. We develop in Sec. 4 an improved normal-ordering framework that allows to include general 3N inter- actions in nuclear-matter calculations, starting from a plane-wave partial-wave-decomposed form. We then make use of this method and include for the first time N3LO 3N interactions in MBPT up to third order and in self-consistent Green’s function (SCGF) theory. Using these two complementary many-body frameworks, we provide improved predictions for the equation of state of neutron matter at zero temper- ature and also analyze systematically the many-body convergence for different chiral EFT interactions. Furthermore, the normal-ordering framework is extended to finite temperatures. We present in Sec. 5 an efficient Monte-Carlo framework for perturbative calculations of nuclear mat- ter based on NN, 3N, and 4N interactions. The method enables the incorporation of all many-body contributions in a straightforward and transparent way. After first fourth-order calculations with chiral low-momentum interactions, we explore new chiral interactions up to N3LO. In Sec. 6, we discuss solutions of the BCS gap equation in the channels 1S0 and 3P2−3F2 in neutron matter based on the mentioned local and semilocal chiral NN interactions. In addition, we investigate for the first time the impact of N3LO 3N forces on pairing gaps and derive uncertainty estimates by taking into account results at different orders in the chiral expansion. We also discuss different methods for obtaining self-consistent solutions of the gap equation, including the first application of the modified Broyden method to the BCS gap equation. Finally, we conclude and give an outlook in Sec. 7. Throughout the thesis we use natural units, i.e., ħh= c = ħhc = kB = 1. 12 1.1 Motivation 1.2 Chiral effective field theory Nuclear interactions derived within chiral EFT are central for the present thesis. Starting with QCD we briefly review in this section the theoretical basics of chiral-symmetry breaking and its consequences for chiral EFT. We also introduce the hierarchy of chiral forces as well as theoretical uncertainty estimates based on the expected scaling behavior of chiral EFT. Comprehensive introductions can be found in Refs. [5–7, 67–70]. 1.2.1 Quantum chromodynamics and chiral symmetry The fundamental theory of strong interactions is quantum chromodynamics (QCD) [71]. It describes hadrons as composite particles made of elementary spin-1 2 quarks which interact via gluon exchanges. There are in total six known quark flavors as summarized in Table 1 including their charge, isospin, and approximate mass. In addition, (anti)quarks carry color charge, (anti)red, (anti)green, or (anti)blue, corresponding to the additive color model RGB, whereas gluons carry both, a color and an anticolor. The additional degree of freedom is introduced because isolated quarks have not been observed: they are con- fined to color-neutral hadrons. At high temperatures and/or high densities deconfinement sets in, forming a so-called quark-gluon plasma or a color-superconductor. Color-neutral hadrons consist of two quarks with color and anticolor (meson) or three different colors (baryon), but also exotic pentaquarks [72] and tetraquarks [73, 74] have been discovered recently at the Large Hadron Collider (LHC). An important property of QCD is that the strong coupling strength [75], αs(Q) = 6π 33− 2N f log−1 � Q ΛQCD � , (1.5) depends significantly on the momentum scale Q (see also Ref. [76]). It is therefore called running coupling. In Eq. (1.5) ΛQCD ∼ 200 MeV denotes the characteristic scale of QCD and N f is the number of active flavors. With increasing energy αs goes to zero such that QCD is perturbative at high-energy scales [77]. This is called asymptotic freedom [78, 79] and its discovery was recognized by a Nobel prize to Gross, Politzer, and Wilczek in 2004. However, at low-energy scales QCD is strongly nonperturbative (i.e., as ¦ 1), which makes direct applications to finite nuclei or nuclear matter very involved, if currently feasible at all. A promising approach to address these challenges is lattice QCD [80]: the space-time is discretized on a four-dimensional grid of finite volume, where quarks sit on the lattice points and gluons on the links between them. Observables are calculated from correlation functions after computing high-dimensional path integrals using Monte-Carlo integration. These extremely-demanding calculations are limited by computational resources. Effective field theories, nevertheless, can be used to extrapolate lattice results to the continuum limit and to physical quark masses. Vice versa, constraining LECs by lattice QCD would allow exciting first-principles studies of nuclear matter as well as finite nuclei (see also Ref. [81]). As preparation for chiral effective field theory, we discuss in the following the symmetries of QCD for the three lightest quark flavors (based on Refs. [6, 83]). The QCD Lagrangian is given by LQCD = ∑ f=u, d, s � q f i /D q f −m f q f q f � − 1 2 Tr GµνG µν , (1.6) with /D = γµDµ, the covariant derivative Dµ = ∂µ + i gsAµ, the strong coupling constant gs associated with αs = g2 s /(4π), the quark fields q f as well as masses m f , the gluon fields Aµ, and the gluon field strength Gµν. Let us first consider the chiral limit of vanishing quark masses, m f → 0. Introducing left- and 1.2 Chiral effective field theory 13 Table 1: Properties of the six quarks (for details see Ref. [82]). Up, down, and strange quarks are much lighter than the others. For example, the quark content of protons and neutrons is (uud) and (udd), respectively. flavor ( f ) charge in e isospin mass (m f ) in MeV up +2 3 +1 2 ≈ 2 down −1 3 −1 2 ≈ 5 strange −1 3 − ≈ 96 charm +2 3 − ≈ 1280 bottom −1 3 − ≈ 4180 top +2 3 − ≈ 173000 right-handed quark fields qL,i = 1/2(1−γ5)qi and qR,i = 1/2(1+γ5)qi, respectively, the Lagrangian (1.6) reads L 0 QCD = ∑ f=u, d, s � qL, f i /D qL, f + qR, f i /D qR, f � − 1 2 Tr GµνG µν , (1.7) which has a global U(3)L × U(3)R symmetry in flavor space, i.e., Eq. (1.7) is invariant under rotations of left- and right-handed quarks by independent unitary matrices. This group can also be written as SU(3)L × SU(3)R × U(1)A × U(1)V . The vector symmetry U(1)V corresponds to the baryon number (an exact symmetry also for nonzero quark masses), whereas the axial symmetry U(1)A is broken by the axial anomaly. Further, chiral symmetry SU(3)L×SU(3)R allows independent rotations of left- and right-handed quarks in terms of SU(3) matrices. For nonzero quark masses, however, it is explicitly broken due to the mass term in the Lagrangian (1.6), ∑ f m f q f q f = ∑ f qR, f m f qL, f + h.c. . (1.8) Assuming equal but nonzero quark masses SU(3)L × SU(3)R reduces to its subgroup SU(3)V . As can be seen in Table 1 the mass difference of up and down quarks is indeed small, leading to approximate isospin symmetry SU(2)V ⊂ SU(3)V . Chiral symmetry is in addition spontaneously broken, even in the chiral limit, as indicated by the absence of parity doublets in nature. For instance, the ρ-meson (J P = 1−) and its partner the a1-meson (J P = 1+) have the same quantum numbers except for parity but fairly different masses. Isospin symmetry, on the other hand, is approximately observed by comparable masses of the three charged states ρ0 and ρ± [6]. A spontaneously-broken continuous symmetry gives rise to a so-called Goldstone boson for each generator [84]. These particles are massless in case the broken symmetry was exact. Since chiral symmetry is also explicitly broken we therefore deal with so-called pseudo-Goldstone bosons that have finite mass: pions, kaons, and the η-meson. Pions are the lightest because of they have no strange-quark content. 1.2.2 Hierarchy of chiral many-body forces Chiral effective field theory (EFT) was pioneered by Steven Weinberg in a series of publications in the early 1990’s [85–87] (see also Ref. [88]) and has become since then the modern approach to nuclear forces [5, 6]. It considers nucleons as well as pions as the relevant degrees of freedom at the low- 14 1.2 Chiral effective field theory energy scales of nuclear physics, where the fundamental quarks and gluons are not resolved. Different choices may be efficient as well (or even more efficient) depending on the energy scale of interest, e.g., below the pion mass, also pions are integrated out, called pion-less EFT [89] (see Ref. [90] for a recent application). Another example is ∆-full EFT that is expected to improve the convergence of chiral EFT by using nucleons, pions, and ∆-resonances as degrees of freedom [5, 6]. Hence, EFTs inherently predict their own breakdown. For chiral EFT, the breakdown scale roughly corresponds to the mass of heavier mesons, in particular, the ρ-mesons: Λb ∼ 500 MeV < mρ is a more conservative value [9]. Above the breakdown scale additional degrees of freedom have to be added to the theory. Although not treated explicitly those contributions are still encoded in LECs fit experimental data. Using these effective degrees of freedom, one determines the most general Lagrangian, LEFT = Lππ +LπN +LNN + . . ., consistent with the exact and broken symmetries of QCD [88]. Note that this connection to the fundamental theory distinguishes chiral EFT from phenomenological approaches (e.g., Refs. [91–93]). In the Lagrangian,Lππ contains the dynamics between pions, LπN the pion-nucleon interactions, and LNN the nucleon-nucleon contact interactions, while the ellipsis accounts for contributions with two or more nucleons and pions, or A-nucleon contact interactions [6]. The infinitely many terms of the effective Lagrangian can be organized in powers of the expansion parameter, i.e., Qν with Q(p) =max � p Λb , mπ Λb � . (1.9) Equation (1.9) is the ratio of a typical nucleon momentum p or the pion mass mπ (soft scale) and Λb (hard scale). The value Q(p) ∼ 1/3 suggests that higher-order terms are sufficiently less important, if the coefficients of the expansion (including the LECs) are of natural size. In analogy to the multipole expansion in electrodynamics, the idea is to successively work out different orders Qν until the desired accuracy is obtained. This can be done because only finite number of terms and LECs contributes at each order (see below), being crucial for the predictive power of chiral EFT. Furthermore, the systematic expansion naturally allows to estimate the neglected contributions in form of theoretical uncertainties (see Sec. 1.2.3). One derives, practically, nuclear forces in a diagrammatic representation based on the effective Lagrangian at a given order. To associated these diagrams with their corresponding order ν a power-counting scheme is required. Weinberg power counting [85, 86] is based on dimensional analysis and obtains for a diagram involving N nucleons [5, 6] ν= −2+ 2N + 2 (L − C) + ∑ i ∆i , with ∆i = di + ni 2 − 2¾ 0 , (1.10) Here, L specifies the number of loops and C denotes the number of separately connected pieces. The sum goes over all vertices i, where di represents the number of derivatives or pion mass insertions and ni is the number of nucleon fields. Based on the power counting (1.10) Fig. 5 depicts the hierarchy of chiral forces up to N4LO. The orders ν= 0, 2, 3, . . . are referred to as leading order (LO), next-to-leading order (NLO), next-to-next-to-leading order (N2LO), respectively, and so on. Notice that ν= 1 is forbidden due to parity as well as time-reversal invariance [6]. We give in the following a brief overview of the chiral forces up to N4LO based on Refs. [6, 7]. Section 2 is then dedicated to a detailed discussion including the operatorial expressions relevant for our work, antisymmetrization, and regularization. One expects the dominant contributions in the chiral expansion from LO NN forces. They are given by momentum-independent (i.e., S-wave) contact interactions plus the long-range 1π-exchange potential. Nucleons (pions) are depicted in Fig. 5 by solid (dashes) lines. At NLO, momentum-dependent contacts (up to P-waves) as well as the 2π exchange for the intermediate-range attraction start contributing. Many-body forces are suppressed, VNN� V3N� V4N . . ., specifically, the first nonvanishing 3N forces appear at N2LO in three topologies (orange-shaded). The 3N 2π exchange is governed by the LECs c1, 3,4, similar to the 2π exchange in the NN forces at this order. Different determini- 1.2 Chiral effective field theory 15 N2LO (Q3) N3LO (Q4) N4LO (Q5) NLO (Q2) LO (Q0) 3N force 4N forceNN force Figure 5: Hierarchy of chiral nuclear forces up to N4LO. The diagrams are organized according to Wein- berg power counting (1.10). Dashed lines denote pions, solid lines denote nucleons. Solid dots, filled circles, filled rectangles, filled diamonds, and open rectangles depict vertices with ∆i = 0, 1, 2, 3, 4, respectively. Blue-, orange-, and green-shaded diagrams are available for cal- culations, whereas the others (gray-shaded) are currently under development (3N forces) or have not been worked out yet (4N forces). The figure has been modified from Ref. [94], see also Refs. [6, 69]. ations of these show significant deviations of the order of 30 % but agree within the uncertainties [70]. Additionally, the 3N LECs cD and cE emerge from the 1π-exchange-contact and the pure 3N contact in- teraction, respectively, which can be fit to few-body observables. There are no additional contacts due to N2LO NN forces. At N3LO, one has subleading 3N interactions as well as the first nonvanishing 4N forces (green-shaded). Both are free of unknown parameters. Moreover, contact interactions (up to D-waves), 2π exchanges, and the appearance of 3π exchanges in the NN forces complete this order. At N4LO, further corrections to these pion exchanges occur. The derivation of the 3N interactions at this order is currently ongoing [95, 96], while the 4N forces have not been worked out yet (gray-shaded). During the course of the thesis we study NN forces up to N4LO and perform nuclear-matter calculations with consistent many-body contributions up to N3LO. 1.2.3 Theoretical uncertainties As discussed in the previous section, the EFT framework allows to estimate theoretical uncertainties for an observable X (p) at a given momentum scale p [9]. In this work, X (p) corresponds to the energy per particle of nuclear matter as well as the pairing gap, associated with the Fermi momentum (or its average). Uncertainties arise in chiral EFT from several sources [9], e.g., from the truncation of the chiral expansion [9, 15, 19, 97], or from the determination of πN and contact LECs [15, 21, 98–100]. The cited references already indicate that nowadays advanced statistical techniques are being used to study 16 1.2 Chiral effective field theory theoretical uncertainties in more detail. We refer the reader to Ref. [20] for a comprehensive recipe for uncertainty quantification using Bayesian methods [101]. The employed framework to calculate the observable leads to additional uncertainties. These will be addressed in Sec. 1.3.1 for many-body perturbation theory. Cutoff variation in the regulator function (see Sec. 2.1.1) measures the sensitivity to neglected contact interactions and thus systematic uncertainties due to the truncation of the chiral expansion. However, the residual cutoff dependence can be quite misleading as the example in Ref. [9] demonstrates. Regarding NN forces, contacts appear only at even orders, so cutoff variation probes similarly N3LO short-range physics at both, NLO and N2LO. The uncertainty attached to the observable at NLO (or likewise N3LO) consequently will likely be underestimated. Furthermore, the range of applicable cutoff values excludes in practice too high momenta, whereas for low-momentum scales finite-cutoff effects distort extractions of realistic uncertainties. A new uncertainty estimate has recently been proposed in Refs. [9, 10] and applied to few-body calcu- lations [102–105] as well as to nuclear matter [32, 106]. The corresponding statistical interpretation is provided in Ref. [19] in terms of Bayesian uncertainty quantification. Instead of cutoff variations at fixed chiral order, the new approach considers results for X := X (p) at different orders while keeping cutoff values constant. This allows to assess the order-by-order convergence of the chiral expansion related to that observable. As we will report in Sec. 2.1.1, there are several complete sets of NN potentials available by now (e.g., from LO up to N4LO) that can be studied accordingly. To be specific, the prediction is given by the telescoping series X (ν) = ∑ν i=0 dX (i) at chiral order ν= 0, 2, 3, . . ., where dX (i) =    X (0) i = 0 , X (2) − X (0) i = 2 , X (i) − X (i−1) i ¾ 3 , (1.11) are the individual corrections from each order (see also Ref. [102]). Assuming the expected scaling in powers of the expansion parameter (1.9), i.e., dX (i) = O (Qi X (0)), Refs. [9, 10] defined the theoretical uncertainty as δX (ν) = ( Q2 � �X (0) � � ν= 0 , max 2¶ j¶ν � Qν+1 � �X (0) � � , Qν+1− j � �dX ( j) � � � ν¾ 2 . (1.12) That is, the range X (ν) ± δX (ν) specifies a band for the prediction. To furthermore ensure a systematic rate of convergence within the uncertainty estimate the additional constraint δX (ν) ¾max j,k �� �X ( j¾ν) − X (k¾ν) � � � (1.13) was imposed which accounts for higher-order results. Notice that the above discussion is based on consistent calculations at each order, however, it should be modified for ν¾ 3 if many-body contributions are neglected. In these cases, Ref. [102] avoids Eq. (1.13) and adopts in place of Eq. (1.12) the constraints δX (ν) =max � Qν+1 � �X (0) � � , Qν−1 � �dX (2) � � , Qν−2 � �dX (3) � � � , (1.14a) δX (2) ¾QδX (0) , and (1.14b) δX (ν¾3) ¾QδX (ν−1) . (1.14c) 1.2 Chiral effective field theory 17 1.3 Infinite nuclear matter From the perspective of infinite nuclear matter, we discuss in this section MBPT and the relationship to other nonperturbative methods. In particular, we elaborate on the inclusion of 3N contributions beyond Hartree-Fock using normal ordering. Furthermore, we provide all energy relations in MBPT that are relevant for our improved matter calculations. 1.3.1 Many-body perturbation theory Perturbation theory is a well-known framework in physics and at the heart of the present thesis. In particular, we are interested MBPT for infinite matter [18, 27, 28, 31, 32, 37, 48, 107–113] but the same concept also applies to finite nuclei [114, 115], see Refs. [116, 117] for recent applications. To familiarize ourselfs with notations, we briefly introduce the formalism at zero temperature, and refer to the standard literature [118–121] for extensive discussions. For now, we consider the generic nuclear Hamiltonian, H = T + V = H0 +λH1 , (1.15) where T denotes the kinetic energy operator, V contains the interactions, and λ is the expansion parameter. To apply perturbation theory, we have partitioned Eq. (1.15) in an unperturbated part H0 and a (small) perturbation λH1. Then, perturbation theory solves, if it converges, the Schrödinger equation of the many-body system, (H0 +λH1) |Ψi〉= Ei |Ψi〉 , with i ¾ 0 , (1.16) in terms of a perturbation expansion about the unperturbed system, i.e., Ei = ∑∞ n=0λ n E(n)i and similarly for the eigenstates |Ψi〉. It is more likely an approximation in practice since the series has to be truncated at a given finite order. We focus here on ground states, hence, i = 0 in Eq. (1.16). The coefficients E(n)i=0 are order-by-order determined by the known solutions of the unperturbed Schrödinger equation H0 |Φi〉 = E(0)i |Φi〉, as discussed in the following. More strictly, this corresponds to the diagonal case of perturbation theory, allowing simplifications, while generally it is only required that |Φ0〉 denotes an eigenstate in the arbitrary orthonormal basis |Φi¾0〉 [119]. Multiplying the Schrödinger equation (1.16) from left by 〈Φ0| and using the intermediate normalization 〈Φ0|Ψ0〉= 1 leads to ∆E := E0 − E(0)0 = 〈Φ0|λH1|Ψ0〉 . (1.17) It remains to derive an expression that relates |Ψ0〉 with the known |Φ0〉. To this end, we define the projection operator onto the unperturbed ground state P = |Φ0〉 〈Φ0| and its complement Q = 1− P, such that we obtain |Ψ0〉= (P +Q) |Ψ0〉= |Φ0〉+Q |Ψ0〉 . (1.18) In infinite nuclear matter |Φ0〉 represents the free Fermi sea, whereas |Φi¾1〉 includes excitations of particles and holes with respect to this ground state (see also the discussion in Sec. 1.3.3). Adding a term ζ |Ψ0〉 on both sides of the Eq. (1.16) as well as multiplying by Q results in, Q |Ψ0〉= Q ζ−H0 (λH1 − E0 + ζ) |Ψ0〉 , (1.19) 18 1.3 Infinite nuclear matter (a) (b) (c) (d) Figure 6: Hugenholtz diagrams for NN forces at second (a) and third order (b–d). Particles (holes) are indicated by up (down) arrows. The dots correspond to antisymmetrized interaction vertices. At third order, one has contributions from particle-particle (b), particle-hole (c), and hole-hole excitations (d). The figure has been modified from Ref. [122]. which is iterated in Eq. (1.18) to read |Ψ0〉= |Φ0〉+ Q ζ−H0 (λH1 − E0 + ζ) |Ψ0〉= ∞ ∑ n=0 � Q ζ−H0 (λH1 − E0 + ζ) �n |Φ0〉 . (1.20) We employ Rayleigh-Schrödinger perturbation theory, i.e., ζ := E(0)0 , and evaluate Eq. (1.17) ∆E = E0 − E(0)0 = ∞ ∑ n=0 〈Φ0|λH1 (R0 (λH1 −∆E))n |Φ0〉 , with R0 = Q E(0)0 −H0 = ∑ k 6=0 |Φk〉 〈Φk| E(0)0 − E(0)k . (1.21) Furthermore, we expand in a perturbation series, ∆E = ∑∞ n=1λ nE(n)0 . Organizing the terms in powers of λ→ 1 determines the desired coefficients of the expansion, here given up to fourth order [118], E(0)0 = 〈Φ0|H0|Φ0〉 , E(1)0 = 〈Φ0|H1|Φ0〉 , (1.22a) E(2)0 = 〈Φ0|H1R0H1|Φ0〉 , E(3)0 = 〈Φ0|H1R0(H1 − E(1)0 )R0H1|Φ0〉 , (1.22b) E(4)0 = 〈Φ0|H1R0(H1 − E(1)0 )R0(H1 − E(1)0 )R0H1|Φ0〉 − E(2)0 〈Φ0|H1R2 0H1|Φ0〉 . (1.22c) The presented analytic derivation of perturbation theory comes along with an equivalent pictorial ap- proach, so-called diagrammatic perturbation theory, which is conceptionally similar to Feynman tech- niques. One draws all possible (e.g., Hugenholtz) diagrams with n vertices (dots) and connects them by continuous lines, following remarkable simple rules [120, 122]. Because of Goldstone’s linked-diagram theorem [123] only connected diagrams contribute to the expansion. The explicit cancellation of size- inconsistent terms (see also Ref. [124]) related to unlinked diagrams transforms Rayleigh-Schrödinger perturbation theory to MBPT [118, 125]. In practice, the diagrammatic approach is more convienient and thus typically the method of choice, especially, regarding automation on a computer [122, 126–128]. Rules to translate diagrams to analytic expressions are given in the cited literature. It is however a nontrivial assumption that the perturbation series convergences at a useful rate. The efficiency of MBPT clearly depends on the underlying interaction as well as the chosen partition in Eq. (1.15). To quantify the perturbativeness of recent NN potentials in free space, we make use of Weinberg eigenvalues in Sec. 3 as a powerful diagnostic tool. Additionally, in Sec. 4.3 we benchmark the neutron- matter equation of state involving NN plus 3N interactions up to N3LO order-by-order to a nonperturbative method, where calculations using two partitions serve as a many-body uncertainty. 1.3 Infinite nuclear matter 19 One has for NN interactions 1, 3, 39, . . . Hugenholtz diagram(s) at second, third, fourth order, respectively, and so on [122, 129]. Figure 6 shows these at second and third order; for the fourth-order diagrams we refer to Ref. [130]. Below, we summarize the underlying analytic expressions up to fourth order, as this is the highest order considered in the present thesis. They are expressed in the particle-hole formalism, where only the particle (hole) states created above (below) the Fermi surface are considered (see, e.g., Ref. [118]). For brevity we define E(n)NN := E(n)0 . In the following Sec. 1.3.2, we will discuss normal-ordering and the inclusion of many-body forces at and beyond the Hartree-Fock level. Energy relations up to second order The energy contributions up to second order based on antisymmetrized NN interactionsA12VNN are given by the following expressions (see also Refs. [43, 118, 119]) T V = + ∑ i 〈i|T |i〉 , (1.23) E(1)NN V = + 1 2 ∑ i j 〈i j|A12VNN|i j〉 , (1.24) E(2)NN V = + 1 4 ∑ i j ab 〈i j|A12VNN|ab〉 〈ab|A12VNN|i j〉 Di jab . (1.25) Equation (1.25) is associated with diagram (a) in Fig. 6. We use the short-hand notation for the single- particle states |i〉 = |kiσiτi〉, having the momentum ki, the spin and isospin projections σi = ± 1 2 and τi = ± 1 2 , respectively. Furthermore, particles are labeled by a, b, . . . and holes by i, j, . . ., so the sums transform into, e.g., ∑ a −→ ∑ σaτa ∫ dka (2π)3 � 1− nτa ka � , and ∑ i −→ ∑ σiτi ∫ dki (2π)3 nτi ki , (1.26) where nτi ki is the Heaviside step function. Intermediate states beyond first order are weighted in terms of the single-particle energies, Di jk...abc... = εki + εk j + εkk + . . .− εka − εkb − εkc − . . . . (1.27) The trivial partition in Eq. (1.15), H0 = T and H1 = V , corresponds to the free spectrum with εki = k2 i /(2m), whereas H0 = T + VHF and H1 = V − VHF adds first-order self-energy corrections to the kinetic energy, called Hartree-Fock spectrum. Due to translational invariance, both, the kinetic-energy operator T and the Hartree-Fock potential VHF are diagonal in the plane-wave basis (see also Ref. [121]). We have thus the second-quantized Fock operator H0 = ∑ i j fi j a† j ai with fi j = δi j εi and the single-particle energies εi = k2 i 2m + ∑ j 〈i j |A12VNN | i j〉 . (1.28) In our calculations, we average Eq. (1.28) over spin as well as isospin quantum numbers. The two employed spectra lead to the same Hartree-Fock energy, i.e., the sum of all zero- and first-order terms. Adding 3N contributions to Eq. (1.28) will be addressed in Sec. 1.3.2. 20 1.3 Infinite nuclear matter For completeness, we note that there is an additional single-excitation diagram at second order, which is anomalous due to momentum conservation. It cancels using a Hartree-Fock spectrum [118]. Energy relations at third order One has in total three diagrams at third order when using a Hartree-Fock spectrum, depicted by (b–d) in Fig. 6. This results from a fine cancellation of several diagrams driven by the Hartree-Fock Hamiltonian. In a free spectrum, however, eleven additional (partly anomalous) terms arise [118]. We therefore consider solely Hartree-Fock single-particle energies beyond second order. Specifically, one has then contributions from hole-hole, particle-hole, and particle-particle excitations with respect to the Hartree-Fock reference state. These correspond to the following expressions, respectively, (see also Refs. [112, 118, 119]) E(3)1 V = + 1 8 ∑ i jkl ab 〈i j|A12VNN|ab〉 〈kl|A12VNN|i j〉 〈ab|A12VNN|kl〉 Di jabDklab , (1.29a) E(3)2 V = + ∑ i jk abc 〈i j|A12VNN|ab〉 〈ak|A12VNN|ic〉 〈bc|A12VNN| jk〉 Di jabDjkbc , (1.29b) E(3)3 V = + 1 8 ∑ i j abcd 〈i j|A12VNN|ab〉 〈ab|A12VNN|cd〉 〈cd|A12VNN|i j〉 Di jabDi jcd . (1.29c) Time reversal (exchanging holes and particles) relates the hole-hole and particle-particle terms, unfortu- nately, without reducing the number of diagrams to be computed. In total, the third-order contribution is given by E(3)NN V = E(3)1 V � � � � hh + E(3)2 V � � � � ph + E(3)3 V � � � � pp . (1.30) Energy relations at fourth order: an overview The fourth order consists of 39 linked diagrams when using a Hartree-Fock spectrum [118, 122]. They are categorized according to the level of excitations obtained after the second interaction with respect to the Fermi sea [131]. One has thus 4 single-, 12 double-, 16 triple-, and 7 quadruple-excitation diagrams, hence, the total energy of all fourth-order terms is given by E(4)NN V = 4 ∑ i=1 E(4)i V � � � � single + 16 ∑ i=5 E(4)i V � � � � double + 32 ∑ i=17 E(4)i V � � � � triple + 39 ∑ i=33 E(4)i V � � � � quadruple . (1.31) The sum of all contributions is order-by-order real in perturbation theory because the nuclear Hamiltonian is hermitian. This is also the case for each individual term up to third order. However, starting at fourth order we encounter complex-conjugated pairs of diagrams, which in total are again real [130]. Exploiting momentum conservation, in addition, reduces the number of diagrams to be computed to effectively 24. In the following, we carefully provide the complete set of energy expressions (including time-reversed pairs) as preparation for Sec. 5. These are more familiar in quantum chemistry but have not been studied to the best of our knowledge in infinite-matter calculations. 1.3 Infinite nuclear matter 21 Energy relations at fourth order: single excitations There are four contributions having single excitations after the second interaction [130, 131]: E(4)1 V = + 1 4 ∑ abcde i jk 〈i j|A12VNN|ab〉 〈ab|A12VNN|c j〉 〈ck|A12VNN|de〉 〈de|A12VNN|ik〉 Di jabDic Dikde , (1.32a) E(4)2 V = − 1 4 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈ab|A12VNN|c j〉 〈kl|A12VNN|id〉 〈cd|A12VNN|kl〉 Di jabDic Dklcd , (1.32b) E(4)3 V = − 1 4 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kb|A12VNN|i j〉 〈al|A12VNN|cd〉 〈cd|A12VNN|kl〉 Di jabDkaDklcd , (1.32c) E(4)4 V = + 1 4 ∑ abc i jklm 〈i j|A12VNN|ab〉 〈kb|A12VNN|i j〉 〈lm|A12VNN|kc〉 〈ac|A12VNN|lm〉 Di jabDkaDlmac , (1.32d) Note that these (anomalous) diagrams do not contribute at zero temperature because of momentum conservation or chancel in a Hartree-Fock spectrum. The diagrams labeled by (2, 3) are related via complex conjugation. Energy relations at fourth order: double excitations There are twelve contributions having double excitations after the second interaction [130, 131]: E(4)5 V = + 1 16 ∑ abcde f i j 〈i j|A12VNN|ab〉 〈ab|A12VNN|cd〉 〈cd|A12VNN|e f 〉 〈e f |A12VNN|i j〉 Di jabDi jcd Di je f , (1.33a) E(4)6 V = + 1 16 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈ab|A12VNN|cd〉 〈kl|A12VNN|i j〉 〈cd|A12VNN|kl〉 Di jabDi jcd Dklcd , (1.33b) E(4)7 V = + 1 16 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|i j〉 〈ab|A12VNN|cd〉 〈cd|A12VNN|kl〉 Di jabDklabDklcd , (1.33c) E(4)8 V = + 1 16 ∑ ab i jklmn 〈i j|A12VNN|ab〉 〈kl|A12VNN|i j〉 〈mn|A12VNN|kl〉 〈ab|A12VNN|mn〉 Di jabDklabDmnab , (1.33d) E(4)9 V = − 1 2 ∑ abcde i jk 〈i j|A12VNN|ab〉 〈ab|A12VNN|cd〉 〈kd|A12VNN|ie〉 〈ce|A12VNN|k j〉 Di jabDi jcd Djkce , (1.33e) E(4)10 V = − 1 2 ∑ abcde i jk 〈i j|A12VNN|ab〉 〈kb|A12VNN|ic〉 〈ac|A12VNN|de〉 〈de|A12VNN|k j〉 Di jabDjkac Djkde , (1.33f) E(4)11 V = − 1 2 ∑ abc i jklm 〈i j|A12VNN|ab〉 〈kl|A12VNN|i j〉 〈am|A12VNN|cl〉 〈cb|A12VNN|km〉 Di jabDklabDkmbc , (1.33g) 22 1.3 Infinite nuclear matter E(4)12 V = − 1 2 ∑ abc i jklm 〈i j|A12VNN|ab〉 〈ak|A12VNN|c j〉 〈lm|A12VNN|ik〉 〈cb|A12VNN|lm〉 Di jabDikbc Dlmbc , (1.33h) E(4)13 V = + ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈ak|A12VNN|c j〉 〈cl|A12VNN|dk〉 〈d b|A12VNN|il〉 Di jabDikbc Dil bd , (1.33i) E(4)14 V = − ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kb|A12VNN|c j〉 〈cl|A12VNN|id〉 〈ad|A12VNN|kl〉 Di jabDikac Dklad , (1.33j) E(4)15 V = − ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kb|A12VNN|c j〉 〈al|A12VNN|kd〉 〈cd|A12VNN|il〉 Di jabDikac Dilcd , (1.33k) E(4)16 V = + ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kb|A12VNN|ic〉 〈al|A12VNN|d j〉 〈dc|A12VNN|kl〉 Di jabDjkac Dklcd , (1.33l) The complex-conjugation pairs are labeled by (6, 7), (9, 10), and (11, 12). Energy relations at fourth order: triple excitations There are sixteen contributions having triple excitations after the second interaction [130, 131]: E(4)17 V = − 1 2 ∑ abcde i jk 〈i j|A12VNN|ab〉 〈ak|A12VNN|cd〉 〈cb|A12VNN|ek〉 〈ed|A12VNN|i j〉 Di jabDi jkbcd Di jde , (1.34a) E(4)18 V = − 1 2 ∑ abcde i jk 〈i j|A12VNN|ab〉 〈ak|A12VNN|cd〉 〈cd|A12VNN|e j〉 〈eb|A12VNN|ik〉 Di jabDi jkbcd Dikbe , (1.34b) E(4)19 V = − 1 2 ∑ abc i jklm 〈i j|A12VNN|ab〉 〈kl|A12VNN|ic〉 〈mb|A12VNN|kl〉 〈ac|A12VNN|mj〉 Di jabDjklabc Djmac , (1.34c) E(4)20 V = − 1 2 ∑ abc i jklm 〈i j|A12VNN|ab〉 〈kl|A12VNN|ic〉 〈mc|A12VNN|k j〉 〈ab|A12VNN|ml〉 Di jabDjklabc Dmlab , (1.34d) E(4)21 V = + ∑ abcde i jk 〈i j|A12VNN|ab〉 〈ak|A12VNN|cd〉 〈cb|A12VNN|e j〉 〈ed|A12VNN|ik〉 Di jabDi jkbcd Dikde , (1.34e) E(4)22 V = + 1 4 ∑ abcde i jk 〈i j|A12VNN|ab〉 〈ak|A12VNN|cd〉 〈cd|A12VNN|ek〉 〈eb|A12VNN|i j〉 Di jabDi jkbcd Di j be , (1.34f) E(4)23 V = + 1 4 ∑ abc i jklm 〈i j|A12VNN|ab〉 〈kl|A12VNN|ic〉 〈mc|A12VNN|kl〉 〈ab|A12VNN|mj〉 Di jabDjklabc Djmab , (1.34g) 1.3 Infinite nuclear matter 23 E(4)24 V = + ∑ abc i jklm 〈i j|A12VNN|ab〉 〈kl|A12VNN|ic〉 〈mb|A12VNN|k j〉 〈ac|A12VNN|ml〉 Di jabDjklabc Dlmac , (1.34h) E(4)25 V = − 1 4 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈ak|A12VNN|cd〉 〈l b|A12VNN|i j〉 〈cd|A12VNN|lk〉 Di jabDi jkbcd Dklcd , (1.34i) E(4)26 V = − ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈ak|A12VNN|cd〉 〈ld|A12VNN|ik〉 〈cb|A12VNN|l j〉 Di jabDi jkbcd Djl bc , (1.34j) E(4)27 V = − ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|ic〉 〈ac|A12VNN|dl〉 〈d b|A12VNN|k j〉 Di jabDjklabc Djkbd , (1.34k) E(4)28 V = − 1 4 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|ic〉 〈ab|A12VNN|d j〉 〈dc|A12VNN|kl〉 Di jabDjklabc Dklcd , (1.34l) E(4)29 V = + 1 2 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈ak|A12VNN|cd〉 〈l b|A12VNN|ik〉 〈cd|A12VNN|l j〉 Di jabDi jkbcd Djl cd , (1.34m) E(4)30 V = + 1 2 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|ic〉 〈ab|A12VNN|dl〉 〈dc|A12VNN|k j〉 Di jabDjklabc Djkcd , (1.34n) E(4)31 V = + 1 2 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|ic〉 〈ac|A12VNN|d j〉 〈d b|A12VNN|kl〉 Di jabDjklabc Dkl bd , (1.34o) E(4)32 V = + 1 2 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈ak|A12VNN|cd〉 〈ld|A12VNN|i j〉 〈cb|A12VNN|lk〉 Di jabDi jkbcd Dkl bc . (1.34p) Note that the (anomalous) diagrams with index 25 and 28 do not contribute at zero temperature because of momentum conservation or chancel in a Hartree-Fock spectrum. The complex-conjugation pairs are labeled by (25, 28), (26, 27), (29, 30), and (31, 32). Energy relations at fourth order: quadruple excitations There are seven contributions having quadruple excitations after the second interaction [130, 131]: E(4)33 V = + ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|cd〉 〈cb|A12VNN|il〉 〈ad|A12VNN|k j〉 Di jabDi jklabcd Djkad , (1.35a) E(4)34 V = + 1 16 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|cd〉 〈cd|A12VNN|i j〉 〈ab|A12VNN|kl〉 Di jabDi jklabcd Dklab , (1.35b) E(4)35 V = + 1 16 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|cd〉 〈ab|A12VNN|kl〉 〈cd|A12VNN|i j〉 Di jabDi jklabcd Di jcd , (1.35c) 24 1.3 Infinite nuclear matter E(4)36 V = − 1 4 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|cd〉 〈cb|A12VNN|i j〉 〈ad|A12VNN|kl〉 Di jabDi jklabcd Dklad , (1.35d) E(4)37 V = − 1 4 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|cd〉 〈ad|A12VNN|kl〉 〈cb|A12VNN|i j〉 Di jabDi jklabcd Di j bc , (1.35e) E(4)38 V = − 1 4 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|cd〉 〈cd|A12VNN|il〉 〈ab|A12VNN|k j〉 Di jabDi jklabcd Djkab , (1.35f) E(4)39 V = − 1 4 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|cd〉 〈ab|A12VNN|k j〉 〈cd|A12VNN|il〉 Di jabDi jklabcd Dilcd . (1.35g) These expressions can be cast into the following forms E(4)40 V = 1 2 � E(4)33 V + E(4)33 V � , E(4)41 V = E(4)34 V + E(4)35 V , (1.36a) E(4)42 V = E(4)36 V + E(4)37 V , E(4)43 V = E(4)38 V + E(4)39 V , (1.36b) where the energy denominator in E(4)40 is simplified compared to E(4)33 and only the contributions with index 40, 41, 42, and 43 have to be evaluated (four instead of seven). The corresponding analytic expressions read [131] E(4)40 V = + 1 2 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|cd〉 〈cb|A12VNN|il〉 〈ad|A12VNN|k j〉 Di jabDil bc Djkad , (1.37a) E(4)41 V = + 1 16 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|cd〉 〈cd|A12VNN|i j〉 〈ab|A12VNN|kl〉 Di jabDi jcd Dklab , (1.37b) E(4)42 V = − 1 4 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|cd〉 〈cb|A12VNN|i j〉 〈ad|A12VNN|kl〉 Di jabDi j bc Dklad , (1.37c) E(4)43 V = − 1 4 ∑ abcd i jkl 〈i j|A12VNN|ab〉 〈kl|A12VNN|cd〉 〈cd|A12VNN|il〉 〈ab|A12VNN|k j〉 Di jabDjkabDilcd . (1.37d) One proves the relations (1.36) by rewriting the sum of two denominators [131], e.g., δ−1 1 +δ −1 2 = � Di jabDi jcd Dklab �−1 , (1.38a) δ1 = Di jabDi jklabcd Di jcd , (1.38b) δ2 = Di jabDi jklabcd Dklab , (1.38c) combined with the permutation symmetries of Di jk...abc... and, for E(4)40 , one needs to interchange the tuples: (i, j), (k, l), (a, b), and (c, d). 1.3 Infinite nuclear matter 25 1.3.2 Normal ordering with 3N forces In this section, we extend the discussion of MBPT to include 3N forces. Their analytic definition up to N3LO will be given in Sec. 2.2. Normal ordering is a well-known method to account for dominant 3N con- tributions in terms of density-dependent effective two-body interactions (see, e.g., Refs. [28, 43, 132–134] for infinite matter and Refs. [24, 117, 135, 136] for nuclear structure). They are obtained by summing one particle over the occupied states of a finite-density reference state. However, we also consider here the computationally more involved residual 3N term of the normal-ordered Hamiltonian [137]. The new Monte-Carlo framework in Sec. 5 is well-suited for including this (and other) frequently neglected contributions. In second quantization, the general nuclear Hamiltonian with antisymmetrized NN and 3N interactions, H = ∑ 12 T12a† 1a2 + 1 (2!)2 ∑ 1234 〈12|A12V12|34〉 a† 1a† 2a4a3 + 1 (3!)2 ∑ 123456 〈123|A123V123|456〉 a† 1a† 2a† 3a6a5a4 , (1.39) is inherently normal ordered with respect to the vacuum state |0〉. That means, creation operators a† i are located on the left-hand side of annihilation operators ai, so the vacuum expectation value vanishes by construction, i.e., 〈0|a(†)|0〉 = 0. However, this is just a specific choice and, in fact, it is usually more convenient to consider likewise a finite-density reference state. As discussed in Sec. 1.3.1, in infinite matter it is natural to work with the Fermi sea in which all single-particle states |i〉= |kiσiτi〉 are filled up to the Fermi momentum kF, |Φ〉 := ki