FISSION IN EXOTIC NUCLEI USING DENSITY FUNCTIONAL THEORY By Zachary Matheson A DISSERTATION Michigan State University in partial fulfillment of the requirements Submitted to for the degree of Computational Mathematics, Science and Engineering — Dual Major Physics — Doctor of Philosophy 2019 FISSION IN EXOTIC NUCLEI USING DENSITY FUNCTIONAL THEORY ABSTRACT By Zachary Matheson Historically, most experimental and theoretical studies of fission have centered around the region of actinides near 235U, which includes isotopes of uranium, plutonium, and thorium relevant for reactor physics and stockpile stewardship. Isotopes in this region tend to fission asymmetrically, with the larger prefragment influenced by the doubly-magic shell structure of 132Sn and resulting in a heavy fragment distribution centered around 140Te. However, fission is a common decay mode of nuclei, both lighter and heavier than the actinides. Given the nuclear physics community’s interest in rare and exotic nuclei, we have applied nuclear density functional theory to study spontaneous fission primary fragment yields in exotic systems found in other regions of the nuclear chart, including neutron-deficient 178Pt, superheavy 294Og, and neutron-rich nuclei with relevance to the r-process such as 254Pu. Dedicated to Eli. iii ACKNOWLEDGMENTS I first and foremost wish to express my gratitude to Witek, for providing direction, perspec- tive, and so many wonderful opportunities. I consider myself very lucky to have been one of Witek’s graduate students. I would like to acknowledge my fission collaborators, Nicolas, Jhilam, and Samuel. Thank you for your patience and support. I can’t stress that enough. I have been fortunate to rub shoulders with several terrific young scientists as part of Witek’s research group, many of whom have helped me along the way. Here I would like to recognize Erik, Kevin, Bastian, Yannen, Chunli, Simin, Jimmy, Samuel, Nicolas, Rolo, Yue, Nobuo, Futoshi, Mao, Maxwell, Mengzhe, Tong, Leo, ... I am grateful for the many graduate students I’ve encountered along the way. The community of graduate students at the NSCL is, in my opinion, one of its greatest assets. Perhaps it’s unfair of me to single out any individual; nevertheless, there are a few who made specific, important contributions at critical moments of my PhD career. I would like to recognize: Amy, for going first; John, for thinking different; Terri, for being reliable; Juan, for being relatable; Wei Jia, for being empathetic; and Mao, for being sincere. Special shout-out to the members of the East Lansing University Ward. Double shout-out to Ethan for being such a good roommate. And finally, love and gratitude to my wife and family. Melon, Beaky, KJ+N8, Krystal, Xander, Ninrick, Rick Stick, Mom and Dad: Thank you, and I love you. iv TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Overview of fission theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Liquid drop model 1.1.2 Microscopic-macroscopic approach . . . . . . . . . . . . . . . . . . . 1.1.3 Self-consistent models and the supercomputing era . . . . . . . . . . 1.2 The scope of this dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 2 Describing fission using nuclear density functional theory . . . 2.1 Nuclear density functional theory . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Density functional theory . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1.1 Kinetic energy term . . . . . . . . . . . . . . . . . . . . . . 2.1.1.2 Coulomb interaction term . . . . . . . . . . . . . . . . . . . 2.1.1.3 Nuclear interaction term . . . . . . . . . . . . . . . . . . . 2.1.1.4 Pairing interaction term . . . . . . . . . . . . . . . . . . . . 2.1.2 Hartree-Fock-Bogoliubov method . . . . . . . . . . . . . . . . . . . . 2.1.2.1 Kohn-Sham and Hartree-Fock methods . . . . . . . . . . . 2.1.2.2 Bogoliubov transformation . . . . . . . . . . . . . . . . . . 2.1.2.3 HFB equations . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.3 Nucleon localization function . . . . . . . . . . . . . . . . . . . . . . 2.2 Microscopic description of nuclear fission . . . . . . . . . . . . . . . . . . . . 2.2.1 Potential energy surfaces . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Collective inertia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3 WKB approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.4 Langevin dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.5 Fragment identification via localization function . . . . . . . . . . . . 3.1 Calculating the potential energy surface Chapter 3 Numerical implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 PES Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Calculating the collective inertia . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Minimum action path . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Stochastic Langevin calculations . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 4 Two fission modes in 178Pt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Asymmetric fission in the region of 180Hg 4.2 Multimodal fission of 178Pt 1 1 2 2 3 4 6 7 8 9 10 10 11 12 12 13 14 16 17 18 20 21 22 23 25 25 26 26 28 29 30 30 31 v 4.3 The physical origin of fragment asymmetry in the neutron-deficient sub-lead region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 5 Cluster decay in 294Og . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Cluster emission of superheavy nuclei . . . . . . . . . . . . . . . . . . . . . . 5.2 Predicted spontaneous fission yields of 294Og . . . . . . . . . . . . . . . . . . 5.3 Fragment formation in 294Og . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Experimental search for cluster emission in 294Og . . . . . . . . . . . . . . . Chapter 6 Fission yields for r-process nuclei . . . . . . . . . . . . . . . . . . 6.1 The role of fission in the astrophysical r process . . . . . . . . . . . . . . . . 6.2 Fission fragment yields for r-process nuclei . . . . . . . . . . . . . . . . . . . 6.3 Future survey-level calculations for the r process . . . . . . . . . . . . . . . . Chapter 7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Improved model fidelity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 7.2 Computational tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 37 37 39 45 45 48 48 50 53 56 57 60 61 APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 Appendix A Prescription for fragment identification via localization function 64 Appendix B Finite-temperature ATDHFB collective inertia . . . . . . . . . 70 78 Appendix C List of my contributions . . . . . . . . . . . . . . . . . . . . . . . . BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 vi LIST OF TABLES Table 2.1: Table B.1: A list of collective coordinates used throughout this dissertation to generate potential energy surfaces. . . . . . . . . . . . . . . . . . . . As explained in the text, a cutoff for fb(cid:0)fa may be needed to reduce numerical noise in the collective inertia. This table shows how the calculated inertia changes depending on the size of the cutoff. . . . . 16 77 vii LIST OF FIGURES Figure 2.1: Proton and neutron densities (cid:26)p; (cid:26)n (in nucleons/fm3) and spatial lo- calizations Cp;Cn for 240Pu at several configurations along the most- probable path to fission. Figure from [1]. . . . . . . . . . . . . . . . Figure 2.2: Schematic overview of the Langevin framework for fission calcula- tions used in this dissertation. Figure from [2]. . . . . . . . . . . . . Figure 3.1: Mq20q20 calculated for some configuration of 240Pu as a function of finite-difference spacing (cid:14)q. . . . . . . . . . . . . . . . . . . . . . . . Figure 3.2: Norm of the difference matrix between subsequent iterations of the density. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 4.1: Fragment yields as a function of N and Z for several nuclei ranging from actinides, where primary fission yields tend to be asymmetric, down to near-thorium, where yields become more symmetric, and finally to the region near neutron deficient 180Hg, where asymmetry . . . . . . . . . . . . . . . . . . . . . . . . returns. Figure from [3]. Figure 4.2: After projecting the fission fragment mass vs total kinetic energy (TKE) correlation (b) onto the TKE axis, the TKE distribution is deconvoluted into high- and low-energy components, centered around the values TKEhigh and TKElow (a). The mass distributions at these particular energies are fitted to a double (c) and single (d) Gaussian, respectively. This procedure is repeated for three different compound . . . . . . . . nucleus excitation energies E* (e-g). Figure from [4]. Figure 4.3: UNEDF1HFB potential energy surface for 178Pt as a function of mul- tipole moments Q20, representing elongation, and Q30, representing mass asymmetry. The static pathway is marked in purple. Note the two different trajectories ABCD and ABcd and their corresponding . . . . . . . . . . . . . . . . . . . . . . . . . NLFs. Figure from [4]. Figure 4.4: a) D1S potential energy surface for 178Pt as a function of multi- pole moments Q20, representing elongation, and Q30, representing mass asymmetry. The static pathway is marked in red. b) With the quadrupole moment fixed to Q20 = 190b, a PES is constructed as a function of Q40, which relates to the size of the neck, and Q30. This PES shows the existence of two symmetric modes, separated by a saddle point with a height of (cid:24)4 MeV. Figure from [4]. . . . . . . . 18 19 27 27 31 32 33 34 viii Figure 5.1: Calculated (left) and experimental (right) decay modes for heavy and superheavy elements. The theoretical estimates are based on an analysis of lifetimes calculated via empirical formulae. The boxed isotopes in the left panel are those which have been measured exper- Isotopes falling inside the 1(cid:22)s contour are predicted to imentally. live longer than 1(cid:22)s. Figure adapted from [5]. . . . . . . . . . . . . Figure 5.2: Dominant decay modes for superheavy elements are indicated based on predictions obtained in a nuclear DFT-based framework are in- dicated. The label “RS” stands for “reflection symmetric”, meaning the least-action wavefunction at the saddle is reflection symmetric, and “NRS” stands for “non-reflection symmetric.” The legend indi- cates predicted half-lives on a logarithmic scale. Figure from [6]. . . Figure 5.3: Comparison of the PESs for 294Og in the (Q20; Q30) collective plane obtained in UNEDF1HFB (a), D1S (b), and SkM* (c) EDFs. The ground-state energy Egs is normalized to zero. The dotted line in each figure corresponds to E0 (cid:0) Egs = 1 MeV, which was used to de- termine the inner and outer turning points. The local energy minima . . . . . at large deformations are marked by stars. Figure from [7]. 38 38 40 Figure 5.4: Fission fragment distributions for 294Og obtained in UNEDF1HFB (a), D1S (b), and SkM* (c) EDFs using the non-perturbative crank- ing ATDHFB inertia and the baseline dissipation tensor (cid:17)0 = ((cid:17)11; (cid:17)22; (cid:17)12) = (50(cid:22)h; 40(cid:22)h; 5(cid:22)h). Known isotopes are marked in gray [8]. Magic num- bers 50, 82, and 126 are indicated by dotted lines. Figure from [7]. . 43 Figure 5.5: Upper panel: Predicted heavy fragment mass (left) and charge (right) yields of 294Og using different functionals (top, linear scale). Lower panels: collective inertias (middle, logarithmic scale) and dissipa- tion tensor strengths (bottom, logarithmic scale). The baseline cal- culation was performed using the UNEDF1HFB functional in a 4D space with non-perturbative cranking ATDHFB inertia and dissipa- tion tensor strength (cid:17)0. Figure from [7]. . . . . . . . . . . . . . . . . 44 Figure 5.6: Nucleon localization functions for several deformed configurations of 294Og. For comparison, localizations are also shown for the fragments 208Pb and 86Kr on the left side of each subplot. The configurations shown correspond to Fig. 1 from the paper, with multipole moments 2 ) ((cid:24)outer turning line); (Q20; Q30) = (a) (75 b, 0); (b) (120 b, 17 b (c) (140 b, 24 b 3 2 ); and (d) (264 b, 60 b 3 3 2 ) ((cid:24)scission). Figure from [7]. 46 ix Figure 6.1: Figure 6.2: Schematic overview of the r process. The competition between neu- tron capture and beta decay increases the mass of a seed nucleus until it reaches the region where fission dominates (spontaneous, beta-delayed, or neutron-induced). . . . . . . . . . . . . . . . . . . . Final abundance patterns from four different simulations of neutron star merger ejecta, each using a different model for the fission frag- ment yields. The results are presented in two graphs solely for visual clarity. The dots represent the known solar r-process abundance pattern. Figure from [9]. . . . . . . . . . . . . . . . . . . . . . . . . Figure 6.3: Localization function, PES, and yield distribution for the sponta- neous fission of 254Pu. The stars in the PES indicate outer turning points which were used to compute the cumulative yields. . . . . . . Figure 6.4: Similar to Figure 6.3, but for 290Fm. . . . . . . . . . . . . . . . . . Figure 6.5: Heatmap showing isotopes whose fission yields are especially relevant to the r process. This combines the results from four different nuclear structure models, and counts how many of the models found each isotope to have an integrated fission flow above a certain threshold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure from [10]. Figure A.1: Localization function for 240Pu with lines drawn to represent pre- fragment identification . . . . . . . . . . . . . . . . . . . . . . . . . Figure A.2: PES and yield distribution for spontaneous fission of 240Pu. On the left, the stars in the PES indicate outer turning points which were used to compute the cumulative yields. On the right, the black curve in the yields corresponds to the Langevin result, and experimental data comes from [11, 12]. . . . . . . . . . . . . . . . . . . . . . . . . Figure A.3: Similar to Figure A.2 but for 294Og. . . . . . . . . . . . . . . . . . . 49 50 52 54 55 66 67 69 x Chapter 1 Introduction 1.1 Overview of fission theory Nuclear fission is the fundamental physical process by which a heavy nucleus splits into two smaller nuclei. Fission was first observed by Hahn and Straßmann in 1939 [13] when they bombarded uranium atoms with neutrons and detected barium, a paradoxically lighter element. At the time of their discovery, they were unable to explain this puzzling outcome. An explanation came shortly thereafter in letters to Nature by Meitner and Frisch [14] and by Bohr [15]. Meitner and Frisch’s described the findings in the following way: On account of their close packing and strong energy exchange, the particles in a heavy nucleus would be expected to move in a collective way which has some resemblance to the movement of a liquid drop. If the movement is made sufficiently violent by adding energy, such a drop may divide itself into two smaller drops… It seems therefore possible that the uranium nucleus has only small stability of form, and may, after neutron capture, divide itself into two nuclei of roughly equal size (the precise ratio of sizes depending on finer structural features and perhaps partly on chance). In addition to neutron-induced fission, spontaneous fission, which is a different form of fission so-dubbed because it occurs without bombardment by neutrons or any other projectiles, was reported by Flerov and Petrjak in a single-paragraph letter to Physical Review in 1940 [16]. For the remainder of this dissertation, I will be referring mainly to spontaneous fission unless stated otherwise. 1 Following the explanation of Meitner and Frisch, fission is relatively simple to concep- tualize, but remarkably difficult to explain quantitatively. Quantitative fission predictions based on microscopic nuclear theory are challenging because of the large number of particles involved, along with the complex many-body dynamics, which takes place when the nu- cleus deforms and splits. Historically, one could argue that theoretical attempts to describe nuclear fission have leapt forward in three waves. 1.1.1 Liquid drop model The first attempt to describe fission goes back to the very beginning of the nuclear age in the 1930s with the liquid drop model. The liquid drop model was first developed by Weizsäcker in 1935 [17] as a way to describe collective properties of nuclei. It was later adapted by Bohr and Wheeler to quantitatively describe nuclear fission in terms of bulk properties of nuclei [18]. This model was able to successfully explain nuclear binding energies and the energetics of nuclear fission. The liquid drop model has its weaknesses, however. For example, it could not explain the fission fragment mass asymmetry which characterizes spontaneous fission in many actinides. Indeed, until the 1960s, little attention was given to the quantum nature of the fissioning nucleus. 1.1.2 Microscopic-macroscopic approach A breakthrough came during a time of renewed interest in fission, triggered by the discovery of fission isomerism in 1962 by Polikanov, et al [19]. This was understood as a manifesta- tion of nuclear shape deformation [20, 21]. Recognizing the important connection between 2 shell effects, collective shape deformations, and the fission process, Strutinsky added a quan- tum mechanical correction to the liquid drop energy in 1967 in order to account for the added stability that occurs when a nucleus contains a “magic number” of protons and/or neutrons [22, 23, 24]. The models based on such an approach go by the name “microscopic- macroscopic” because they combine the “macroscopic” bulk properties of the liquid-drop model with the “microscopic” quantum mechanical Strutinsky shell correction. Microscopic-macroscopic (“micmac”) fission models are computationally fairly inexpen- sive, and can achieve quite satisfactory results (some recent highlights include Refs. [25, 26, 27, 28, 29]). In this approach, for instance, the fragment mass asymmetry that occurs in actinides is understood in terms of strong shell effects within the fragments during the split. However, since the model is based on a phenomenological description of what is actually a quantum many-body system, its predictive power is limited, with no clear way of mak- ing systematic improvements. Micmac models also rely on many parameters which must be adjusted to data, and are therefore subject to all the pitfalls that come with parameter fitting. A more fundamental approach would be to consider individual nucleon states using a consistent quantum many-body method. 1.1.3 Self-consistent models and the supercomputing era The third phase of fission theory is taking place now, heralded by the age of supercomputers. In fact, fission was listed as an exascale problem in a 2017 technical report to the Department of Energy [30] - that is, one of the problems which motivates the drive towards exascale computing. For large systems with many particles, density functional theory (DFT) is a way to recast the (cid:24)200 particle Schrödinger equation into a simpler problem involving only a few local densities and currents (see Section 2.1.1 as well as [31]). State-of-the-art 3 calculations such as the ones described in this dissertation and others reviewed in [32] allow us to combine nuclear DFT techniques with modern high-performance computing platforms to predict fission properties, such as lifetimes and fragment yields. These advances in computing come simultaneously with technological advances which allow experimental nuclear physics to reach far beyond what has been achieved previously. New facilities such as the upcoming Facility for Rare Isotope Beams (FRIB) at Michigan State University, which is projected to nearly double the number of isotopes which can be produced synthetically [33], and the Superheavy Elements Factory in Dubna, Russia [34], which focuses on the synthesis of new superheavy elements and detailed study of known superheavy isotopes, will provide crucial missing information. Many current facilities will soon receive upgrades which should allow them to produce experimental fission data that is much more detailed and precise than ever before [3]. Together, state-of-the-art experi- ments and high-performance computing are expected to lead to rapid advancement in our understanding of atomic nuclei and their decays, including fission. 1.2 The scope of this dissertation Microscopic models of fission (as self-consistent models are often called) are increasingly able to predict properties of fission fragments; however, a comprehensive microscopic description of fission fragments (including mass and charge distributions, excitation and kinetic energy distributions, angular dependence, spin, neutron emission) still remains elusive. Chapter 2 will discuss our approach for estimating the mass and charge of fission fragments. The next challenge is to do these calculations efficiently. In every theoretical calculation, one must ask the questions “What approximations can I safely make?” and “What are 4 the important degrees of freedom for this problem, and which ones can I ignore?” These simplifications, together with improvements to the computational workflow itself, such as better file handling and parallelization, can significantly reduce the total time-to-answer. Computational details will be discussed in Chapter 3. After describing the model in Chapters 2 and 3, it is applied to several isotopes in Chap- ters 4- 6 to compute primary fragment yields. Historically, most experimental and theoreti- cal studies of fission have centered around the region of actinides near 235U, which includes isotopes of uranium, plutonium, and thorium relevant for reactor physics and stockpile stew- ardship. Isotopes in this region tend to fission asymmetrically, with the larger prefragment influenced by the doubly-magic shell structure of 132Sn and resulting in a heavy fragment distribution centered around 140Te. However, given the aforementioned interest in rare and exotic nuclei, we have applied our model to study fragment distributions in exotic systems found in other regions of the nuclear chart. First, in Chapter 4 we discuss bimodal fission of the neutron-deficient isotope 178Pt, which until recently was expected to fission sym- metrically. This region is a good one in which to test fission models because of the small isospin asymmetry (N=Z (cid:25) 1:3 in this region, compared to N=Z (cid:25) 1:5 near the valley of stability). Then in Chapter 5 we discuss cluster radioactivity in 294Og, the heaviest element ever produced in a laboratory. In Chapter 6 we move to the neutron-rich side of the nuclear chart (N=Z > 1:7) to study isotopes which are expected to play a significant role in the astrophysical r process. This dissertation concludes in Chapter 7 with a discussion of our results and their signif- icance. Suggestions are then made for future model developments, computational improve- ments, and physical applications. 5 Chapter 2 Describing fission using nuclear density functional theory There are basically two microscopic frameworks in which to study fission: time-dependent and static (time-independent). Since fission is an inherently time-dependent process, time- dependent methods offer deep insight into the fission process and the characteristics of the fragments, especially kinetic and excitation energies [35, 36, 37, 38, 39]. However, they can only treat a single event at a time and are quite expensive, making them impractical for fission yield predictions. Furthermore, and most important to a discussion of spontaneous fission, is the fact that fully time-dependent approaches have no mechanism with which to de- scribe quantum tunneling, making them totally unsuited to spontaneous fission calculations. Finally, even if the tunneling problem was solved, this approach would still be unsuitable for (cid:0)22 any nucleus with a reasonably-long lifetime (compared to a typical time step size (cid:24)10 sec). In the other, so-called static approach, it is assumed that collective motion of the nucleus is slow compared to the intrinsic motion of the nucleons, and therefore that collective and intrinsic degrees of freedom can be decoupled. This assumption of adiabaticity is supported by experimental evidence which suggests a characteristic timescale for fission (from saddle to (cid:0)18sec, compared to typical nuclear scission, the point at which the neck snaps) of (cid:24)10 (cid:0)20(cid:0)10 6 (cid:0)22sec (see [40] for a review of fission timescale experiments). timescales on the order of (cid:24)10 The validity of the adiabatic assumption was further discussed for fission and other nuclear processes in [41]. The assumption of adiabaticity justifies the creation of a potential energy surface (PES) in some space of collective shape coordinates, and the dynamics of fission are then described in terms of trajectories across the PES. Quantum tunneling pops out in a fairly natural way in this formalism using the WKB approximation, and half-life estimates follow. The kinetic and excitation energies can be computed in this framework, but they are extremely sensitive to the characterization of scission, which is not well-defined in static frameworks [42]. However, as we shall show, the static approach is well-suited to estimating fission yields. In an effort to be as self-consistent as possible, the PES is computed in the framework of nuclear DFT, which combines the Hartree-Fock-Bogoliubov variational approximation to the energy with a many-body method inspired by Kohn-Sham DFT. An overview of this self-consistent mean-field framework is described below, followed by a description of the dynamical calculations which are used to calculate fission yields. 2.1 Nuclear density functional theory As with all nonrelativistic quantum systems, nuclei can be described using the many-body Schrödinger equation. However, one often finds this type of description difficult or impossible in practice, for two reasons: • In order to use the Schrödinger equation, one needs to describe the interactions between nucleons. However, protons and neutrons interact via the strong nuclear force, which is non-perturbative at low energies and has a complex spin and isospin dependence. 7 • Even if an interaction was known precisely, nuclei are large systems made up of many protons and neutrons. Solving the Schrödinger equation directly quickly becomes computationally intractable as the number of particles increases. Nuclear density functional theory is one of several methods which has been developed to address these challenges. It is particularly useful for heavy nuclei, where other approaches such as configuration interaction and ab initio methods become prohibitively expensive. 2.1.1 Density functional theory Nuclear DFT is rooted in the Hohenberg-Kohn theorems [43], which are briefly described in the following. The first Hohenberg-Kohn theorem states that the energy of the system is a uniquely- defined functional E((cid:26)) of the particle density (cid:26)(⃗r). This means one can describe a compli- cated system of N interacting particles using a single density which depends on 3 coordinates, rather than many-body wavefunction depending on 3N coordinates - a huge simplification! The second Hohenberg-Kohn theorem states that the functional which gives the energy of the system will give the ground state energy if, and only if, it acts on the true ground state density: E((cid:26)) = Egs () (cid:26) = (cid:26)gs. Thus, given a particular functional E((cid:26)), one can vary the input density (cid:26) to minimize the total energy and be assured that one is approaching the ground state energy of the system. In the case of electronic systems, a variational prescription to compute the electron density was subsequently proposed by Kohn and Sham in [44]. However, extensions of this prescription to nuclei, which are self-bound objects governed by a poorly-known Hamiltonian, led to a Kohn-Sham scheme that was too difficult to be applied [45, 46, 47]. Instead, the practical alternative is to introduce mean-field 8 variational methods such as Hartree-Fock (HF) and Hartree-Fock-Bogoliubov (HFB), as we do in Section 2.1.2. The next step is to find the correct energy density functional. Unfortunately, neither the Hohenberg-Kohn theorems nor the Kohn-Sham method specify how this is to be done. In anticipation of mean-field plus pairing methods such as HFB (see Section 2.1.2), the total energy will instead be assumed to be a sum of several contributions, each of which can be treated individually in terms of the particle density (cid:26)(⃗r; ⃗r ′ ): E((cid:26)) = Ekin + ECoul + Enuc + Epair; (2.1) where Ekin is the kinetic energy term, ECoul comes from the Coulomb interaction between protons, Enuc is a nucleon-nucleon interaction term, and Epair describes the tendency of nucleons to form pairs, a behavior which would otherwise be smeared out in non-interacting mean-field models. Each of the terms in (2.1) is described below. 2.1.1.1 Kinetic energy term Defining the kinetic density (cid:28)(cid:11) = ∇ (cid:1) ∇′ ( tion is 1 (cid:0) 1 A Ekin = (cid:22)h2 2m (cid:12)(cid:12) ′ )∫ ( ) d3⃗r (cid:28)n(⃗r) + (cid:28)p(⃗r) ; (cid:26)(cid:11)(⃗r; ⃗r ) ⃗r=⃗r′, (cid:11) = p; n, the kinetic energy contribu- ( 1 (cid:0) 1 A ) (2.2) term where A is the number of nucleons and m is the mass of a single nucleon. The is a simple center-of-mass correction. 9 2.1.1.2 Coulomb interaction term The repulsive Coulomb interaction between protons is divided into a direct term and an exchange term, which is related to the antisymmetry of proton wave functions: ECoul = ECoul;dir + ECoul;exch; ∫ ∫ ECoul;dir = ECoul;exch = e2 2 e2 2 d3⃗r1d3⃗r2 d3⃗r1d3⃗r2 (cid:26)p(⃗r1)(cid:26)p(⃗r2) j⃗r1 (cid:0) ⃗r2j (cid:26)p(⃗r2; ⃗r1)(cid:26)p(⃗r1; ⃗r2) ; j⃗r1 (cid:0) ⃗r2j (2.3) (2.4) (2.5) : The direct term uses local densities, which are related to the nonlocal densities found in the exchange term: (cid:26)(⃗r) = (cid:26)(⃗r; ⃗r). Often, the exchange term is computed in the Slater approximation [48, 49]: ∫ ) 1 3 ( 3 (cid:25) d3⃗r(cid:26) 4 3 p (⃗r): ECoul;exch (cid:25) (cid:0) 3e2 4 (2.6) 2.1.1.3 Nuclear interaction term Describing the interaction energy between nucleons Enuc (and to a lesser extent, Epair) continues to be an active area of research in nuclear theory today [50, 51, 52, 53, 54]. For nuclear DFT, the most commonly-used strong-interaction energy density functionals belong to the Skyrme, Gogny, or covariant family of functionals (each of which is discussed in [31]). We will primarily be using Skyrme functionals, which can be written as a sum of time-even 10 and time-odd terms [55]: ∫ ESkyrme = d3⃗r ) ; Heven t + Hodd t ∑ ( t=0;1 t = C(cid:26) Heven Hodd t = Cs t (cid:26)2 t + C∆(cid:26) t (cid:26)t∆(cid:26)t + C(cid:28) t ⃗s2 t + C∆s t ⃗st∆⃗st + CT t + C t J2 t (cid:26)t(cid:28)t + CJ t ⃗st (cid:1) ⃗Tt + Cj t j2 t + C ∇J (cid:26)t∇ (cid:1) ⃗Jt; t ∇j t ⃗st (cid:1) (∇ (cid:2) ⃗jt); (2.7) (2.8) (2.9) where (cid:28)t is the kinetic energy density; Jt and ⃗J(cid:20);t are, respectively, the scalar and vector parts of the spin current density; ⃗st is the spin density; ⃗Tt is the spin kinetic density; and ⃗jt is the momentum density (to see how these each relate to (cid:26), see e.g. [31]). The index t = 0(1) refers to isoscalar(isovector) energy densities, e.g. (cid:26)0 = (cid:26)n + (cid:26)p ((cid:26)1 = (cid:26)n (cid:0) (cid:26)p). Note that Heven depends only on time-even densities (and similarly for Hodd ). t t Since this energy density functional is phenomenological, rooted in a zero-range contact force between nucleons, the coefficients must be determined from experiment and/or ab initio theory. There are dozens of Skyrme parameterizations on the market, each one optimized to a particular observable or set of observables. The parameter sets SkM* [56] and UNEDF1 [57] (along with its derivative, UNEDF1HFB [58]) have been optimized to datasets which include fission data, making them suitable for fission calculations. 2.1.1.4 Pairing interaction term The simplest mean field approximation fails to take into account some correlations between nucleons, especially correlations between nucleons which occupy nearby states. Such nucleons (for example, those occupying orbitals with equivalent orbital quantum numbers but opposite spins) have a tendency to form pairs, similar in mechanism to BCS superconductivity or 3He superfluidity [59]. To help account for these correlations, we introduce a so-called pairing 11 density (cid:20)(⃗r) = (cid:20)(r; (cid:27); (cid:28) ) (which will be defined in Section 2.1.2) and use a density-dependent pairing functional [60]: ∑ ∫ ) ) (cid:12) ( ( 1 (cid:0) (cid:26)(⃗r) (cid:26)0 Epair = V(cid:11) d3⃗r f ((cid:20)(⃗r)) ; (2.10) (cid:11)=p;n ∑ (cid:3) (cid:27)(cid:28) (cid:27)2(cid:20)(cid:11)(r; (cid:27)(cid:28) ; r;(cid:0)(cid:27); (cid:28) )(cid:20) (cid:11)(r; (cid:27)(cid:28) ; r;(cid:0)(cid:27); (cid:28) ) connects states with opposite spins. As with the nuclear interaction term, the pairing term contains adjustable parameters where f ((cid:20)(⃗r)) = V0; (cid:12); and (cid:26)0. The denominator (cid:26)0 determines whether pairing is concentrated within the volume ((cid:26)0 ! 1), around the surface ((cid:26)0 (cid:25) 0:16 fm(cid:0)3), or somewhere in between. The ex- ponent (cid:11) can partially account for the formation of surface effects, such as neutron skins and halos, but is usually set to (cid:11) = 1. The pairing strength V0 may be adjusted to experimental odd-even mass differences. 2.1.2 Hartree-Fock-Bogoliubov method 2.1.2.1 Kohn-Sham and Hartree-Fock methods Having now chosen an energy density functional, the variational principle is invoked to find the density which minimizes the total energy, approximating the ground state of the system. The Kohn-Sham variational method was proposed soon following the Hohenberg- Kohn theorems in [44]. Given the exact energy density functional, Kohn-Sham actually solves the many-body problem exactly; however, as stated before, its use in nuclear physics is limited, principally because nuclei are self-bound systems with no external constraining potential. Kohn-Sham has been extended to self-bound systems (for instance in [45, 46, 47]), but the resulting expressions are complicated and in practice it is common to fall back on a mean-field approximation, such as Hartree-Fock. 12 The Hartree-Fock method, which actually predates Kohn-Sham by over 30 years, treats the system as a set of independent particles, each interacting with a mean field created by the other particles [61]. The mean-field approximation is well-justified in nuclear physics, its primary experimental evidence being the existence of magic numbers. The downside to Hartree-Fock is that some correlations, such as pairing correlations, are lost. The Hartree- Fock-Bogoliubov method introduces pairing correlations in the following way. 2.1.2.2 Bogoliubov transformation In anticipation of the HFB formalism, we define the so-called Bogoliubov transformation. The fundamental entities in the Bogoliubov transformed basis are ‘quasiparticle’ states, defined by quasiparticle creation and annihilation operators given by i i or in block matrix notation, 0B@ (cid:12) y (cid:12) 1CA = 0B@ U y 1CA y V 0B@ c V T U T c y 1CA (cid:17) Wy 0B@ c y c y where c i ; ci are single particle creation and annihilation operators, and where the trans- y obey the fermion commutation formation matrix W must be unitary to ensure that (cid:12); (cid:12) relations [61]. The HFB ground state j(cid:8)0⟩ is defined to be a quasiparticle vacuum, such that (cid:12)(cid:22) j(cid:8)0⟩ = 0 for all (cid:22). As before, the system is described using a density matrix, which can be written in terms 13 ∑ ∑ i ∑ ∑ i (cid:3) i(cid:22)c y i ; V Vi(cid:22)ci; (cid:3) i(cid:22)ci + U y i + Ui(cid:22)c (cid:12)(cid:22) = y (cid:22) = (cid:12) (2.11) (2.12) (2.13) 1CA ; y of single particle operators and an HFB vacuum state: (cid:26)ij = ⟨(cid:8)0jc jcij(cid:8)0⟩. In the HFB case, one must also define a second density matrix, called the pairing density: (cid:20)ij = ⟨(cid:8)0jcjcij(cid:8)0⟩. In coordinate space, the density matrix (cid:26) and the pair tensor (cid:20) take the form ′ ′ y ) = ⟨(cid:8)0jc ⃗r′c⃗rj(cid:8)0⟩ ; ) = ⟨(cid:8)0jc⃗r′c⃗rj(cid:8)0⟩ : (cid:26)(⃗r; ⃗r (cid:20)(⃗r; ⃗r R = (cid:20) (cid:3) 1 (cid:0) (cid:26) The densities (cid:26) and (cid:20) are combined to form a single HFB generalized density matrix: 0B@ (cid:26) (cid:3) (cid:0)(cid:20) 1CA : (2.14) (2.15) (2.16) Note that in nuclei, there are two R’s: one for describing neutrons and another for protons. 2.1.2.3 HFB equations The ground state configuration of the system with a particular energy density functional is described by the density R which minimizes E(R). The solution can be found through the variational principle, in which the energy is minimized with respect to the generalized density, subject to the constraint that R2 = R. Defining the HFB Hamiltonian density Hba (cid:17) 2@E=@Rab, this variation leads to the expression [H;R] = 0, which is called the Hartree-Fock-Bogoliubov equation. The HFB equation is not typically solved in its commutator form, but rather it is recast in the following way: Recalling that two Hermitian operators whose commutator is zero can be simultaneously diagonalized, we choose to diagonalize H using the same Bogoliubov 14 transformation W which diagonalizes R: WyHW (cid:17) E or HW = WE; where E = 0B@ E(cid:22) 0 0 (cid:0)E(cid:22) 1CA (2.17) (2.18) is a matrix of quasiparticle energies. In this form, the problem can then be solved iteratively: 1. Choose some ansatz to estimate the density; 2. Construct the HFB Hamiltonian density H; 3. Solve the eigenvalue problem: HW = WE; 4. Extract the densities corresponding to the eigenfunctions (which are related to W); 5. Update H; 6. Repeat starting from 3 until some predetermined convergence criterion is met. Often one will want to minimize the energy of the system subject to a particular con- straint. Constraints can be introduced via the method of Lagrange multipliers. Some com- mon examples of constraints include multipole moments representing nuclear shape defor- mations. In this case, a particular multipole moment might be constrained to the value (cid:22)Q(cid:21)(cid:22): ∑ (⟨ ⟩ ) ′ E = E (cid:0) C(cid:21)(cid:22) ^Q(cid:21)(cid:22) (cid:0) (cid:22)Q(cid:21)(cid:22) 2 : (2.19) (cid:21)(cid:22) Since the Bogoliubov transformation breaks particle number symmetry, an important con- straint on average particle numbers represents the first step towards particle number restora- 15 Coordinate Units Description Axial quadrupole moment Roughly corresponds to elongation Triaxial quadrupole moment Axial octupole moment Roughly corresponds to mass asymmetry Dynamic pairing fluctuations b b b 3 2 MeV Q20 Q22 Q30 (cid:21)2 Table 2.1: A list of collective coordinates used throughout this dissertation to generate potential energy surfaces. tion: ⟨ ^Nn ′ E = E (cid:0) (cid:21)n where (cid:21)(cid:11) is determined by the condition that ⟩ (cid:0) (cid:21)p ⟩ ⟨ ^N(cid:11) ⟨ ⟩ ^Np ; (2.20) = N(cid:11) (more sophisticated particle number restoration schemes also exist, such as the Lipkin-Nogami method [62, 63, 64, 65, 66]). The collective coordinates employed in this dissertation are described in Table 2.1. 2.1.3 Nucleon localization function The major advantage nuclear DFT offers over microscopic-macroscopic models is a direct link to the underlying many-body structure. This can be visualized using the nucleon localization function (NLF), originally introduced for the identification of localized electronic groups in atomic and molecular systems in [67] and then applied to nuclei in [68, 1]. The NLF is defined using the single particle densities (cid:26); (cid:28), and j, and is related to the probability of finding two nucleons of the same spin (cid:27) and isospin q at locations ⃗r and ⃗r′: Pq(cid:27)(⃗r; ⃗r′) = (cid:26)q(cid:27)(⃗r)(cid:26)q(cid:27)(⃗r′) (cid:0) : (2.21) (cid:12)(cid:12)(cid:12)(cid:26)q(cid:27)(⃗r; ⃗r′) (cid:12)(cid:12)(cid:12)2 Suppose we fix the location ⃗r of one nucleon, and consider only a small spherical region 16 of radius (cid:14) around ⃗r. Then we can expand the conditional probability Rq(cid:27) = Pq(cid:27)=(cid:26)q(cid:27)(⃗r) in terms of (cid:14) to find [67]: Rq(cid:27)(⃗r; (cid:14)) (cid:25) 1 3 ( (cid:28)q(cid:27) (cid:0) 1 4 j∇(cid:26)q(cid:27)j2 (cid:26)q(cid:27) (cid:0) j2 q(cid:27) (cid:26)q(cid:27) ) (cid:14)2 + O((cid:14)3); (2.22) The term in parentheses is neither dimensionless nor normalized, so we divide by the Thomas- Fermi kinetic density (cid:28) T F 5 q(cid:27) to convert the measure to a proper scale. Finally, 3 note that Rq(cid:27) goes to zero whenever two like-spin nucleons are well-separated. In order to 2 3 (cid:26) q(cid:27) = 3 5(6(cid:25)2) have a measure which takes on its maximal value whenever a nucleon is well-localized, the following expression defines the NLF: ( 241 + Cq(cid:27) = ) 2 35(cid:0)1 : (2.23) (cid:28)q(cid:27)(cid:26)q(cid:27) (cid:0) 1 j∇(cid:26)q(cid:27)j2 (cid:0) j2 q(cid:27) 4 (cid:26)q(cid:27)(cid:28) T F q(cid:27) A localization value C (cid:25) 1 means that nucleons are well-localized; that is, the probability of finding two nucleons of equal spin and isospin in the same spatial region is low. A value of C = 1 2 corresponds to a Fermi gas of nucleons, as found in nuclear matter. As demonstrated in Figure 2.1, the NLF tends to sharpen the internal structure of the nucleus compared to the density. When applied to fission as in [69], it enables one to see the formation of well-defined prefragments whose shell structure is responsible for the peak of the fragment distribution. 2.2 Microscopic description of nuclear fission With a nuclear many-body method now in hand, it can be used as a tool to describe fission. Recently in [2], a nuclear DFT-based approach based on this assumption was used to compute 17 Figure 2.1: Proton and neutron densities (cid:26)p; (cid:26)n (in nucleons/fm3) and spatial localizations Cp;Cn for 240Pu at several configurations along the most-probable path to fission. Figure from [1]. fragment yields from a PES that was computed self-consistently, using the WKB approxi- mation to describe tunneling and Langevin dynamics to describe post-tunneling dissipation. The half-life can be computed as in [70]. The entire procedure, illustrated schematically in Figure 2.2, is described below. 2.2.1 Potential energy surfaces The primary degrees of freedom in the adiabatic approximation are nuclear collective shapes, and the basic ingredient to fission calculations is a PES. One could describe any nuclear con- figuration using a sufficiently-detailed PES; however, for practical computations one must use a truncated set of only a few collective coordinates. Thus, an important challenge for researchers is to select the most relevant collective coordinates, ideally while demon- 18 -10010(30,0)ρn(89,0)(185,18)(245,27)(336,41)-10010Cn↑-10010ρp-505-10010Cp↑-505-505-505-50500.040.0800.250.5000.030.0600.250.50r(fm)z(fm) Figure 2.2: Schematic overview of the Langevin framework for fission calculations used in this dissertation. Figure from [2]. strating that others can be safely neglected. Often, one will use the first few lowest-order moments, but there is some discussion within the community about whether these are the best for describing configurations on the way to fission, especially near scission (see, for in- stance, [71]). Other collective coordinates which are sometimes used in the literature include (cid:0) (⃗r(cid:0)⃗rneck)2 a2 qN = e for the range parameter), and non-geometric variables which control pairing and pairing , which is related to the size of the neck (a = 1 fm is a common value ⟨ fluctuations, such as the pairing gap ∆ or the mean value of particle number fluctuations ⟩ ∆ ^N 2 Once the appropriate collective-coordinate constraints ⃗q (cid:17) (q1; q2; : : : ) are chosen, the PES is computed on a mesh: one DFT calculation per grid point. The density at each point . ⃗q is then used to compute the HFB energy E ′ (⃗q) (2.19). 19 2.2.2 Collective inertia Just as important as potential energy to the fission dynamics is the collective inertia, which describes the tendency of the system to resist configuration changes (such as shape changes) [20]. The form of the collective inertia used here is the non-perturbative adiabatic time-dependent HFB (ATDHFB) inertia with cranking [72], which takes the form 1A ; (2.24) 0@@R21 M(cid:22)(cid:23) = (cid:22)h2 2 1 (Ea + Eb) @R12 (0);ba @q(cid:23) @R12 (0);ab @q(cid:22) @R21 (0);ba @q(cid:22) + (0);ab @q(cid:22) The subscripts and superscripts are explained in the full temperature-dependent derivation of the collective inertia found in Appendix B, but an important feature to note is that computing the inertia requires differentiating the density matrix with respect to a pair of collective coordinates. In expressions where the collective coordinates are shape degrees of freedom, the collective inertia acts as a masslike term to resist changes to the collective motion. A perturbative expression for the ATDHFB inertia also exists, which allows one to esti- mate the inertia without taking derivatives of the density [72]. It is computationally much faster and easier to implement, but it loses many of the important features of the inertia compared to the non-perturbative case, as we shall see in Chapter 5. Nevertheless, it is commonly-used in calculations and later on we shall show how this choice affects fission yields compared to the non-perturbative inertia. Another common expression for the collective inertia comes from the Generator Coordi- nate Method (GCM). The GCM inertia also exists in two varieties: perturbative and non- perturbative [73]. Like the ATDHFB inertia, the perturbative GCM inertia is smoothed-out compared to the non-perturbative inertia. Both the perturbative and non-perturbative GCM 20 inertias are found to be smaller in magnitude than their ATDHFB counterparts [74, 73]. 2.2.3 WKB approximation Spontaneous fission is an extreme case of quantum tunneling. If the wave function of the fissioning nucleus is assumed to be slowly varying inside the (large) potential barrier (which condition is satisfied under the adiabatic assumption), then the WKB approximation allows us to estimate the probability of tunneling through a classically-forbidden region in the PES. sin in the col- Under the WKB approximation, the most-probable tunneling path L(s)jsout lective space is found via minimization of the collective action ∫ √ 2M(s) (E′(s) (cid:0) E0)ds; S(L) = 1 (cid:22)h sout sin (2.25) where s is the curvilinear coordinate along the path L, M(s) is the effective collective inertia given by [70] M(s) = M(cid:22)(cid:23) dq(cid:22) ds dq(cid:23) ds (2.26) ∑ (cid:22)(cid:23) ′ and E (s) is the potential energy along L(s). E0 stands for the collective ground-state energy. The calculation is repeated for different outer turning points, and each of these points is then assigned an exit probability P (sout) = [1 + expf2S(L)g] (cid:0)1 [75]. The half-life corresponds to the pathway of minimum action, and the expression for the half-life is T1=2 = ln(2)=nP (smin), where the parameter n is the number of assaults (cid:25) on the fission barrier per unit time and which we set to the standard value n = 1MeV (cid:0)1 [75]. However, because of the exponential dependence of the half-life on the action, half-life calculations suffer from large uncertainties due to the uncertainties in theoretical 1020:38s h 21 ingredients. 2.2.4 Langevin dynamics Although the link between collective and intrinsic degrees of freedom was assumed away in the adiabatic approximation, it is necessary to reintroduce some connection between the two, especially in the region between the outer turning points and scission where the adiabatic assumption is weakest. Through the semiclassical Langevin approach [76], we introduce a ∑ dissipation tensor (cid:17) that mimics energy exchange between intrinsic and collective modes. The dissipation tensor is related to a random force with strength g via the fluctuation- dissipation theorem [77, 76]: k gikgjk = (cid:17)ijkBT. The fluctuation-dissipation theorem √ effectively couples collective and intrinsic degrees of freedom via the system temperature T, E(cid:3)=a where a = A=10MeV(cid:0)1 parameterizes the level density and the ij pipj (see [78, 79, 2] given by kBT = excitation energy is given by E ∑(M(cid:0)1 (sout) (cid:0) E (⃗q) = E ′ (⃗q) (cid:0) 1 2 ) (cid:3) ′ for details). After emerging from the classically-forbidden region of the PES, fission trajectories pro- ceed across the PES from the outer turning line, evolving according to the Langevin equa- tions: ( M(cid:0)1 ) jk dpi dt = (cid:0)pjpk 2 @ @qi ( M(cid:0)1 ) pk + gij(cid:0)j(t) ; jk ′ ) (cid:0) (cid:17)ij ( (cid:0) @E @qi M(cid:0)1 pj ; ij dqi dt = (2.27) (2.28) where pi is the collective momentum conjugate to qi, and (cid:0)j(t) is a Gaussian-distributed, time-dependent stochastic variable. From here one can see that the dissipation term (cid:0)(cid:17)ij jk pk represents a transfer of energy from the collective motion of the system to (M(cid:0)1 ) 22 intrinsic single particle degrees of freedom, while the random fluctuation term gij(cid:0)j(t) does the opposite. Dissipation is treated in our work phenomenologically, since a self-consistent description of dissipation has not yet been fully developed [80]. In the meantime, we use the values of the dissipation tensor from [2], which were fitted to reproduce the 240Pu spontaneous fission fragment yields. 2.2.5 Fragment identification via localization function The Langevin approach represents the best of what static approaches currently have to offer in terms of quantitative reproduction of the yields; however, it can be expensive and time- consuming to perform a full Langevin calculation with non-perturbative collective inertia in a multidimensional collective space. For applications where the level of quantitative accuracy is somewhat flexible (an example would be r-process network calculations, which are explained in Chapter 6), one might ask whether it is necessary to do all these calculations if simple predictions are required (e.g., most-probable fission yields), or if there is some kind of physical insight that can be leveraged to reduce the workload. By studying nucleon localization functions (Section 2.1.3) for deformed configurations of fissioning nuclei, we observe that, in many cases, two distinct fragments start to form independently and then separate well before scission is reached (see Figure 2.1, taken from [1], as well as Sections 4.2 and 5.3). Thus, a fissioning nucleus can be thought of as two well- formed prefragments connected by a neck. At scission, the two prefragments separate and the remaining neck nucleons are distributed to either of the two prefragments in some predictable way. Based on this, we developed a prescription for identifying prefragments and then estimating fragment yields at the outer turning line, using a statistical technique based on 23 the grand canonical ensemble. This prescription is described in Appendix A. 24 Chapter 3 Numerical implementation In this modern era of fission theory which is dependent upon high-performance computers, it is worth discussing some of the numerical and computational issues associated with com- puting the potential energy surface, the collective inertia, and the collective dynamics of a fissioning nucleus. 3.1 Calculating the potential energy surface By far, the most time-consuming part of a microscopic fission calculation is computing the PES. For this we use a pair of DFT solvers, HFODD [81] and HFBTHO [82], which solve the HFB equations in a deformed harmonic oscillator basis. The solver HFBTHO is limited to shapes with axial symmetry, while HFODD allows for the breaking of any symmetry needed (for instance, rotational invariance or parity). Since the major bottleneck of each of these programs involves constructing matrices which represent nucleonic densities and currents and then diagonalizing the HFB Hamiltonian matrix, the flexibility of symmetry breaking drastically increases the time-to-solution. Thus, HFBTHO is used for fast calculations, while HFODD is slower but more flexible. Fortunately, the problem of calculating a PES is embarrassingly-parallel. So while an in- dividual point in the PES may be time-consuming to compute, many points can be computed on a large platform with many processors simultaneously. This does have its limitations; 25 highly-deformed configurations may be numerically-unstable or ill-defined. In such cases, one fix can sometimes be to use a nearby point which converged successfully as a seed function. 3.1.1 PES Tools Apart from the issue of walltime, generating a PES creates a lot of data, which quickly becomes unwieldy. To help manage all this data, a Python module called PES Tools was developed for manipulating, extracting, interpolating, and plotting [83]. Aside from the parser scripts, which collect data from the output of a particular DFT solver and store it in an XML data structure, PES Tools is not dependent on a specific DFT solver. In particular, a submodule was created to interface between PES Tools and Fission Tools [84], a set of codes used for fission calculations. These codes are described in the following sections. 3.2 Calculating the collective inertia The partial derivatives from equation (2.24) are computed using the Lagrange three-point formula [72]: ( ) @R @q ′ (cid:0)(cid:14)q (cid:14)q ((cid:14)q + (cid:14)q′) R(q0 (cid:0) (cid:14)q) + ′ (cid:14)q (cid:0) (cid:14)q (cid:14)q(cid:14)q′ R(q0) + (cid:14)q (cid:14)q′ ((cid:14)q + (cid:14)q′) R(q0 + (cid:14)q ′ ): (3.1) (cid:25) q=q0 The accuracy and precision of the collective inertia M are therefore functions of the spacings ′, and of R. An accurate value of the collective inertia is especially important for (cid:14)q and (cid:14)q computing half-lives, for which there is an exponential dependence on M. Finite-difference calculations such as (3.1) are dependent on the size of the finite-difference 26 Figure 3.1: Mq20q20 calculated for some configuration of 240Pu as a function of finite- difference spacing (cid:14)q. Figure 3.2: Norm of the difference matrix between subsequent iterations of the density. ′ spacing (cid:14)q. Figure 3.1 shows the effect of different values of (cid:14)q(= (cid:14)q ) on the collective inertia for some configuration of 240Pu. After some initial fluctuations for large (cid:14)q, the value of M appears to be nearly constant over a range of several values of (cid:14)q. The abrupt jump at the end of Figure 3.1 is related to the level of convergence of the density matrix R, illustrated in Figure 3.2. Figure 3.2 shows the norm of the matrix which corresponds to the difference between the density matrix at the last iteration and the second-to-last iteration. This can be thought of 27 as numerical noise. Predictably, the noise decreases as the convergence parameter becomes tighter. This gives a sense of the uncertainty (cid:14)R associated with R, which in turn should be propagated through equations (3.1) and (2.24). Additionally, if the amount of noise in R is much larger than, or comparable in size to (cid:14)q, then M will suffer from jumps like (cid:0)6 b in Figure 3.1. To avoid such problems in our calculations, our the one near (cid:14)q (cid:25) 10 finite-difference spacing (cid:14)q should be larger than our HFB convergence parameter. There are additional complications which arise in the finite-temperature formalism. These are discussed in Appendix B. 3.3 Minimum action path For the tunneling part of the collective evolution described in Section 2.2.3, the dynamic programming method [85] was used to minimize the action. The dynamic programming method proceeds inductively: once the minimum action is known for all grid points up to a certain value of Q20, say qn layer with Q20 = qn+1 20 n and each grid point in layer n + 1, and then selecting the path which minimizes the total 20, then the minimum action at each grid point in the (n + 1)th is obtained by computing the action between each grid point in layer action at each grid point in layer n + 1: ( (cid:12)(cid:12)(cid:12) )(cid:12)(cid:12)(cid:12) Smin( ⃗Q) Q20=qn+1 20 = argmin ⃗Q′ Smin( ⃗Q′) + ∆Smin( ⃗Q; ⃗Q′) ′ 20=qn Q 20;Q20=qn+1 20 : (3.2) The inductive step which connects layers n and n + 1 involves several small, independent calculations which lend themselves well to shared-memory parallelism. This was implemented in the code using OpenMP, resulting in a walltime reduction from order O(nD) to O(nD=m), 28 where m is the number of processors (with the usual caveats that parallelization requires additional overhead time, and that the actual speedup might be somewhat less when the number of processors m reaches the same order as the number of points per layer, nD(cid:0)1). In a 2D calculation, where the runtime is on the order of seconds, the difference is inconsequential. However, this speedup was essential for the analysis of a 4D PES, as described in Chapter 5. 3.4 Stochastic Langevin calculations Once the action and relative probability are known for a set of points along the outer turning line, Langevin trajectories are computed starting from each outer turning point. These are straightforwardly evaluated at discrete time steps over a discretized PES mesh. Because of the random force term in equation (2.27), a large number of trajectories per outer turning point must be computed to reduce statistical uncertainty. Fortunately, each trajectory is completely independent of every other trajectory, lending the code readily to shared-memory parallelism. However, for the Fission Tools Langevin code, distributed- memory parallelism was chosen over shared-memory parallelism in order to simplify access to shared resources, such as variables and output files. 29 Chapter 4 Two fission modes in 178Pt 4.1 Asymmetric fission in the region of 180Hg As mentioned in the introduction, fission has been studied most carefully in the region of the actinides (Z=90 to Z=103) [86], since naturally-occurring isotopes in this region are fissile. Within this region, there is a characteristic tendency for fission fragment yields to be asymmetric (that is, one light fragment and one heavy fragment), with the heavy peak centered around A (cid:25) 140 [87]. This has been understood as a manifestation of nuclear shell structure in the prefragments: doubly-magic 132Sn drives the nucleus towards scission, and once the neck nucleons are divided up between the two fragments, the heavy fragment dis- tribution peaks near A=140 [88]. As one moves to the lighter nuclei, however, this tendency becomes less and less pronounced as yields tend to become more symmetric [89, 90]. For sub-thorium isotopes, there is no doubly-magic nucleus candidate that could drive the sys- tem toward asymmetry as there is with actinides. To the contrary: neutron-deficient 180Hg, for instance, might be expected to split evenly into two 90Zr fragments, each with N = 50 closed neutron shells. However, it was reported in a 2010 study [91] that neutron-deficient 180Tl undergoes beta-delayed fission, leading to the excited intermediate state 180Hg which in turn decays into two fragments of unequal mass. This finding triggered a flurry of theoretical papers trying to explain this new and unexpected phenomenon (for instance, see [6, 92, 93, 94]). 30 Figure 4.1: Fragment yields as a function of N and Z for several nuclei ranging from ac- tinides, where primary fission yields tend to be asymmetric, down to near-thorium, where yields become more symmetric, and finally to the region near neutron deficient 180Hg, where asymmetry returns. Figure from [3]. Meanwhile, a follow-up study using 178Tl [95] further established that this was a region of asymmetric fission, and not just a one-time occurrence. Since then, other nuclei in the region have been studied using a variety of reactions and techniques, each with mass asymmetric fragments. An overview of fission fragment distributions in the region of 180Hg which have been experimentally studied (as of 2016), marked by the experimental technique used, is shown in Figure 4.1. 4.2 Multimodal fission of 178Pt One particular follow-up experiment was performed to investigate the fission of 178Pt [4], which differs from 180Hg by 2 protons. This system was studied at various excitation en- ergies and found to fission consistently in a bimodal pattern, as shown in Figure 4.2. Of 31 Figure 4.2: After projecting the fission fragment mass vs total kinetic energy (TKE) correla- tion (b) onto the TKE axis, the TKE distribution is deconvoluted into high- and low-energy components, centered around the values TKEhigh and TKElow (a). The mass distributions at these particular energies are fitted to a double (c) and single (d) Gaussian, respectively. This procedure is repeated for three different compound nucleus excitation energies E* (e-g). Figure from [4]. the sample measured, roughly 1/3 of cases were found to fission symmetrically, while the other 2/3 fissioned asymmetrically with a light-to-heavy mass ratio of approximately 79/99. Furthermore, it was observed that mass-asymmetric fragments tended to have higher kinetic energies than symmetric fragments. To better interpret the results of this experiment, we performed DFT calculations us- ing the Skyrme functional UNEDF1HFB [58] and the Gogny functional D1S [96]. These calculations involved computing a PES using the collective coordinates Q20 and Q30. The UNEDF1HFB PES is shown in Figure 4.3, and the D1S PES is in Figure 4.4. A calculation with full Langevin dynamics was not performed because the framework was still being de- veloped at the time; however, the static (minimum-energy) pathways shown in the figures correspond to a fragment split AL=AH (cid:25) 80=98, in good agreement with the experimental ratio AL=AH (cid:25) 79=99. Also shown in Figure 4.3 are nucleon localization functions (recall Section 2.1.3) corre- sponding to various marked configurations in the PES. Along the symmetric path (ABcd in 32 Figure 4.3: UNEDF1HFB potential energy surface for 178Pt as a function of multipole mo- ments Q20, representing elongation, and Q30, representing mass asymmetry. The static pathway is marked in purple. Note the two different trajectories ABCD and ABcd and their corresponding NLFs. Figure from [4]. the figure), the fragments appear highly-elongated, including shortly before scission. Since elongation tends to minimize the Coulomb repulsion between fragments, then this configu- ration might be expected to lead to fragments with relatively low kinetic energies. On the other hand, compact fragments such as those in ABCD will tend to have a larger Coulomb repulsion, resulting in compact asymmetric fragments with a higher kinetic energy, which is in qualitative agreement with experiment. Comparing the UNEDF1HFB PES in Figure 4.3 to the D1S PES in Figure 4.4a, one may note that, despite the inherent differences between the functionals, the overall topography of the PES is similar in both cases. In fact, the topography in both cases is quite flat, which differs drastically from the familiar double-humped structure commonly seen in actinides [97] and suggests (in agreement with the observed data) the possibility of a competition between symmetric and asymmetric fission modes. Additionally, as shown in Figure 4.4b, there appears an additional channel corresponding to compact symmetric fragments. This channel 33 Figure 4.4: a) D1S potential energy surface for 178Pt as a function of multipole moments Q20, representing elongation, and Q30, representing mass asymmetry. The static pathway is marked in red. b) With the quadrupole moment fixed to Q20 = 190b, a PES is constructed as a function of Q40, which relates to the size of the neck, and Q30. This PES shows the existence of two symmetric modes, separated by a saddle point with a height of (cid:24)4 MeV. Figure from [4]. is higher in energy than the elongated symmetric pathway and is blocked by a (cid:24)4 MeV saddle (cid:3). point, but it is conceivable that this channel might be more easily accessed at higher E 4.3 The physical origin of fragment asymmetry in the neutron-deficient sub-lead region By now it is now well-established that many isotopes in the neutron-deficient sub-lead region will fission asymmetrically. Microscopic-macroscopic calculations correctly predict asymme- try in the yields, and static DFT calculations do even better by predicting the peak of the distribution to within 1-2 nucleons. In that sense, it may be said that the phenomenon is understood. That said, a simple, intuitive interpretation or explanation of that phenomenon is still a matter of some debate. After the discovery of asymmetric fission in 180Hg, one of the earliest theoretical re- 34 sponses [92] argued that no simple explanation or intuitive interpretation is possible, and that the best explanation for the observation, however unsatisfying, is simply the sum total of all the physics that gets wrapped together to compute a PES. Fragment yields are de- termined by subtle interplays in local regions of the potential energy surface, while magic number-based arguments are just a lot of fortuitous hand-waving that happens to work well for actinides. The fact of asymmetric fragments, in this explanation, is something that hap- pens in spite of, and not because of, fragment shell effects (which, two of the authors later note [98], are rather weak in this region). By contrast, it is argued in [99] that the fragment masses are determined by the shell structure of two well-formed prefragments connected by a thin neck. The nucleons in the neck are then redistributed to either of the prefragments somewhere along the way from saddle point to scission. This argument was further reinforced in a follow-up paper [93], in which it is shown that the asymmetric path is associated with stronger shell effects than the symmetric path. Two more recent papers [100, 35] fundamentally agree with this second perspective, except they argue that one should consider the shell structure of octupole-deformed prefrag- ments instead of familiar spherical magic numbers. By their reasoning, octupole softness would favor fragment distributions with Zlight (cid:25) 34 and Nheavy (cid:25) 56. Most recently, an argument is presented in [94] which likewise invokes deformed shell gaps; however, in this case the argument is made for octupole-deformed (or rather, mass asymmetric) parent nuclei, as opposed to octupole-deformed fragments. Taken together, these suggest that fragment asymmetry is the result of a coupling between highly-deformed shell gaps in the parent nucleus and (possibly-deformed) shell gaps in the fragment nuclei. The separation process might also introduce some additional complications, 35 which may involve dissipation and particle emission. So far, no microscopic dynamical description has been given for fission in this region; it is impossible to follow a fission trajectory all the way to scission, since a combination of very high fission barriers and relatively low-energy “fusion” configurations lead to sharp discontinuities in the PES. It would be interesting to see how the collective inertia behaves for these nuclides, and how they affect the fission dynamics. 36 Chapter 5 Cluster decay in 294Og 5.1 Cluster emission of superheavy nuclei The region of superheavy elements, defined as the set of nuclides with Z > 104, is an interesting one for the study of spontaneous fission because the liquid drop model predicts that all isotopes with Z > 104 are unstable with respect to spontaneous fission. These nuclei receive additional stability due to shell effects, but they nevertheless remain short-lived and many of them will decay by alpha decay and spontaneous fission. Experimentally, spontaneous fission has been observed from several superheavy isotopes (those marked in green in the right panel of Figure 5.1). Other observed superheavy elements undergo a chain of several alpha decays followed by spontaneous fission. Furthermore, a variety of models predict regions where spontaneous fission dominates. One example is shown in the left panel of Figure 5.1, in which branching ratios were estimated by computing and comparing different decay lifetimes obtained from empirical formulas. Figure 5.2 is similar, except that the spontaneous fission lifetimes were computed theoretically, as were the Q(cid:11) values used to estimate alpha decay lifetimes. This still leaves open the question of fission fragments of superheavy elements, however. Several authors have predicted, using models both phenomenological [101, 102, 103, 104, 105, 106, 107] and microscopic [108], that some superheavy elements will undergo a highly- asymmetric form of fission in which the parent nucleus decays into two fragments, with 37 Figure 5.1: Calculated (left) and experimental (right) decay modes for heavy and superheavy elements. The theoretical estimates are based on an analysis of lifetimes calculated via empirical formulae. The boxed isotopes in the left panel are those which have been measured experimentally. Isotopes falling inside the 1(cid:22)s contour are predicted to live longer than 1(cid:22)s. Figure adapted from [5]. Figure 5.2: Dominant decay modes for superheavy elements are indicated based on predic- tions obtained in a nuclear DFT-based framework are indicated. The label “RS” stands for “reflection symmetric”, meaning the least-action wavefunction at the saddle is reflection symmetric, and “NRS” stands for “non-reflection symmetric.” The legend indicates predicted half-lives on a logarithmic scale. Figure from [6]. 38 the heavy fragment near the doubly-magic nucleus 208Pb. This particular decay mode is significant enough to have earned its own name in the literature, where it is known variously as cluster emission, cluster radioactivity, or lead radioactivity [109, 110, 111, 112, 113]. The phenomenon of cluster emission was first observed in 1984 in the decay 223Ra!209Pb + 14C [114] and has since been observed in several actinides. In all cases seen so far, it is a rare event with a small branching ratio [112]. The mechanism of cluster emission is based on the stability of the doubly-magic nucleus 208Pb. 294Og is a particularly excellent candidate for cluster emission because the cluster it is predicted to emit, 86Kr, receives additional stability due to its having a magic number of neutrons N = 50. Semiempirical arguments based on the nuclear symmetry energy lend additional support to this candidate, since 294Og and 208Pb have a similar N=Z ratio [108]. We took this prediction one step further by calculating the full spontaneous fission frag- ment distribution of 294Og using a microscopic approach. As will be shown, the distribution is sharply-peaked around 208Pb, and is quite robust with respect to various inputs. A vi- sualization of the process using the nucleon localization function shows strong evidence of a 208Pb prefragment being formed early in the post-saddle stage of the evolution. 5.2 Predicted spontaneous fission yields of 294Og Potential energy surfaces were computed for each of three EDFs: UNEDF1HFB [58], a Skyrme functional which was optimized to data for spherical and deformed nuclei, including fission isomers; SkM* [56], another Skyrme functional designed for fission barriers and surface energy; and D1S [96], a Gogny parametrization fitted on fission barriers of actinides. A comparison of the different PESs as a function of multipole moments Q20 and Q30 is shown 39 Figure 5.3: Comparison of the PESs for 294Og in the (Q20; Q30) collective plane obtained in UNEDF1HFB (a), D1S (b), and SkM* (c) EDFs. The ground-state energy Egs is normalized to zero. The dotted line in each figure corresponds to E0 (cid:0) Egs = 1 MeV, which was used to determine the inner and outer turning points. The local energy minima at large deformations are marked by stars. Figure from [7]. in Figure 5.3. Despite the quantitative differences between the functionals, one can see that the surfaces strongly resemble one another. Some common features we highlight are: a symmetric saddle point occurring around Q20 (cid:25) 40 b; a second barrier beginning around Q20 (cid:25) 100 (cid:0) 120 b along the symmetric fission path; the presence of local minima at large deformations (marked by stars in the figure); a deep valley that leads to an highly-asymmetric split; and a secondary, less-asymmetric fission valley that emerges at large elongations. But there are differences as well, such as the height of the first saddle point, the depth of the highly-asymmetric fission valley, and the height of the ridge separating the two fission valleys. As a result, the outer turning points are pushed to larger elongations in D1S and SkM* as compared to UNEDF1HFB. These differences in the PES topography strongly affect the predicted spontaneous fission half-lives (cid:28)SF, which in the case of UNEDF1HFB, SkM* (cid:0)2 s, respectively (see also [115, 116] and D1S are 9:1 (cid:2) 10 (cid:0)5 s and 3:2 (cid:2) 10 (cid:0)9 s, 4:0 (cid:2) 10 40 153045UNEDF1HFB153045MassasymmetryQ30(b32)D1S050100150ElongationQ20(b)153045SkM*0369Energy(MeV)(a)(b)(c)294Og for a detailed discussion of half-lives). These large variations of (cid:28)SF reflect the well-known exponential sensitivity of spontaneous fission half-lives to changes in the quantities entering the collective action (2.25). The (cid:28)SF predictions of UNEDF1HFB and, to a lesser degree, SkM* are incompatible with experiment, as 294Og is known to decay by (cid:11)-decay with a half- life of 0.58 ms [117]. The D1S result should not be taken too seriously, either, since they were performed in a smaller collective space leading to overestimation of the half-lives [118, 119]. It is to be noted, however, that while half-lives are very sensitive to details of the calcula- tions, the models used here are very consistent with each other and with experiment when it comes to global observables, such as alpha-decay energies, deformations, and radii [120, 121]. It will be demonstrated below that spontaneous-fission mass and charge yields are also ro- bustly predicted (though experimental confirmation has not yet occurred). This robustness will be shown by varying the EDF, the collective space, the collective inertia, and the strength of the Langevin dissipation tensor (see Section 2.2). The EDF and the collective space were varied together; a different collective space was used for each EDF in the tunneling portion of the PES. The UNEDF1HFB calculation was carried out in a four-dimensional collective space consisting of the collective coordinates (Q20; Q30; Q22; (cid:21)2). By examining this PES, we were able to reduce the dimensionality of the PES for the other two functionals. The SkM* calculation was performed in a piecewise- continuous space. Between the ground state and the isomer state, Q30 does not play a significant role and so the system was described using (Q20; Q22; (cid:21)2). Likewise, Q22 (cid:25) 0 once the isomer state is reached and reflection symmetry is broken (Q30 ̸= 0), suggesting use of the coordinates (Q20; Q30; (cid:21)2). Finally, for the functional D1S we decided to see how our yields would be affected if we did not consider dynamical pairing fluctuations or axial symmetry breaking at all, which results in a drastic reduction to computation time. Those 41 calculations were performed in the collective space (Q20; Q30). In the region of Langevin dynamics, we found again that Q22 (cid:25) 0. Consequently, this region was limited to two dimensions (Q20; Q30) in all cases. The resulting yields in the (N,Z)-plane are shown in Figure 5.4. As expected, the yields are peaked in the region of 208Pb with a sharp fall-off. Likewise, the projected distributions onto the mass and charge axes show a clear preference for cluster emission, regardless of functional used, as seen in the top panels of Figure 5.5. As discussed in Chapter 2, the collective inertia can also have a large impact on the fission dynamics. Using the UNEDF1HFB functional, we compared our result with the non-perturbative cranking ATDHFB inertia, MA, to two other expressions for the inertia which are computationally less expensive to compute and therefore commonly-used in large- scale calculations: the perturbative ATDHFB inertia MAP (which appears smoothed-out compared to MA [73]), and the perturbative GCM inertia MGCM (which is smooth and also lower in magnitude than MA or MAP by roughly a factor of 1.5 [73]). The result, shown in the middle panels of Figure 5.5, show that the distribution has shifted slightly, but that MAP and MGCM yield identical, or nearly-identical results. The smoothness of the perturbative inertias apparently allow fluctuations to drive the system to more extreme fragment configurations. This suggests that the magnitude of the inertia matters less than its topography for computing fission yields (though we note that this would not be true for calculating half-lives, which depend exponentially on the magnitude of the inertia). We also vary the strength of dissipation tensor by adjusting the parameter (cid:17). Our starting point (cid:17)0 = ((cid:17)11; (cid:17)22; (cid:17)12) = (50(cid:22)h; 40(cid:22)h; 5(cid:22)h) is taken from [2], where it was obtained by adjusting (cid:17) to match the experimental fragment distribution of 240Pu. Shown in the bottom panels of Figure 5.5, we find that fluctuations do not affect the peak of the distribution, 42 Figure 5.4: Fission fragment distributions for 294Og obtained in UNEDF1HFB (a), D1S (b), and SkM* (c) EDFs using the non-perturbative cranking ATDHFB inertia and the baseline dissipation tensor (cid:17)0 = ((cid:17)11; (cid:17)22; (cid:17)12) = (50(cid:22)h; 40(cid:22)h; 5(cid:22)h). Known isotopes are marked in gray [8]. Magic numbers 50, 82, and 126 are indicated by dotted lines. Figure from [7]. 43 %yield294Og fission yieldsUNEDF1HFBD1SSkM*(a)(b)(c)0.20.40.60.81.01.21.41.6>1.802550750255075Z050100150N0255075 Figure 5.5: Upper panel: Predicted heavy fragment mass (left) and charge (right) yields of 294Og using different functionals (top, linear scale). Lower panels: collective inertias (middle, logarithmic scale) and dissipation tensor strengths (bottom, logarithmic scale). The baseline calculation was performed using the UNEDF1HFB functional in a 4D space with non-perturbative cranking ATDHFB inertia and dissipation tensor strength (cid:17)0. Figure from [7]. 44 consistent with the results of [122, 123, 69]. The primary effect is in the tails, suggesting that fluctuations affect only a relatively-small number of nucleons. 5.3 Fragment formation in 294Og The Langevin dynamics approach is useful for calculating fragment yield distributions, but it does not offer much physical insight into the process of fragment formation. In order to better understand fragment formation for the most-probable fragments, we computed the nucleon localization function (section 2.1.3 and References [1, 69]) along the cluster emission path. This is shown in Figure 5.6, where it is compared to the spherical fragments 208Pb and 86Kr. We found that the lead prefragment is well-localized just outside the outer turning line. The N (cid:25) 50 neutrons belonging to krypton are also well-localized at early stages of the evolution, but the Z (cid:25) 36 protons are not, highlighting the importance of shell closures in prefragment formation. 5.4 Experimental search for cluster emission in 294Og The effort to detect cluster emission from superheavy elements is underway. However, one of the biggest challenges standing in the way of observation of cluster emission from 294Og is the problem of statistics. There have only been 5 recorded instances of 294Og, from the first observations in 2005 [124] to the most recent in 2018 [117].1 To maximize the possibility of detecting spontaneous fission events as they happen, there is some discussion in the experimental community about building an ionization chamber for superheavy element 1In fact, there were reports as early as 2004 [125] that fission fragments resulting from the spontaneous fission of 294Og may have been detected, but these remain unconfirmed. 45 Figure 5.6: Nucleon localization functions for several deformed configurations of 294Og. For comparison, localizations are also shown for the fragments 208Pb and 86Kr on the left side of each subplot. The configurations shown correspond to Fig. 1 from the paper, with multipole 2 ) ((cid:24)outer turning line); (c) (140 b, moments (Q20; Q30) = (a) (75 b, 0); (b) (120 b, 17 b 24 b 2 ) ((cid:24)scission). Figure from [7]. 3 3 2 ); and (d) (264 b, 60 b 3 46 research with the ability to distinguish fragments by their Z-value. This is discussed at some length in [117]. 47 Chapter 6 Fission yields for r-process nuclei 6.1 The role of fission in the astrophysical r process One of the outstanding mysteries in astrophysics is the origin of heavy elements. It is thought that heavy elements are formed in violent, neutron-rich astrophysical environments, such as core-collapse supernovae or neutron star mergers, as the result of a rapid neutron capture process. This process, generally shortened to “r process”, involves bombardment of small seed nuclei with neutrons followed by beta decay toward the valley of stability. This process is illustrated schematically in Figure 6.1. Fission plays an important role in the r process. The competition between neutron capture and beta decay determines the manner in which the r process proceeds. Fission places an upper limit on the mass of isotopes that can be produced in an r-process scenario, as fission lifetimes for heavy and superheavy isotopes become small enough to compete with neutron capture and beta decay for very heavy systems. Also, as will be discussed shortly, fission fragments of isotopes in this mass region are thought to shape the rare earth peak around A (cid:25) 165. Elemental abundance patterns in some old, metal-poor stars outside our solar system have been found to resemble that of our own [126], and it has been suggested that fission may be responsible for this “universal r process” via a mechanism called fission cycling [127]. Once a heavy nucleus fissions, its fragments can then absorb more neutrons and make their 48 Figure 6.1: Schematic overview of the r process. The competition between neutron capture and beta decay increases the mass of a seed nucleus until it reaches the region where fission dominates (spontaneous, beta-delayed, or neutron-induced). way back up the r-process chain. In a neutron star merger environment, the system may undergo somewhere between 1-4 fission cycles before settling into its final state [128]. In order to guide experimental and theoretical efforts and to reduce uncertainties in the abundance pattern, sensitivity studies (in which model inputs are tweaked to measure their effect on the final yield) can estimate the impact of a particular set of observables, such as fission fragment yields, on the r-process abundance pattern. Sensitivity studies indicate that fission yields primarily affect the elemental abundances between the second r-process peak (A (cid:25) 130) and the rare earth peak (A (cid:25) 160) [129, 9]. This can be seen in Figure 6.2, in which four different sets of phenomenological fission fragment yields were used to compute abundances. That this region is sensitive to fission yields should come as no great surprise, since most naturally fissile nuclei produce fission fragments in this region. R-process network calculations combine nuclear physics inputs, such as decay lifetimes 49 fissionbeta decayneutron capture Figure 6.2: Final abundance patterns from four different simulations of neutron star merger ejecta, each using a different model for the fission fragment yields. The results are pre- sented in two graphs solely for visual clarity. The dots represent the known solar r-process abundance pattern. Figure from [9]. and fission fragment yields, with astrophysical inputs, such as temperature and neutron den- sity, to simulate the complex competition between various nuclear decays. Such calculations require quality inputs from a variety of sources. Many nuclei of interest to the r process are too neutron-rich or too heavy to be studied in a laboratory, which means that theoretical calculations provide essential model inputs. These calculations should be performed using models that are known to extrapolate reliably. As a global microscopic model that is de- signed to reproduce data across the nuclear chart, nuclear DFT is well-suited for this type of problem. 6.2 Fission fragment yields for r-process nuclei The predictive power of r-process network calculations could be improved with a global sur- vey of microscopically-calculated fission fragment yields, including yields from spontaneous fission, beta-delayed fission, and neutron-induced fission. However, the computational ex- pense associated with full Langevin calculations makes this method impractical for large-scale 50 survey calculations. Instead, we have developed an efficient model for estimating fragment yields using the nucleon localization function which is designed for large-scale calculations like this. Our model is introduced in Section 2.2.5, with a detailed description in Appendix A. Initial calculations were performed using the isotopes 254Pu, a neutron-rich isotope of plutonium predicted to have a large fission flow during neutron star mergers [10], and 290Fm, an extremely neutron-rich isotope of fermium. For each isotope, a 1D fission trajectory was calculated up to the outermost isomer state by constraining Q20 and leaving the other degrees of freedom unconstrained. Between the isomer state and the outer turning line, the collective space was expanded to include Q30 in order to account for mass asymmetry between the fragments. The WKB method was used to estimate relative probabilities along the outer turning line, at which stage fragment pairs were identified using the localization- based method described in Appendix A. Sample localizations, 2D PESs in the region of the outer turning line, and the resulting distributions are shown in Figures 6.3 and 6.4. Looking first at the 254Pu localizations, we see a spherical heavy prefragment and a highly-elongated lighter prefragment connected by a thick neck. The spherical prefragment is well-localized, leading to a narrow uncertainty band in the overall yield. However, because of the large number of nucleons remaining in the neck, the statistical redistribution of neck particles produces a broad distribution of potential fragments. It would be interesting to see whether the peak will narrow and the lighter fragment will become better localized well-beyond the outer turning line, but that is beyond the scope of this project. Compared to 254Pu, the 290Fm nucleus is not as well-localized at the outer turning line, and a neck is only just starting to form in the protons. The neutron localization shows a ring- or disk-like of neutrons forming around the neck, which may lead to a large number of prompt neutrons being emitted. Because of this, the error bands associated with 51 (a) 254Pu neutron localization (b) 254Pu proton localization (c) 254Pu PES (d) 254Pu yield Figure 6.3: Localization function, PES, and yield distribution for the spontaneous fission of 254Pu. The stars in the PES indicate outer turning points which were used to compute the cumulative yields. 52 150175Mass0.010.1110% yieldMin energy path onlySum over many OTPs506070Charge0.010.1110 prefragment identification and the yields are much wider than for 254Pu. On the other hand, the overall distribution is much narrower, suggesting that symmetric fission will dominate for this nucleus. Yield contributions from mass-asymmetric outer turning points slightly broaden the neutron distribution, but they have almost no effect on the proton distribution, which remains peaked at Z = 54. 6.3 Future survey-level calculations for the r process Even with methods such as ours, which can reduce the computational burden, a survey- level calculation of fission fragment distributions across the nuclear chart, or even across the region relevant for r-process astrophysics, would be a major undertaking. A more efficient approach would be to isolate a small subset of fissioning nuclei which most affect the final abundance pattern, and make yield predictions for these first. A modified sensitivity study was performed in [10] in which, instead of assessing uncer- tainties, the goal was to isolate “hot spots” of nuclei that would have the greatest impact on the final r-process abundance pattern. Four different nuclear structure models were used to estimate these hot spots, with fission fragment yields provided by the phenomenological GEF model [130]. The neutron-induced fission hot spots are shown in Figure 6.5, with the four different results superimposed on top of one another. It should be noted that Figure 6.5 shows isotopes relevant specifically for their neutron- induced fission yields; however, the GEF yields on which the calculation was based do not show a strong dependence on excitation energy for many nuclei of interest. Therefore, as a first attempt, one might start simply by calculating the spontaneous fission yields for the isotopes indicated in Figure 6.5. A proper treatment of neutron-induced fission using a 53 (a) 290Fm neutron localization (b) 290Fm proton localization (c) 290Fm PES (d) 290Fm yield Figure 6.4: Similar to Figure 6.3, but for 290Fm. 54 145155165175Mass0.010.1110% yieldMin energy path onlySum over many OTPs506070Charge0.010.1110 Figure 6.5: Heatmap showing isotopes whose fission yields are especially relevant to the r process. This combines the results from four different nuclear structure models, and counts how many of the models found each isotope to have an integrated fission flow above a certain threshold. Figure from [10]. finite temperature formalism is under development, but many challenges still remain (see Appendix B as well as [93, 131, 132]). 55 12413414415416417418419420421482879297102107Z (Proton Number)1234(n,f) Chapter 7 Conclusions The overarching goal of the project, in which this dissertation is a part, is to describe sponta- neous fission yields in a self-consistent framework based on microscopic nuclear physics. We highlight spontaneous fission in three distinct and relatively-unexplored regions of the nuclear chart. New phenomena such as the observed asymmetry of fragments in the neutron-deficient sub-lead region, superheavy element synthesis and decay, and r-process nucleosynthesis force us to grapple with our understanding of fission and ultimately lead to greater insight and understanding of the fission process. New computational tools are supporting these devel- opments by allowing us to use more efficient, more detailed models in our calculations. In Chapter 4, we estimated the peak of the fragment distribution corresponding to spon- taneous fission of 178Pt. In addition, we were able to qualitatively characterize properties of the fragments using properties of the PES, in a way that was compatible with experimental data. In Chapter 5, we used Langevin dynamics to estimate the fission yields belonging to the superheavy element 294Og. We predicted that this nucleus may fission according to an exotic form of fission called cluster emission, which is driven by the strong shell closures associated with N = 126 and Z = 82. We also took this calculation as an opportunity to test the robustness of microscopic fission calculations using Langevin dynamics, and we found that the peak was robustly predicted to within a few nucleons even as the inputs changed. In Chapter 6, we began using a new method for estimating fragment yields (see Ap- 56 pendix A) to calculate fission fragment distributions that are relevant to the astrophysical r process. This method is designed to be efficient, which will allow it to be used for survey-level calculations that might later be used as inputs in r-process network calculations. With that said, there is much still to accomplish on the way to a comprehensive micro- scopic description of spontaneous fission. This includes additional physics modeling as well as basic technical developments. The following few pages give an overview of some of these future prospects. 7.1 Improved model fidelity Qualitatively and quantitatively, microscopic models are useful for understanding the com- plex physics of the fission process. To be useful for many applications, however, the level of quantitative agreement with experiment needs to improve significantly. Using the case of 240Pu as a benchmark [2], the Langevin approach can predict spontaneous fission fragment yields to within around 30% accuracy. Depending on the application, this number must be brought down to (cid:24)5% to be useful. Half-lives, which are an essential defining element of spontaneous fission, are especially difficult to estimate. Assuming that the minimum action is computed to the same 30% accuracy as the yields (which is admittedly a rough estimate, since these quantities depend on two different regions of the PES with drastically different dynamics), then, because of the exponential dependence of the half-life on the action, the half-life may be off by several orders of magnitude. Since fission lifetimes are largely determined by the topography of the PES through which the nucleus must tunnel, it is important to improve predictions of both the potential energy 57 and the collective inertia. The perturbative inertia is easy and fast to compute, but not expected to be as reliable as the non-perturbative inertia. And as discussed in Section 3.2, the finite difference formulae from which the non-perturbative inertia is often computed are subject to errors due to level crossings and numerical artifacts. One tool that would eliminate the need for finite-difference derivatives is automatic differentiation. Automatic differentiation essentially keeps track of arithmetic operations as they are performed and combines their derivatives according to the chain rule. Because this is performed as the code runs, it adds almost no additional runtime and no numerical error. Once implemented, automatic differentiation would allow derivatives to be computed with machine precision, without the need to compute any additional nearby points. The chief drawback is that automatic differentiation can be difficult to implement. As for the potential energy, the current state-of-the-art functional used for fission - UN- EDF1 - tends to underestimate the height of saddle points on the PES, which in turn has the effect of artificially reducing half-lives. One might suggest refitting the functional using addi- tional fission data, but unfortunately, the UNEDF2 project [133] concluded essentially that Skyrme-like functionals have been pushed to their limits. A next-generation set of functionals based on the UNEDF project has been proposed which uses an expansion of the density ma- trix and a link to modern effective forces, both of which can be systematically improved [134]. This set of functionals has not been extensively tested, save for some initial benchmarking, so it is not clear whether there will be any added benefit to using them. It may be that new functionals such as these will better reproduce fission data, but it should be noted that, while the calculations presented here would not be affected, introducing an explicit density-dependence can cause problems for beyond-mean field corrections [135, 136, 137]. It was assumed throughout this dissertation that fissioning nuclei behave as superfluids 58 maintained at temperature T = 0. However, a proper description of fission from an excited compound nucleus as in the case of Chapter 4, or neutron-induced fission as in Chapter 6, requires a finite-temperature formalism. This is discussed in Appendix B. A known problem when constructing PESs is the existence of discontinuities, which occur when one reduces a complex system with hundreds of degrees of freedom into only a few collective coordinates [138]. This can obscure important physics, such as saddle point heights and the complex physics of scission. Brute force computation with a larger number of collective coordinates (such as we did in Chapter 5) will become more accessible as computers continue to become more powerful. However, we might see a faster turnaround on investment if we start using collective coordinates that better describe the shape of fissioning nuclei. For instance, a set of coordinates D; (cid:24) was proposed as an alternative to the multipole moments Q20 and Q30 in [71], where D is the distance between prefragment centers of mass and (cid:24) = (AR (cid:0) AL)=A characterizes the fragment asymmetry. Using these coordinates, the authors found that they were better able to describe scission configurations than if they had used multipole moments instead. Finally, to make the calculations described in this dissertation fully self-consistent will require the development of a self-consistent description of dissipation. In addition to con- tributing to a more reliable model, it may be that a better understanding of the exchange mechanisms between intrinsic and collective degrees of freedom could allow one to more reliably model fragment kinetic and excitation energies. The topic of dissipative motion in quantum systems is being studied theoretically [139, 140, 141, 142, 80], but the ideas are not yet sufficiently well-developed for quantitative fission calculations. 59 7.2 Computational tools As we look for and find new ways to study and apply fission, we might consider the availability of new computational tools, such as machine learning and GPUs. As stated earlier, our biggest bottleneck is the calculation of the PES. The expense of calculations has the side effect of forcing us to use small collective spaces with relatively-few collective variables, which may mean a loss of important physics. Therefore, it is worth investing in methods to reduce the burden of PES calculations. One potential solution is a machine learning method called speculative sampling. Speculative sampling works well for problems that have an expensive computational model, and calculations that must be performed for many inputs or in a large parameter space. The essence of this technique is to create an emulator for your model using machine learning, then initiate many calculations using the expensive model and terminate early any calculations which do not appear to diverge strongly from the emulator. In this way, we don’t waste time calculating things that the emulator can already predict, and instead focus our efforts on those portions of the PES which provide new information. After this set of calculations converges, the emulator can be retrained using the new data and the procedure is repeated. This could be implemented for our problem in the following way: Start by using Monte Carlo sampling to begin constructing a PES in a large collective space using a traditional DFT solver, and then use those results to train a PES emulator. Then use Monte Carlo sampling to start a new round of PES calculations, only this time compare the results from the DFT solver in real-time to the emulator’s prediction. As long as the DFT solver diverges from the prediction it should be allowed to continue, but once it appears that the DFT calculation agrees with the prediction the calculation terminates. The emulator should then 60 be retrained using the new data, and the feedback cycle continued until the PES is predicted sufficiently-well by the emulator. Although traditional simulations will likely continue to be performed on traditional CPUs, machine learning performs well on GPUs. That means this technique will be able to take advantage of many of the newest architectures. For microscopic calculations to become practical for r-process calculations, though, would require an efficient way of estimating fission properties on a large scale, ideally covering the whole nuclear chart. To do this using conventional calculations will still be a tremendous undertaking, even with the potential speedup from speculative sampling. Another machine learning paradigm called “transfer learning” might be useful here. In transfer learning, one starts with a set of data which is large, but imperfect. For fission, this might mean yields computed globally via a phenomenological model such as GEF [130]. This data is used to train a model, such as an artificial neural network, so that essentially one has an emulator that understands the physics of the cheap model. This knowledge is then transferred (hence the name) to a new model with new data, which in our case might be experimental data or even the results from detailed microscopic calculations. In a neural network with many layers, this transfer might take place by retraining the model while holding all but the last two or three of layers constant. In our example, the first several layers contain most of the basic physics of fission, and then the last few layers contain corrections and improvements to the simple model. 7.3 Outlook It is an exciting time to study fission. New phenomena are being discovered as fission exper- iments are performed in unfamiliar regions of the nuclear landscape, and there is a healthy 61 feedback loop between theory and experiment thanks in great part to modern computational tools. It is, in many ways, the culmination of years of theoretical development, but there is still much to be done. 62 APPENDICES 63 Appendix A Prescription for fragment identification via localization function As mentioned in Section 2.2.5, we have developed a method for estimating primary fission yields using the nucleon localization function and a PES which has been computed out to the outer turning line. In fact, a 1D trajectory might be sufficient for some purposes, such as to estimate the location of the peak of the distribution. Our prescription is the following: • Identify possible prefragments using the density and nucleon localization function (see Figure A.1): – Along the primary axis of the nucleus, find the two locations (one for each frag- ment) at which the density is widest. This is taken to be the approximate center of each fragment. – For each fragment, start from the central value and locate the two nearest optima in the localization (maximum or minimum; one on either side). This gives you a range of possible central values (see Figure A.1), which leads to a built-in uncertainty band in the calculated yields. – Sum up the nucleons above (below) the central value and multiply by two to determine the identity of the upper (lower) prefragment. 64 √jEneck=aj, where • The remaining neck nucleons are assumed to form into clusters: alpha particles first, then a deuteron if there is one of each nucleon left, plus some remaining free protons or neutrons. • The total binding energy of the prefragments and any particles/clusters is subtracted from the total energy of the system to get the energy of the neck, Eneck. – Use this energy to define the temperature of the system: T = the level density parameter a is taken to be a = Ap=10 MeV(cid:0)1 and with Ap the mass number of the parent nucleus. • Use brute force combinatorics to find every possible way of assigning the neck particles to either of the prefragments (including the possibility that a particle/cluster might be emitted instead of flowing to a prefragment). • For each fragments + particles combination, look up the binding energy of each of the final fragments in a table (such as MassExplorer [143]). • The Coulomb repulsion energy between the fragments plus the total binding energy of the fragments and any remaining clusters is subtracted from the total energy of the system to get the residual (kinetic + intrinsic) energy, Eres. • This energy defines a state in the grand canonical ensemble (note that the chemical potentials are hidden in the binding energies): Z = ∑ (cid:0) Eres kB T . e (cid:0) Eres(n1;z1;n2;z2) • Calculate probabilities: P (n1; z1; n2; z2) = 1 Z e kB T . • Fold those probabilities with a Gaussian smoothing function. 65 (a) 240Pu neutron localization function (b) 240Pu proton localization function Figure A.1: Localization function for 240Pu with lines drawn to represent prefragment iden- tification As a test case to benchmark this method, we estimate the 240Pu spontaneous fission yield, which was also estimated theoretically using Langevin dynamics in [2]. Experimental data exists for this nucleus as well, in [11, 12]. These results are shown in Figure A.2. Comparing the localization-based results with those from experiment and from Langevin dynamics, the peak of the yield appears to be identified in all cases to within just a few nucleons of one another. That the tails are poorly-reproduced or that the uncertainty bands cover roughly an order of magnitude should not be surprising, as it was shown in [2, 69] that the tails of the distribution are due mainly to a dynamic interplay between fluctuations and the collective inertia which takes place beyond the outer turning line. As an additional test, we also show results computed for the nucleus 294Og in Figure A.3. Like the Langevin results, the yield is sharply peaked near 208Pb. The discrepancy of several nucleons in the peak when comparing between Langevin and the localization-based method proposed here is comparable to the discrepancies found in Chapter 5 when comparing 66 (a) 240Pu PES (b) 240Pu yield Figure A.2: PES and yield distribution for spontaneous fission of 240Pu. On the left, the stars in the PES indicate outer turning points which were used to compute the cumulative yields. On the right, the black curve in the yields corresponds to the Langevin result, and experimental data comes from [11, 12]. 67 125150175Mass0.010.1110% yieldNew method - Min action path onlyNew method - Sum over many OTPsLangevinExperiment506070Charge0.010.1110 between functionals, collective inertias, and values of the dissipation tensor strength. As in the 240Pu case, the tails of the distribution are poorly-reproduced compared to the Langevin result, especially on the more-symmetric side. The fact that yield peaks can be estimated at the outer turning line instead of at scission leads to a tremendous reduction in the number of PES points which must be computed. Furthermore, our experience with the WKB approximation shows that most trajectories leading to the outer turning line actually follow the minimum action path, branching away only just before reaching the outer turning line. Since the yields computed in this manner depend not on the absolute WKB probabilities, but on the relative probability of points along the outer turning line, we propose further reducing the number of PES calculations by computing a 1-dimensional fission trajectory between the ground state and the outer turning line (either the static path or the minimum action path), and a 2-dimensional PES in the region of the outer turning line. 68 (a) 294Og PES (b) 294Og yield Figure A.3: Similar to Figure A.2 but for 294Og. 69 150200Mass0.010.1110% yieldNew method - Min action path onlyNew method - Sum over many OTPsLangevin75100Charge0.010.1110 Appendix B Finite-temperature ATDHFB collective inertia The calculations presented in this dissertation assumed that the system was maintained at temperature T = 0 and that the nucleus behaved as a superfluid. However, in future calculations, in order to properly describe fission from an excited compound nucleus as in the case of Chapter 4, or neutron-induced fission as in Chapter 6, one should technically use a finite-temperature formalism. In the T ̸= 0 case, pairs may be broken and various excited states are accessible, resulting in changes to the PES and the collective inertia. A finite-temperature formalism for calculat- ing the PES has already been developed (see [132, 93]), as have perturbative expressions for the finite-temperature adiabatic time-dependent Hartree-Fock-Bogoliubov (FT-ATDHFB) inertia with cranking [144, 145, 146]. Since the yields in Chapter 5 showed a dependence on whether the collective inertia was computed perturbatively or non-perturbatively, this chapter is dedicated to deriving a non-perturbative finite-temperature expression for the collective inertia. We will begin with a brief overview of ATDHFB theory and finite temperature HFB. Then we present the resulting collective inertia without detailed derivation, which is straightfor- ward. 70 B.1 Adiabatic time-dependent HFB Adiabatic time-dependent Hartree-Fock, which is generalized to adiabatic time-dependent Hartree-Fock-Bogoliubov, is presented carefully in [147, 148]; we merely summarize the rele- vant results here. The fundamental assumption of Time-Dependent Hartree-Fock-Bogoliubov (TDHFB) theory is that a system is described by a Slater determinant. This assumption allows us to write the TDHFB equation: i(cid:22)h _R = [H;R] ; where ~H = 0B@ h (cid:0) (cid:21) ∆ (cid:3) (cid:0)h (cid:3) (cid:0)∆ + (cid:21) 1CA ; 0B@ (cid:26) (cid:3) (cid:0)(cid:20) ~R = 1CA : (cid:20) 1 (cid:0) (cid:26) (cid:3) (B.1) (B.2) One additional assumption, which has already been discussed in Chapter 2, is that col- lective motion is slow compared to single particle motion of the system. This is called the adiabatic approximation, and the consequent model is called Adiabatic Time-Dependent Hartree-Fock-Bogoliubov (ATDHFB). Both here and in Chapter 2, the adiabatic approxi- mation is invoked in order to describe the dynamics of a fissioning nucleus in terms of a small set of collective variables and their derivatives (analogous to coordinates and veloci- ties): suppose there is a set of collective variables (q1; q2; : : : ; qn) which describes changes in the density R(t) at all times t. Then n∑ (cid:22)=1 _R = @R @q(cid:22) : _q(cid:22) (B.3) Once the system is described in terms of collective coordinates and velocities, the energy 71 can be expressed as the sum of a “potential” term (which depends on the coordinates) and a “kinetic” term (which depends on the velocities). Our goal is to understand the kinetic part of the energy, which in some sense describes the dynamics of a deformed nucleus. A key component of this will be the inertia tensor M, which plays the role of the “mass”: Ekin (cid:24) 1 M _q2. 2 B.2 Finite-temperature density As in any statistical theory, one first must determine which sort of ensemble properly de- scribes the system. Nuclei are finite systems which have a conserved number of particles; however, HFB theory explicitly breaks particle number symmetry. In principle, we should use a microcanonical ensemble to describe a nucleus as a closed, isolated system, but that turns out to be challenging to solve because it requires a full knowledge of the eigenspec- trum of the nucleus. Based on the particle non-conserving nature of HFB theory, we instead describe our system using the grand canonical ensemble. The formalism to do this is shown in [149]. In the end, finite-temperature HFB looks almost identical to standard zero-temperature HFB, except the density in the quasiparticle (B.4) 0B@ 0 0 1CA ! 0B@ f 1CA ; 0 0 1 (cid:0) f basis is replaced by R = 0 1 with the Fermi factor f given by f(cid:22) = 1 1+e(cid:12)E(cid:22) . 72 B.3 ATDHFB equations With the adiabatic assumption in place, we can write the density as an expansion around some time-even zeroth-order density: R(t) = ei(cid:31)(t)R0(t)e (cid:0)i(cid:31)(t) = R0 + R1 + R2 + : : : ; where (cid:31) is assumed to be “small” (which is explained more rigorously in [148]) and R1 = i [(cid:31);R0] ; R2 = [[(cid:31);R0] ; (cid:31)] : 1 2 The HFB matrix, being a function of R, is likewise expanded: H = H0 + H1 + H2 + : : : ; (B.5) (B.6) (B.7) (B.8) (B.9) and R and H are then substituted into the TDHFB equation (B.1). Gathering terms in powers of (cid:31) yields i(cid:22)h _R0 = [H0;R1] + [H1;R0] ; i(cid:22)h _R1 = [H0;R0] + [H0;R2] + [H1;R1] + [H2;R0] : (B.10) (B.11) These two equations are the ATDHFB equations. They can be solved self-consistently to find both (cid:31) and R0; however, this is rarely done in practice. A more common approach is 73 to exploit the fact that solutions to the ATDHFB equations are (by design) close to true HFB solutions. One can then take an HFB solution and compute its time derivative by the first ATDHFB equation (B.10) to get ATDHFB-like behavior instead of a true solution to the ATDHFB equations. Working in the HFB quasiparticle basis, in which 0B@ E 0 0 (cid:0)E 1CA ; H0 = 0B@ f 0 0 1 (cid:0) f 1CA ; (B.12) R0 = one finds by combining equations (B.7) and (B.10) that the ATDHFB equations can be written as: ab + (fb (cid:0) fa)H11 (1);ab; (cid:22)h _R11 (cid:22)h _R12 (cid:22)h _R21 (cid:22)h _R22 (0);ab = (Ea (cid:0) Eb)(fb (cid:0) fa)(cid:31)11 (0);ab = (Ea + Eb) (1 (cid:0) (fa + fb)) (cid:31)12 (0);ab = (Ea + Eb) (1 (cid:0) (fa + fb)) (cid:31)21 (0);ab = (Ea (cid:0) Eb)(fb (cid:0) fa)(cid:31)22 ab ab + (1 (cid:0) (fa + fb))H12 (cid:0) (1 (cid:0) (fa + fb))H21 ab (1);ab; (1);ab; (B.13) (cid:0) (fb (cid:0) fa)H22 (1);ab; where the superscripts refer to N 2 (cid:2) N 2 subblocks of the full N (cid:2) N matrix. It is com- mon (the so-called “cranking approximation”) to assume that changes in the density have approximately no effect on the mean field, in which case these relations reduce to (cid:22)h _R11 (cid:22)h _R12 (cid:22)h _R21 (cid:22)h _R22 (0);ab = (Ea (cid:0) Eb)(fb (cid:0) fa)(cid:31)11 ab; (0);ab = (Ea + Eb) (1 (cid:0) (fa + fb)) (cid:31)12 ab; (0);ab = (Ea + Eb) (1 (cid:0) (fa + fb)) (cid:31)21 ab; (0);ab = (Ea (cid:0) Eb)(fb (cid:0) fa)(cid:31)22 ab: 74 (B.14) Finally, by expanding the total HFB energy up to second-order in (cid:31), the total energy of the system is found to be E(R) = EHF B + Tr (H0R1) + 1 2 Tr (H0R2) + 1 2 Tr (H1R1) : 1 4 (B.15) The “kinetic energy” of the system is given by the latter two terms, which are both second order in (cid:31): Ekin(R) = Tr (H0R2) + 1 2 Tr (H1R1) : 1 4 (B.16) B.4 Finite-temperature ATDHFB inertia If we solve the ATDHFB equations (B.10) and (B.11) and then substitute equations (B.14) and (B.3) into the resulting expression for the kinetic energy (B.16), we end up with the following result: 2 where the collective inertia M is given by ∑ (cid:22)(cid:23) Ekin (cid:25) 1 _q(cid:22) _q(cid:23)M(cid:22)(cid:23); 0@@R11 0@ @R21 (0);ab @q(cid:22) + @R11 (0);ba @q(cid:23) @R12 (0);ba @q(cid:23) (B.17) (B.18) 1A + 1A35 : @R22 (0);ab @q(cid:22) @R12 (0);ab @q(cid:22) @R22 (0);ba @q(cid:22) @R21 (0);ba @q(cid:22) + 24 M(cid:22)(cid:23) = (cid:22)h2 2 1 (Ea (cid:0) Eb)(fb (cid:0) fa) 1 (Ea + Eb)(1 (cid:0) fa (cid:0) fb) (0);ab @q(cid:22) In the zero temperature case, the first piece (involving R11 and R22) goes to zero. 75 B.5 Numerical implementation Numerical implementations of the finite-temperature ATDHFB inertia may encounter some challenges. First, as discussed in Section 3.2 we know there is an ideal range of values (cid:14)q to use when computing derivatives of R using finite differences. Too large, and the derivative is artificially-smoothed out; too small, and densities of adjacent points can become indistinguishable (dependent on the HFB convergence parameter). Additionally, there may be problems when dealing with small, non-zero temperatures (or, more generally, when fb(cid:0)fa is small). We might consider introducing a cutoff on the difference fb (cid:0) fa. If we make our cutoff too tight, we start cutting off physics (the tail of the Fermi distribution). If we leave our cutoff too small, we start dividing by numbers that are smaller than the noise in the density matrix which will lead to the inertia blowing up. As an example, suppose we take a finite-difference spacing (cid:14)q = 10 (cid:0)3 b and set our HFB (cid:0)7. Roughly-speaking, then, we would expect that we can set convergence parameter to 10 a cutoff for fb (cid:0) fa of something around 10 (cid:0)4 or so and still obtain reasonable results. This was done for a 240Pu calculation (Q20 = 77 b, Q22 = 1 b, chosen at random from the 240Pu PES) with a few different cutoff values. First the exact T=0 inertia was computed. Then the inertia was computed again using the same T=0 densities, but with a fake temperature T=0.05 MeV. The results, shown in Table B.1, indicate the presence of an ideal cutoff parameter that is in the neighborhood of our prediction. 76 Convergence: (cid:14)q: Actual T=0 inertia: Cutoff (cid:0)4 10 (cid:0)5 10 (cid:0)6 10 (cid:0)7 10 (cid:0)8 10 ) (cid:0)7 (cid:0)3 10 10 ( (cid:0)2 1.585632(cid:2)10 (cid:22)h2 Inertia 1.585205(cid:2)10 1.585622(cid:2)10 1.585697(cid:2)10 1.587771(cid:2)10 1.612200(cid:2)10 MeV(cid:1)b2 (cid:0)2 (cid:0)2 (cid:0)2 (cid:0)2 (cid:0)2 MeV b (cid:22)h2 MeV(cid:1)b2 % error (cid:0)2 -2.692933(cid:2)10 (cid:0)4 -6.306634(cid:2)10 (cid:0)3 4.099312(cid:2)10 (cid:0)1 1.348989(cid:2)10 1.675546 Table B.1: As explained in the text, a cutoff for fb (cid:0) fa may be needed to reduce numerical noise in the collective inertia. This table shows how the calculated inertia changes depending on the size of the cutoff. 77 Appendix C List of my contributions • Observation of the competing fission modes in 178Pt – Phys. Lett. B 790, 583-588 (2019) – Performed PES and localization calculations – Created graphics for Figure 3 – Suggested several revisions to the text of the article • Colloquium: Superheavy elements: Oganesson and beyond – Rev. Mod. Phys. 91, 011001 (2019) – Generated Figure 12 – Helped write Section VI: Fission – Suggested comments and revisions to the text • Cluster radioactivity of 294Og – Phys. Rev. C 99, 041304 (2019) – Lead author on paper/wrote first draft – Did all UNEDF1HFB calculations and SkM* Langevin calculations – Generated all figures 78 • Talks – “Cluster formation and emission in 294Og”, Nuclear Fission and Structure of Exotic Nuclei, March 2019, Japan Atomic Energy Agency (JAEA), Tokai, Japan – “Fission in Exotic Nuclei Using Density Functional Theory”, FIRE Collaboration Meeting, May 2018, North Carolina State University – “Cluster Emission in Oganesson-294”, Spontaneous and Induced Fission of Very Heavy and Super Heavy Nuclei Workshop, Apr 2018, ECT*, Trento, Italy – “A ‘Microscopic’ Description of Nuclear Fission”, Contemporary Physics Seminar, Sep 2016, Indiana University South Bend • Fission Tools – Suite of codes which extend the functionality of HFODD, HFBTHO, and other DFT solvers to the problem of nuclear fission. – https://gitlab.com/zachmath/fission_tools – Maintainer, 65 commits – Converted old codes from Fortran77 to Fortran90 – Implemented shared memory (OpenMP) and distributed memory (MPI) paral- lelism – Improved documentation, flexibility, and user-friendliness – Created several Python-based utilities for plotting and file handling – Increased functionality, such as by increasing from 2D to 3D or 4D • PES Tools 79 – Python framework for handling potential energy surfaces as XML files – https://gitlab.com/schuncknf/pes_tools – Developer, 25 commits – Added a point class with various methods to use on a single point of a PES – Interface to fission_tools – Various updates bugfixes • DFTNESS – DFTNESS (Density Functional Theory for Nuclei at Extreme ScaleS) is a com- putational framework to solve the equation of nuclear density functional theory. It is based on the two solvers HFBTHO and HFODD. – https://gitlab.com/schuncknf/dftness – Developer, 60 commits – OpenMP parallelization of subroutine QMULCM – Various updates and bugfixes – NOTE: Most of these commits were related to fission_tools, before that was all moved to a separate repository 80 BIBLIOGRAPHY 81 BIBLIOGRAPHY [1] Chun Li Zhang, Bastian Schuetrumpf, and Witold Nazarewicz. Nucleon localization and fragment formation in nuclear fission. Phys. Rev. C, 94(6):064323, December 2016. [2] Jhilam Sadhukhan, Witold Nazarewicz, and Nicolas Schunck. Microscopic modeling of mass and charge distributions in the spontaneous fission of 240pu. Phys. Rev. C, 93(1):011304, January 2016. [3] A. N. Andreyev, K. Nishio, and K-H Schmidt. Nuclear fission: a review of experimental advances and phenomenology. Rep. Prog. Phys., 81(1):016301, January 2018. [4] I. Tsekhanovich, A.N. Andreyev, K. Nishio, D. Denis-Petit, K. Hirose, H. Makii, Z. Matheson, K. Morimoto, K. Morita, W. Nazarewicz, R. Orlandi, J. Sadhukhan, T. Tanaka, M. Vermeulen, and M. Warda. Observation of the competing fission modes in 178pt. Phys. Lett. B, 790:583–588, March 2019. [5] A V Karpov, V I Zagrebaev, Y Martinez Palenzuela, and Walter Greiner. Super- heavy nuclei: Decay and stability. In Excit. Interdiscip. Phys., pages 69–79. Springer International Publishing, Heidelberg, 2013. [6] M. Warda and J. L. Egido. Fission half-lives of superheavy nuclei in a microscopic approach. Phys. Rev. C, 86(1):014322, July 2012. [7] Zachary Matheson, Samuel A. Giuliani, Witold Nazarewicz, Jhilam Sadhukhan, and Nicolas Schunck. Cluster radioactivity of og294. Phys. Rev. C, 99(4):041304, April 2019. [8] NNDC Interactive Chart of Nuclides. nndc.bnl.gov/nudat2/, 2018. Interactive chart of nuclides, nudat 2.7. [9] M Eichler, A Arcones, A Kelic, O Korobkin, K Langanke, T Marketin, G Martinez- Pinedo, I Panov, T Rauscher, S Rosswog, C Winteler, N T Zinner, and F.-K. Thiele- mann. The role of fission in neutron star mergers and its impact on the r-process peaks. Astrophys. J., 808(1):30, 2015. [10] Nicole Vassh, Ramona Vogt, Rebecca Surman, J. Randrup, Trevor M. Sprouse, Matthew R. Mumpower, P. Jaffke, D. Shaw, Erika M. Holmbeck, Yong-Lin Zhu, and 82 G C McLaughlin. Using excitation-energy dependent fission yields to identify key fis- sioning nuclei in r -process nucleosynthesis. J. Phys. G Nucl. Part. Phys., 46(6):065202, June 2019. [11] J.B. Laidler and F. Brown. Mass distribution in the spontaneous fission of 240pu. J. Inorg. Nucl. Chem., 24(12):1485–1492, December 1962. [12] H. Thierens, A. De Clercq, E. Jacobs, D. De Frenne, P. D’hondt, P. De Gelder, and A. J. Deruytter. Kinetic energy and fragment mass distributions for pu240 (s.f.), 239pu(nth,f), and 240pu((cid:13),f). Phys. Rev. C, 23(5):2104–2113, May 1981. [13] O Hahn and F Strassmann. Uber den nachweis und das verhalten der bei der bestrahlung des urans mittels neutronen entstehenden erdalkalimetalle. Naturwis- senschaften, 27(1):11–15, January 1939. [14] Lise Meitner and O. R. Frisch. Disintegration of uranium by neutrons: a new type of nuclear reaction. Nature, 143(3615):239–240, February 1939. [15] N. Bohr. Disintegration of heavy nuclei. Nature, 143(3617):330–330, February 1939. [16] G.N. Flerov and K.A. Petrjak. Spontaneous fission of uranium. Phys. Rev., 58(1):89– 89, July 1940. [17] C. F. v. Weizsäcker. Zur theorie der kernmassen. Zeitschrift für Phys., 96(7-8):431–458, July 1935. [18] Niels Bohr and John Archibald Wheeler. The mechanism of nuclear fission. Phys. Rev., 56(5):426–450, September 1939. [19] S M Polikanov, V A Druin, V A Karnaukhov, V L Mikheev, A A Pleve, N K Skobelev, V G Subbotin, G M Terakopyan, and V A Fomichev. Spontaneous fission with an anomalously short period. Sov. Physics11 JETP-USSR, 15(6):1016–1021, 1962. [20] David Lawrence Hill and John Archibald Wheeler. Nuclear constitution and the in- terpretation of fission phenomena. Phys. Rev., 89(5):1102–1145, March 1953. [21] S G Nilsson. Binding states of individual nucleons in strongly deformed nuclei. Mat. Fys. Medd. Dan. Vid. Selsk., 16, 1955. [22] V.M. Strutinsky. Shell effects in nuclear masses and deformation energies. Nucl. Phys. A, 95(2):420–442, April 1967. 83 [23] V.M. Strutinsky. “shells” in deformed nuclei. Nucl. Phys. A, 122(1):1–33, December 1968. [24] M. Brack, Jens Damgaard, A. S. Jensen, H. C. Pauli, V. M. Strutinsky, and C. Y. Wong. Funny hills: The shell-correction approach to nuclear shell effects and its applications to the fission process. Rev. Mod. Phys., 44(2):320–405, April 1972. [25] Hans J Krappe and Krzysztof Pomorski. Theory of Nuclear Fission, volume 838 of Lecture Notes in Physics. Springer Berlin Heidelberg, Berlin, Heidelberg, 2012. [26] Peter Möller and Jnorgen Randrup. Calculated fission-fragment yield systematics in the region 74