NOVEL COMPUTATIONAL APPROACHES FOR NUCLEAR FISSION THEORY By Daniel Lay A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Physics — Doctor of Philosophy Computational Mathematics, Science and Engineering — Dual Major 2025 ABSTRACT Nuclear fission is important for energy production, medicinal applications, nonprolifer- ation efforts, and nucleosynthesis studies. The rapid neutron capture process (r process) contributes to observed abundances of medium- and heavy-mass nuclei, and requires fission data for hundreds of nuclei, most of which lie outside of experimental reach. Therefore, predictive fission models with quantified uncertainties are required. In this thesis, spontaneous fission is described as tunneling through an effective barrier defined using a set collective coordinates, called the potential energy surface (PES), which is computed using nuclear density functional theory (DFT). The half-life is then determined by the tunneling pathway, and the primary fragment yields are approximately determined by its endpoint. Computing uncertainties for these quantities in a Bayesian framework requires tens- to hundreds-of-thousands of calculations, making it computationally infeasible. This problem is exacerbated by the large number of nuclei that participate in the r process. This thesis is divided in two parts. The first discusses the formalism necessary for com- puting spontaneous fission observables. Improvements to the tunneling pathway calculation are presented, and the improved methodology is applied to nuclei with competing fission modes. The second part discusses two strategies for approximating these observables across the r process region of the chart, with quantified uncertainties. The first strategy uses neu- ral networks to emulate the PES. The second uses dimensionality reduction techniques to propagate statistical uncertainties from the energy density functional posterior parameter distribution. ACKNOWLEDGMENTS I would first like to thank my advisor, Witek Nazarewicz. Thank you for your guidance, support, and tolerance, and all of the wonderful opportunities you have provided. Working with you has played an enormous role in who I am today, both as a researcher and as a person. Next, I would like to thank the many collaborators I have had the opportunity to publish with: Sylvester Agbemava, Eric Flynn, Pablo Giuliani, Samuel Giuliani, Kyle Godbey, L´eo Neufcourt, and Jhilam Sadhukhan. It has been a pleasure to collaborate on many interesting projects, and I am excited to continue working together in the future. I would also like to thank the other members of my research group, past and present: Maxwell Cao, Xingze Mao, Simin Wang, Tong Li, Menghzi Chen, Josh Wylie, Ante Ravli´c, Sudhanva Lalit, Josh Nicholson, Josh Belieu, Richard Gumbel, Yukiya Saito, Bailey Knight, Yuanzhuo Ma, Mookyong Son, Aaron Philip, An Le, Lauren Jin, Landon Buskirk, and Andrew Yeomans-Stephenson. Simply being in the same research group has been an amazing experience, and I look forward to seeing what everyone does next. I have had the good fortune to be surrounded by quite a number of other graduate students. A number of them joined me at the start of my graduate school journey, and I’ve met many others along the way. They include Grayson Perez, Hieu Le, Mo Hassan, Jordan Purcell, Sheng Lee, Bella Molina, Jason Gombas, Andy Smith, Andrew Sanchez, Peter Farris, Cavan Maher, Julia Hinds, Hannah Berg, Patrick Cook, Danny Jamoa, Nick Cariello, and Yani Udiani. Without you all, my graduate school could not possibly have been as memorable (dare I say, fun) as it was. I would like to express gratitude to my thesis committee members: Heiko Hergert, Ryan LaRose, Hendrik Schatz, and Yang Yang. iii I would like to thank all of the staff who have helped me in my time here, especially Elizabeth Deliyski, Kim Crosslan, and Gillian Olson. I would like to specifically thank Yukiya Saito, Ante Ravli´c, and Sudhanva Lalit for reviewing my thesis and providing excellent feedback. And of course, I want to thank my family. First, I would like to thank my extended family: all of my Aunts and Uncles, all of my cousins and cousins-in-law, and my Grandparents. Thank you for providing a constant backdrop of love and support, and for indulging my clumsy attempts at explaining my work. Next, I would like to thank my brother, Peter, and my sister, KC. Thank you for being there whenever I look up from my work, and keeping me connected to the outside world. Finally, I would like to thank my mom, Jeanne, and my dad, Stephen. Your love and support has made me who I am today. iv TABLE OF CONTENTS Chapter 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Nuclear Fission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Large Scale Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Organization of this Dissertation . . . . . . . . . . . . . . . . . . . . . . . . 2.1 The HFB Equations 2.2 The Energy Density Functional Chapter 2. Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 The Skyrme EDF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 The Gogny EDF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Constrained Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Particle Number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 The Nuclear Shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Dynamical Pairing . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 3. Spontaneous Fission . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 The Tunneling Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 The Collective Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 WKB Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3 An Illustrative Example . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Determining the Least Action Path . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Grid-Based Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 The Nudged Elastic Band . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Initial Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Primary Fission Fragment Yields Chapter 4. Multimodal Fission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 Phenomenological Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.2 Time-Dependent Approaches 4.1.3 A Hybrid Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 The Fermium Chain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Two Collective Coordinates 4.2.2 Three Collective Coordinates Chapter 5. Neural Network Emulation . . . . . . . . . . . . . . . . . . . . . . . 5.1 The Rapid Neutron Capture Process . . . . . . . . . . . . . . . . . . . . . . 5.2 The Neural Network Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Neural Network Quality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.1 The Emulated Quantities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.2 Relation to Observables 5.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v 1 1 2 3 4 4 7 11 13 14 14 17 19 21 21 22 26 28 29 30 32 34 37 38 38 38 40 41 44 44 48 50 52 52 53 55 56 60 64 Chapter 6. Reduced Order Model Emulation . . . . . . . . . . . . . . . . . . . 6.1 Reduced Order Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Illustrative Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 6.2.1 The Schr¨odinger Equation . . . . . . . . . . . . . . . . . . . . . . . . 6.2.2 A Modified Gross-Pitaevskii Equation . . . . . . . . . . . . . . . . . 6.3 The Axial HFB Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.1 Subtleties with the HFB Equations . . . . . . . . . . . . . . . . . . . 6.3.2 The Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.3 The Payoff . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 65 66 66 70 74 74 76 79 81 Chapter 7. Conclusions and Outlook . . . . . . . . . . . . . . . . . . . . . . . . 83 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 vi Chapter 1. Introduction 1.1 Nuclear Fission Nuclear fission, discovered by Otto Frisch, Otto Hahn, Lise Meitner, and Fritz Strassmann in 1938 [1, 2], is the splitting of a heavy nucleus into multiple smaller nuclei. Today, it is important for energy production, medicinal applications, and defense initiatives, especially nonproliferation efforts. It is also a useful testing ground for approaches to the nuclear quantum many-body problem, and plays an important role in describing the rapid neutron capture process (the r process) abundances [3]. The fission process takes place on multiple different timescales, each with a set of asso- ciated observable quantities, see Fig. 2 of Ref. [4]. The first stage is the entrance channel, which is different for induced, β-delayed, and spontaneous fission. The entrance channel for induced fission is a particle (e.g. a neutron or a photon) colliding with a heavy nucleus; for β-delayed fission, it is a heavy nucleus β decaying; and for spontaneous fission, it is a nucleus sitting in its ground state. Next, the nucleus rapidly deforms from a compact shape into a peanut-like shape, and splits into two (rarely, three) fragments. This motion takes place between 10− 21 − 19 10− seconds. The resulting distribution is called the primary fragment yields, and may be recon- structed from experimental data. The kinetic energy of the fragments is also an important observable. The emitted fragments are highly excited. Their excitation energy is carried off by emission of neutrons and gamma radiation, taking place between 10− 18 10− 7 seconds. The − neutron and gamma ray spectrum can both be measured, along with correlations between emitted particles. Finally, the fragments (now in their ground state) are unstable and β 1 decay towards stability. This last process takes place on timescales from microseconds up to hours, and leads to cumulative fission yields. All told, a unified description of fission would require a model that is descriptive from 18 10− − 103 s. Thus, each stage is treated using a different model, typically with inputs from the previous stage. This work focuses solely on the collective motion leading to the primary fragment yields. 1.2 Large Scale Calculations Fission properties of hundreds of nuclei are important for r process nucleosynthesis [5]. A large fraction of the relevant nuclei are short-lived, and in fact outside of experimental reach, even with state-of-the-art facilities such as the Facility for Rare Isotope Beams [6, 7]. It is desirable to model fission based on (effective) nucleon-nucleon interactions that are able to describe a wide swath of nuclear observables, rather than phenomenological approaches tailored to fission and nuclei near stability. Hence, we describe fission using density functional theory (DFT). The tradeoff is that DFT is computationally more expensive than phenomenological ap- proaches. Thus, important observables such as the primary fragment yield distributions tend to be taken from phenomenological calculations [8, 9]. To make matters worse, meaningful predictions require theoretical error bars, requiring tens- to hundreds-of-thousands of DFT calculations for each calculation previously required [10]. Note that this challenge persists for all applications, including comparison with experimental data [11]. Thus, model emulators are necessary. Model emulators may be understood as approx- imate, yet fairly accurate, models that are quicker to evaluate than the true model. Ma- chine learning tools have been successfully used as model emulators in all areas of nuclear 2 physics [12]. Recent developments in nuclear theory, including in DFT, have focused on reducing the number of model calculations required to train emulators. However, continued development is necessary for realistic calculations of fission properties. 1.3 Organization of this Dissertation This dissertation is organized as follows. DFT is reviewed in Ch. 2. Chapter 3 connects DFT to spontaneous fission, and Ch. 4 applies this approach to nuclei with competing fission modes. These two chapters follow my publications, Refs. [13, 14], respectively. Chap- ter 5 discusses fission emulation using neural networks, following my publication, Ref. [15]. Ch. 6 discusses DFT emulation using reduced order models, which is ongoing work. Finally, conclusions and prospects are presented in Ch. 7. 3 Chapter 2. Density Functional Theory Our approach to fission is based on nuclear density functional theory. DFT was introduced in computational chemistry, see e.g. Refs. [16, 17]. In nuclear DFT, one expresses the total energy of the nucleus as a functional of local densities and currents: the energy density functional (EDF), then uses the variational principle to minimize the total energy. Its use is justified by the Hohenberg-Kohn theorems, which prove that, for every many-body interaction, there exists an EDF whose lowest-energy configuration has the same energy as the ground state of the many-body Hamiltonian [18, 19]. DFT has had great success describing many phenomena, including ground-state properties such as masses and charge radii, excited states and decay rates, collisions, and fission, see Ref. [20] and references within. In this chapter, I will provide a brief overview of DFT, focused on the description relevant for fission. First, I will give a condensed derivation of the Hartree-Fock-Boguliubov (HFB) equations in terms of a generic EDF. Second, I will describe the EDFs used throughout this work. And third, I will connect to fission by discussing constrained calculations. 2.1 The HFB Equations The Hartree-Fock-Boguliubov equations have been studied in many works in nuclear physics, including but not limited to Refs. [21–28]. Here, I will summarize the features relevant for this thesis. Consider a second-quantized interaction, composed of single-particle creation and anni- hilation operators a†i , aj: (cid:88) ˆH = tija†i aj + ij 1 4 (cid:88) ijkl ¯vijkla†i a†jalak + · · · (2.1) The one-body term tij is typically the kinetic energy, and ¯v is the antisymmetrized two-body 4 term. The creation and annihilation operators obey the canonical Fermionic anticommuta- tion rules ai, a†j} { = δij and ai, aj} { = 0. Rather than computing the full spectrum of the Hamiltonian, we use the Ritz variational principle [22]. For any trial wave function is not orthogonal to the true ground state wave function, Eϕ ≡ ˆH ϕ ϕ | | ⟩ ⟨ ϕ ϕ ⟩ | ⟨ Egs, ≥ that ϕ ⟩ | (2.2) where Egs is the true ground-state energy of the Hamiltonian. The independent-particle picture is a zero-order depiction of nuclear motion, hence the success of shell-model-type approaches [29, 30]. As such, an appropriate ansatz for a wave function is the Slater determinant, which describes nucleons moving independently within a mean field. Following this line of thought leads to the Hartree-Fock (HF) equations. However, pairing superfluidity is important for many nuclear properties, such as the well-known odd- even mass staggering in nuclear binding energies [31]. The simplicity of the independent particle picture can be retained, even while largely capturing the pairing effects, by instead discussing the independent quasi -particle picture. The following derivation largely follows from Ref. [23]. The HFB equations may also be derived using the Thouless representation of the HFB wave function; see the discussion in [22]. In mathematical terms, the HFB wave function HFB ⟩ | is a Slater determinant expressed in terms of quasiparticle operators α†i , αj via the Boguliubov transformation: αk = (cid:88) n (U ∗nkan + V ∗nka†n), α†k = (cid:88) (Vnkan + Unka†n). n (2.3) 5 Note that these operators also obey the anticommutation rules αi, α†j} { = δij, αi, αj} { = 0 (equivalently, the Boguliubov transformation is unitary), so U, V are not independent. The HFB wave function HFB | ⟩ is a quasiparticle vacuum, i.e. it satisfies αj| HFB ⟩ = 0 for all j. The variational energy can be computed straightforwardly using Wick’s theorem: EHFB = ⟨ HFB ˆH HFB | | ⟩ HFB HFB | ⟨ ⟩ (cid:20)(cid:18) = tr t + (cid:19) (cid:21) Γ ρ 1 2 1 2 − tr[∆κ∗]. The particle, ρ, and pairing, κ, one-body densities are defined as ρij = (V ∗V T )ij, κij = (V ∗U T )ij, and the mean-field potential and pairing field are Γik ≡ δE δρki (cid:88) = jl ¯vijklρlj, ∆ij ≡ δE δκ∗ji = 1 2 (cid:88) kl ¯vijklκkl. (2.4) (2.5) (2.6) As the name suggests, the pairing field ∆ describes pairing phenomena, and is not present in HF theory. The HFB equations are then derived by varying EHFB with respect to ρ and κ indepen- dently, with a constraint to preserve unitarity of the Boguliubov transformation. This yields the self-consistent eigenvalue problem       U V H µ = Eµ       U V µ    , = H h ∆ ∆∗ − h∗ −    . (2.7) The matrix H is called the HFB matrix, and it is a functional of one-body densities ρ, κ. Eµ 6 is the quasiparticle energy of the µ-th eigenstate. Since time-reversal symmetry is imposed, each positive eigenvalue Eµ has a corresponding state with negative eigenvalue, Eµ [24]. − In practice, then, only the positive-energy states must be solved for, and the corresponding negative-energy states are fully occupied in the HFB vacuum. Typically, one considers one set of equations for protons and neutrons separately, as is done in this work (thus, isoscalar pairing that mixes proton and neutron channels [32] is ignored). There are two schematic approaches to solving the HFB equations that are relevant for this thesis. The first iteratively diagonalizes the HFB matrix, and the eigenvectors at the n-th iteration are used to construct the HFB matrix at the (n+1)-th iteration. This is the method used in the axial solver HFBTHO [33] and the 3D solver HFODD [34], for example. The second approach is to directly minimize the HFB energy using a gradient-based approach, as has been used in the axial solver HFBAxial [35]. As a terminology aside, DFT may be understood as HF(B) theory where the EDF is postulated instead of the many-body interaction. Properly, one should denote whether DFT is carried out at the HF or the HFB level; throughout, we consider DFT at the HFB level. 2.2 The Energy Density Functional In position space, the total energy EHFB can be written as an integral over the energy density, EHFB[ρ, κ] = (cid:90) d3r E [ρ(r), κ(r)], where E is the EDF. Quite generically, the EDF can be written as (r) = E Ekin + ECoul + Eph + pp, E 7 (2.8) (2.9) where Ekin and ECoul are the kinetic and Coulomb energy densities, respectively. correspond to particle-hole (ph) and particle-particle (pp) interaction terms. Eph and Eph depends pp depends on both ρ and κ. Each term will be E only on ρ and its derivatives, while pp E explored further below. Note that this separation may lead to divergences in beyond-mean- field approaches in the event that Eph and pp are not derived from the same effective E interaction, see the discussion in Ref. [27]. This thesis does not go beyond the mean field level, and so I will ignore this subtlety. These terms are composed of local densities and currents. Such densities can be cat- egorized based on their behavior under time reversal symmetry, although is even under E time-reversal symmetry. For the ground state of an even-even nucleus, because of time- reversal symmetry, the time-odd densities vanish [28], so we will only introduce the time-even densities. We begin with the one-body density matrix, written in position space as ρq(rσ, r′σ′) = (cid:88) k Vkq(rσ)V ∗kq(r′σ′). (2.10) The index q labels the nucleon species (q = n, p) and σ, σ′ are the spin indices. Unlike in HF theory, single-particle states have partial occupancies, i.e. the norm of the lower components of the eigenvectors, Nkq = (cid:82) d3r (cid:80) Vkq(rσ) σ | 2, can take on any value between | 0 and 1. Therefore, the sum over k is over all states (at least, all those allowed by the pairing regularization procedure - see e.g. Ref. [36] for a discussion with a zero-range pairing interaction). The particle density is further decomposed as ρq(rσ, r′σ′) = 1 2 ρq(r, r′)δσσ′ + 1 2 (cid:88) µ (σ σµ | σ′)ρµ,q(r, r′), | (2.11) 8 where the σi are the usual Pauli matrices. The local particle density ρq(r), kinetic energy density τq(r), and spin-current density Jµν,q(r) are then defined as ρq(r) = ρq(r, r), τq(r) = Jµν,q(r) = r ∇ 1 2i ∇r′ ρq(r, r′) |r′=r, ( ∇ µ − ∇′µ)ρν,q(r, r′) |r′=r. (2.12a) (2.12b) (2.12c) Upon varying the EDF, one finds a single-particle potential from terms involving ρq, an effec- tive mass from terms involving τq, and a spin-orbit potential from terms involving Jµν,q [28, 37, 38]. These densities have analogues in the particle-particle channel [24], written in terms of the pairing density matrix ˜ρq(rσ, r′σ′) = (cid:88) − k Vkq(rσ)U ∗kq(r′σ′). (2.13) It is related to κ Eq. (2.5) by a flip of one of the spin indices, ˜ρq(rσ, r′σ′) = 2σ′κq(rσ, r′σ′). − In place of ∆ we use ˜hq(rσ, r′σ′) = − is Hermitian and time-even [24, 36]. 2σ′∆q(rσ, r′σ′). This choice is convenient because ˜ρq The kinetic energy density is simply Ekin = (cid:88) ℏ2 2mq q τq(r), (2.14) with mq the nucleon mass. Commonly, one sets ℏ2/(2mq) = 20.73 . . . MeV fm2, to varying precision depending on the EDF fit (see e.g. Table II of Ref. [28]). Subtracting off the center-of-mass correction changes 1/m 1/m(1 − → 1/A) [22, 39], although some calibrations 9 such as UNEDF1 [40] and UNEDF2 [41] do not include this correction. The Coulomb energy density contains the direct and exchange terms, ECoul ≡ E D Coul + E Coul. E The direct term depends on ρp(r), and can be evaluated directly: D Coul = E e2 2 (cid:90) d3r′ ρp(r) ρp(r′) r′| r | − , (2.15) (2.16) where e is the electron charge. Numerical techniques for evaluating this integral have been discussed in Refs. [25, 36, 38]. The exchange term, on the other hand, is [42, 43] E Coul(r) = E e2 2 − (cid:88) kk′σσ′ Vkp(rσ)V ∗k′p(rσ) (cid:90) d3r′ V ∗kp(r′σ′)Vk′p(r′σ′) − r′| r | , (2.17) which depends on ρp(rσ, r′σ′). While some works evaluate this exactly, see e.g. the Fortran code HFBTHO v4, [44], it is somewhat numerically expensive. Instead, one commonly uses the Slater approximation to write E Coul(r) E ≈ − (cid:19)1/3 e2 3 4 (cid:18) 3 π 4/3 p ρ (r), (2.18) and extensions depending on only ρp(r) have been discussed in [45]. Typically, the total error in the exchange term is of order 0.2-1 MeV [46, 47], even for superheavy nuclei such as 306126 [45]. Throughout this work, we use the Slater approximation for the exchange Coulomb term. Nuclear physics is contained solely in the particle-hole and particle-particle terms, Eph 10 and E pp, respectively. Over the last few decades, many EDFs have been developed. The most common nonrelativistic EDFs are the Skyrme [22, 48] and Gogny [49, 50] functionals, both of which have many parameterizations that have been developed. Also of note are the Fayans functionals, which improve on descriptions of surface-based nuclear properties [51, 52] and the BCPM [53] and SeaLL1 [54] functionals, which aim to minimize the total number of adjustable parameters in their respective EDFs. Efforts are also ongoing to inform EDFs based on ab initio theory, see e.g. [55–59]. Finally, covariant EDFs have been developed using the relativistic formulation of DFT [60, 61]. While many of the techniques developed throughout this thesis apply quite generally, we will focus on the Skyrme- and Gogny-type EDFs. They are discussed in the following. 2.2.1 The Skyrme EDF The Skyrme EDF is motivated by the density matrix expansion. For the short ranged nucleon-nucleon interaction, the full one-body density matrix ρ(r, r′) can be expanded about r r′ ≈ − 0, see e.g. Ref. [28]. In the particle-hole channel, the (time-even) nuclear part of the EDF is commonly written in the form [20, 41, 62] Eph,Sk(r) = (cid:88) t=0,1 χt(r), (2.19) which is the sum over the isoscalar (t = 0) and isovector (t = 1) components. These are, in turn, written as χt(r) = Cρρ t ρ2 t + Cρτ t ρtτt + CJJ t (cid:88) µν Jµν,tJµν,t + Cρ∆ρ t ρt∆ρt + Cρ t ∇ J ρt Jt. ∇ · (2.20) 11 Jt is shorthand for (cid:80)3 Here, ∇ · coefficients Cuu′ t µνη=1 εµνη∂ηJµν,t, with ε being the Levi-Civita tensor. The are real constants, except for the traditional density-dependence in Cρρ t , Cρρ t ≡ t0 + Cρρ Cρρ t1 ργ. This term descends from a three-body contact force that is equivalent to a density-dependent two-body force in the case γ = 1, see Ref. [37]. Following Ref. [63], γ is now taken as an adjustable parameter in the theory. The original Skyrme EDF was derived and used with HF theory [22, 37, 38, 48], which does not include pairing effects. Pairing has often been accounted for using simple models in the Bardeen-Cooper-Schrieffer (BCS) approximation, such as in the SkM∗ parameteri- zation [64]. On the opposite extreme, the same density matrix expansion can be carried out in the pairing channel. The end result is that one simply replaces ph-densities with pp-densities everywhere in Eq. (2.20) except in the density-dependent coefficient Cρρ t . Many modern EDFs, such as the UNEDF series of functionals [39–41], use a simplified local density dependent pairing EDF: pp = E (cid:88) q (cid:20) 1 Vq 2 − (cid:18) ρ0(r) ρc 1 2 (cid:19)β(cid:21) ˜ρ2 q(r), (2.21) see Ref. [28] and references therein. The isoscalar particle density is ρ0(r) = ρp(r) + ρn(r). The coefficients Vq, switching density ρc, and exponent β, are additional fit parameters. Note that this contact pairing leads to divergences in the density matrix as the maximum quasiparticle energy is increased; this has been addressed in e.g. [24, 25, 36] simply by cutting off the maximum quasiparticle energy at (typically) 60 MeV. By renormalizing the pairing strength, the dependence of the HFB energy can be made independent of this cutoff strength [65]. Over the years, many different parameterizations of the Skyrme EDF have been devel- 12 oped. For a non-exhaustive list, see the references in Refs. [28, 66]. A common challenge with many parameterizations is that they are optimized to ground-state properties, and are there- fore unsuitable for fission calculations. For instance, both the SkM [67] and UNEDF0 [39] EDFs predicted actinide fission barriers far below their accepted values (note that fission bar- riers are not, strictly, observable quantities; nevertheless, they have commonly been treated as such in the literature). To remedy this, the SkM∗ EDF adjusted the surface parameters from the original SkM EDF to better reproduce fission barrier heights [64]. Similarly, the UNEDF1 EDF was calibrated including excitation energies of fission isomers in actinide nu- clei [40]. Because of these improved fission properties, the SkM∗ and UNEDF1 EDFs are the only two Skyrme-type EDFs that will be discussed in this work. 2.2.2 The Gogny EDF The Gogny EDF is a finite-range EDF. The finite range nature was desired initially so that the same interaction could be used in the ph and the pp channel, without the pairing divergences of the Skyrme EDF. Instead, the pairing is regularized by the finite-range. For further discussion, see Ref. [27]. The ph part of the Gogny EDF is typically written in the literature as an effective two- body interaction, as ˆV (r1, r2) = (cid:88) (r1− e− r2)2/µ2 i (cid:0)Wi + Bi ˆPσ i=1,2 Hi ˆPτ − − Mi ˆPσ ˆPτ (cid:1) + t0 (cid:0)1 + x0 ˆPσ (cid:1)ρα (cid:18) r1 + r2 2 (cid:19) δ(r1 − r2) + iWLS(σ1 + σ2) (←− ∇ 1 − ∇ 2) ←− · × δ(r1 − r2)(−→ ∇ 1 − ∇ 2), −→ (2.22a) (2.22b) (2.22c) see Refs. [28, 33]. The operators ˆPσ, ˆPτ are the spin and isospin exchange operators, respec- 13 tively. The gradients operate to the left or right corresponding to the direction of the arrow above them. The terms µi, Wi Bi Hi Mi t0, x0, α, and WLS are fit parameters. This two-body interaction consists of terms that are similar to the Skyrme effective interaction, which can be found in Ref. [37]. Indeed, expanding the first term, Eq. (2.22a), r2 = 0 precisely gives the contact terms in the Skyrme EDF. The approximate about r1 − three-body force, Eq. (2.22b), matches the density-dependent coefficient Cρρ t [ρ] in the same approximation, and Eq. (2.22c) corresponds to the spin-orbit term. An analogous interaction (switching ρ → potential [28]. ˜ρ and so on) is commonly used in the pp channel, as is the Coulomb As with the Skyrme EDF, a number of Gogny-EDF parameterizations exist. For a recent review, see Ref. [27]. Also as with the Skyrme EDF, early parameterizations such as the original D1 parameterization [68, 69] were fit excluding fission information, leading to a too-high fission barrier in actinide nuclei. As a result, the D1S parameterization was refit to match the 240Pu barrier height [70]. The D1S parameters are given in Ref. [71], although the sign of the spin-orbit coupling WLS is incorrect. Again, as this thesis is primarily concerned with fission, the D1S parameterization is the only Gogny EDF that will be considered. 2.3 Constrained Calculations The HFB equations can be solved with various constraints imposed. Typical constraints relevant for fission include the total number of nucleons, shape deformations such as multipole moments and the neck width [21], and dynamical pairing fluctuations [72]. Below I describe these constraints, why they are used, and how they are implemented in this work. 2.3.1 Particle Number The particle number (PN) operator for nucleonic species q is, in second-quantized form, ˆNq = (cid:80) k a†kqakq. Notice that the HFB wave function is not an eigenstate of this operator, 14 i.e. the HFB wave function does not have definite PN. This is problematic for describing finite nuclei, which have a definite number of nucleons. Exact PN restoration can be carried out, in which the HFB wave function is projected onto a desired eigenstate of ˆN , see the discussion in [22, 23]. Certainly, this has been carried out in a number of works, see e.g. Refs. [44, 73–76]. As discussed in Ref. [76], while exact PN restoration is important for some observ- ables, ground-state properties of even-even nuclei (binding energies, deformations, etc.) are fairly insensitive to this symmetry breaking. Further, exact PN restoration (especially the variation-after-projection approach [22]) is computationally expensive as compared to sin- gle HFB calculations. Third, there exist problems when applying projection techniques to EDFs [27, 77]. As such, no EDFs have been calibrated with restored PN. Since the EDF is calibrated to some experimental data, the choice of many-body method should be considered as part of the calibration [78], meaning that (for instance) an EDF calibrated without PN restoration should not be used with PN restoration. Therefore, in this thesis, we do not consider exact PN restoration. A first approach to approximate PN symmetry is to constrain the average particle num- ber, ˆNq ⟨ , to the desired value. Here and throughout, I use the notation ⟩ ˆ O⟩ ≡ ⟨ ⟨ ˆ HFB O| | HFB ⟩ (2.23) for a generic operator ˆ O . This is done for protons and neutrons separately. This has the effect of replacing the total energy with EHFB[ρ, κ] EHFB[ρ, κ] → (cid:88) − q=n,p λq( ˆNq ⟨ ⟩ − Nq), (2.24) 15 with λq the chemical potential (alternatively called the Fermi energy [76]). Following the discussion in Ref. [24], we consider only cases with negative chemical potential λq < 0, as nuclei are particle-unbound for λq > 0. States with quasiparticle energy Eq,µ | | < λq are − discrete, and the upper component Vq,µ is localized. These in turn correspond to a localized density matrix, which corresponds to a localized nucleus. The continuous spectrum corre- sponding to Eq,µ | | > λq is discretized, either in a basis expansion or by a coordinate mesh. − The chemical potential is updated to adjust ˆNq ⟨ ⟩ to the desired value at each iteration [36]. Note that in the Skyrme case, a correct evaluation of the particle number requires the pairing regularization discussed above. A next approach to approximate PN symmetry is the Lipkin-Nogami (LN) method [36, 76, 79–82]. In this approach, a particle-number dispersion term is added to the total energy: EHFB[ρ, κ] → EHFB[ρ, κ] (cid:88) − q=n,p λ2q( ˆN 2 ⟨ q ⟩ − ⟨ ˆNq 2). ⟩ (2.25) The parameter λ2q is also adjusted iteration-by-iteration, with expressions given in Refs. [36, 81] tending to assume an effective seniority pairing for simplicity. A comparison between LN and variation-after-PN-projection has been carried out in Ref. [76] for the Ca and Sn isotopic chains, where the pairing energy difference between the two methods is observed to be fairly small (between 0.5-1.5 MeV). I note that projection after LN has also been presented in that work, with an overall improvement of the order of hundreds of keV, but again no EDFs have been calibrated with the projected-LN approach. In this sense, LN may be concluded as an adequate approximation of exact PN restoration, especially for ground-state calculations. It is also of note that this is true even when the HFB equations predict a collapse of pairing correlations, see Ref. [28] and references within. 16 The UNEDF1 functional mentioned above has been calibrated at both the HFB [83] and the LN [40] level. The fit at the LN level has marginally smaller errors on the observables used for fitting, although as the authors point out, this is no guarantee that this EDF is strictly better [83]. In this work, I will distinguish the EDF used, denoting by UNEDF1HFB (UNEDF1LN) the calibration done at the HFB (LN) level. 2.3.2 The Nuclear Shape Fission is an example of large-amplitude collective motion, and the motion is assumed to be driven by a small number of collective coordinates. These coordinates typically parameterize the shape of the nucleus in some way, thus the HFB equations are solved when the nucleus is constrained to said shape. It is also worth noting that HFB wave functions corresponding to different shapes are commonly used in the generator coordinate method approach, and indeed this approach can be used to describe much of the fission phenomena that will be described in later chapters, see Ref. [22]. In this thesis, we parameterize the nuclear shape using the multipole moments [84] ˆQλµ(r, θ, ϕ) = rλY ∗λµ(θ, ϕ), (2.26) with Yλµ(θ, ϕ) the usual spherical harmonics. There are a number of other parameterizations used in the literature. Typically these are associated with microscopic-macroscopic models of the nucleus, which explicitly parameterize the total energy of the system in terms of a bulk energy term (the “macroscopic” term), and corrections from shell effects (the “micro- scopic” term). Such parameterizations are often chosen specifically for describing elongated, asymmetric deformations common to the fission process. For a review of such parameter- izations, see Refs. [21, 85]. Also common is the neck width, defined in position space as 17 ˆQN (r) = exp[ (r − − rneck)2/a2], where rneck is the position of the neck and a is an arbitrary length scale. This is useful especially in elongated configurations near scission; see again Ref. [21]. We restrict our discussion to axially symmetric configurations, i.e. ˆQλµ⟩ ⟨ = 0 for µ = 0. Note that the triaxial quadrupole moment ˆQ22 has been shown to reduce the fission barrier between the ground state configuration and the fission isomer (FI) [86–88]. This moment is excluded purely for computational reasons, and as a consequence we do not directly compare to experimental half-life measurements. We also restrict the multipolarity of constrained moments to λ have an intuitive connection to fission. The monopole moment ˆQ00⟩ ⟨ particle number, which is already constrained. The dipole moment 4, as these moments ≤ corresponds to the total ˆQ10⟩ ∝ ⟨ z ⟨ controls the ⟩ center-of-mass of the nucleus; it determines whether the nucleus is centered in the coordinate grid, and should in principle vanish automatically. In practice, sometimes this coordinate must be constrained explicitly to Q10 = 0 at least for initial iterations of the solver; typically, ˆQ20⟩ ⟨ controls the elongation of the nucleus. The nuclear ground state is fairly compact, with once a solution is reached, it is stable without this constraint. The quadrupole moment increasing during the fission process. The octupole moment ˆQ20⟩ ⟨ asymmetry of the nucleus. Most nuclei are predicted to have reflection-symmetric ground controls the mass ˆQ30⟩ ⟨ states, so that ˆQ30⟩ ⟨ = 0, although there are exceptions [89]. During the fission process, actinide nuclei develop nonzero ˆQ30⟩ ⟨ , leading to the well-known asymmetric primary fission fragment yields, see the discussion in Ref. [90]. Finally, the hexadecapole moment ˆQ40 is responsible for the number of nucleons in the neck region during the deformation process, see Refs. [14, 91]. These constraints are implemented via the augmented Lagrangian method, which tends 18 ̸ to be more numerically stable than e.g. a quadratic constraint [92]. The total energy is modified to EHFB[ρ, κ] → EHFB[ρ, κ] (cid:88) − a ˆQa0⟩ − λa( ⟨ Qa0), (2.27) where Qa0 is the desired multipole moment. The Lagrange multiplier λa is updated using linear response theory, as if the deviation from the desired multipole moment is due to an external field [50]. For numerical convenience, the cranking approximation of the quasipar- ticle random phase approximation (QRPA) matrix is typically used [93] (see also Sec. 3.1.1). This form is particularly convenient for multiple constraints. Some implementations also include the Lagrange multipliers in the Broyden scheme implemented in Refs. [94, 95], which updates the multipliers iteration-by-iteration as the HFB equations are solved. 2.3.3 Dynamical Pairing As a consequence of broken U(1) particle number symmetry, the HFB ground state may be viewed as deformed in the corresponding gauge space [96]. This deformation is described by the distance ∆ (known as the pairing gap [22]) and the angle φ in the gauge space. One can consider fluctuations around the pairing gap with the minimal HFB energy as an additional collective coordinate to be constrained. Pathways exploring large pairing gaps are expected to lower the spontaneous fission half-life by decreasing the collective inertia tensor, see the discussion in Ch. 3. This has been explored in microscopic-macroscopic studies, see e.g. Ref. [97] and other references within [72]. As discussed in Refs. [98, 99], ∆ is closely related to the particle number dispersion term (∆ ˆNq)2 = ˆN 2 q − ⟨ ˆNq 2, see Eq. (2.25). In those works, the constraint is ⟩ handled similarly to the shape constraints defined above. Alternative approaches to tuning 19 ∆ are discussed in Ref. [21]. A constraint on this quantity has been used in HFB studies of fission in Refs. [13, 72]. There, rather than constraining the particle number dispersion to a particular value, the constraint λ2q was varied to explore the parameter space. Specifically, it is split into the isoscalar and isovector components λ2n ± pairing. Increasing λ2q increases the dynamical pairing fluctuations. As was demonstrated λ2p. Setting λ2q = 0 yields the regular HFB-level in those works, the collective action is lower along the fission trajectories with increased dynamical fluctuations. The latter work [13] will be discussed further in the next chapter. 20 Chapter 3. Spontaneous Fission It seems reasonable to describe fission in a time-dependent framework. The Runge-Gross theorem [100] justifies the use of time-dependent DFT (TD-DFT), much as the Hohenberg- Kohn theorems justify static DFT. Spontaneous fission, however, is a tunneling process, as evidenced by nuclei being metastable to fission. And, as discussed in Refs. [20, 101], TD- DFT is unable to describe tunneling processes. One option, then, is to extend TD-DFT theory using a path-integral approach, although this approach has yet to be used in realistic calculations; see the discussion in Ref. [21]. Another approach is to reduce the TD-DFT equations to a collective Hamiltonian and compute an effective fission barrier, for which tunneling behavior can be studied. That is the approach followed in this thesis. In this chapter, I will first describe our approach to the tunneling problem, using the Wentzel-Kramers-Brillouin (WKB) approach. I will next discuss two numerical algorithms used to determine the tunneling path(s), following Ref. [13]. Finally, I will briefly compare the two approaches in both toy models and realistic cases. 3.1 The Tunneling Process To describe spontaneous fission (SF), we must be able to describe tunneling within our theory. As discussed in [20, 101], one Slater determinant (HF or HFB) describes a single nuclear shape, and thus the time evolution tracks a wave packet rather than the true nuclear wave function. In this sense, the theory functions as a classical approximation to the quantum dynamics, and should not be able to describe tunneling behavior. Tunneling behavior is described only when this classical theory is re-quantized. A direct re-quantization of the TD-DFT equations is computationally prohibitive, pre- venting a re-quantization of the adiabatic TD-DFT equations. With the view of fission as 21 an extreme example of large-amplitude collective motion, we may instead consider collective degrees of freedom, in which the theory may be re-quantized in a computationally approach- able manner. In this section, I will first derive the collective Hamiltonian. I will then describe tunneling in the WKB approximation, and conclude with an illustrative example. 3.1.1 The Collective Hamiltonian This derivation is described in great detail for the HF case in Refs. [22, 101]. It is gener- alized to the HFB case in Refs. [21, 23, 102]. While in many cases the derivation focuses on two-body Hamiltonians, the extension to density-dependent forces is straightforward, by expanding the density dependence to quadratic order in the vicinity of the HF (or HFB) solu- tion. Note that a similar derivation is given using the framework of the generator coordinate method [22]. The time-dependent HFB equation is iℏ ˙ R = [ H , ], where R    = R    ρ κ κ∗ 1 − ρ∗ − (3.1) is the generalized density matrix1. Now, , R H are time-dependent, and the dot is the time derivative. We assume a perturbation of the density matrix [21] (t) = e− iχ(t) R R0(t)eiχ(t), where χ =    χ11 χ12 χ21 χ22    (3.2) χ = 0 (recall the representation [ is a perturbation that is assumed small. R0 is the solution of the static HFB equations with , R , which may be split based on isospin index. In Sec. 2.1, this is why the HFB equations appear independently for protons and neutrons. 1Note that I do not distinguish between protons and neutrons here. This distinction appears in R0] = 0). The block form of χ is similar to that of H R , 22 so χ11 (χ12) is a perturbation in the ph (pp) channel. This form may be interpreted as the interacting-picture representation of the evolu- tion [103], thus the wave function HFB(t) ⟩ | R0. Expan- sion of the HFB energy to second order in χ is carried out conveniently using the Thouless the solution to the HFB equations corresponding to generalized density may be written as HFB(t) ⟩ | ψI (t) ⟩ | ψI (t) ⟩ | with = eiχ(t) representation of the HFB wave function, see Chapter 9 of Ref. [23]. This leads to the form [21] E[ ] R ≈ E[ R0] + (cid:18) ℏ 2 χ12, † χ12,T (cid:19)       . ˙ 12 0 R ˙ R 12 ∗0 (3.3) Note that only the pp perturbation contributes to the evolution at this order, underscoring the importance of pairing fluctuations. Finally, χ is related to ˙ R via the linear response matrix , as M    ˙ 12 0 R ˙ R 12 ∗0    = M    ,    χ12 χ12 ∗ leading to the quadratic form of the total energy: E[ ] R ≈ V [ R0] + ℏ2 (cid:18) 2 12, †0 ˙ R 12,T ˙ 0 R (cid:19) 1 M−    .    ˙ 12 0 R ˙ R 12 ∗0 (3.4) (3.5) 23 The linear response (equivalently, QRPA) matrix    A B B∗ A∗    = M (3.6) is defined as2 Ars,mn ≡ δrmδsn(En + Em) + ∂2E ∂κ∗rs∂κmn , Brs,mn ∂2E ∂κ∗rs∂κ∗mn , ≡ (3.7) where the derivatives include only terms in the EDF defined in terms of κ. This is written as a block matrix, with indices rs α, mn → → β. Similarly, ˙ R 12 0,rs → 12 0,α [104]. ˙ R Note that this derivation assumes a small perturbation χ. The appropriate nearby HFB configuration R0 is not obvious a priori. Solving the unconstrained HFB equations results in a (local) minimum of the energy landscape, regardless of the initial guess, which may not be similar to a tunneling configuration. Take, for instance, a configuration at the outer turning point as shown in Fig. 3.1 - it is totally unlike the ground-state configuration, and thus χ would be large. The appropriate configuration, therefore, is the solution to the constrained HFB equations. Note also that the total energy is now that of a classical system, with coordinates 12 0 . R Thus, one could in principle re-quantize from here and describe the tunneling process, in effect replacing the following discussion with an extremely high-dimensional analogue. This is computationally prohibitive, and largely unnecessary. Instead, one assumes that a set of 2Note that some authors represent δrmδsn(En + Em). M using the Thouless representation, and so do not separate out 24 collective coordinates qα { } drives the time evolution of the system. Via the chain rule, ˙ R0 = (cid:88) α ∂ R0 ∂qα ˙qα, (3.8) the total energy is E = V (q) + 1 2 (cid:88) αβ ˙qαMαβ(q) ˙qβ, Mαβ ≡ ∂ ℏ2 (cid:88) ab,cd R0,ab ∂qα M− 1 ab,cd ∂ R0,cd ∂qβ . (3.9) The potential energy V (q) is computed using constrained HFB calculations, and is called the potential energy surface (PES). An illustrative example is presented in Sec. 3.1.3. The tensor Mαβ is called the collective inertia tensor. The collective inertia tensor is challenging to compute exactly, as doing so requires invert- ing the linear response matrix M . Recent work has demonstrated a reduced computational cost using the finite amplitude method, see Refs. [105, 106], but as yet this has not been applied to fission calculations. The so-called cranking approximation includes only the quasi- particle energies Eα, writing Mab,cd = (Ea + Eb)δacδbd [102]. The linear response matrix is then diagonal, and its inversion is trivial. For this reason, we exclusively use the cranking approximation. In the quasiparticle basis, we use the shorthand from Refs. [15, 102] to write the collective inertia tensor as Mαβ(q) = ℏ2 (cid:88) ij ∗ij F β F α ij + F α Ei + Ej ijF β ∗ij . (3.10) The matrix elements F α ij are computed from derivatives of R0 with respect to the collective 25 coordinates, as [104] F α = U † ∂ρ ∂qα V ∗ + U † ∂κ ∂qα U ∗ − V † ∂ρ∗ ∂qα U ∗ − V † ∂κ∗ ∂qµ V ∗. (3.11) The derivatives are traditionally computed using a three-point Lagrange formula [107–109], as constrained calculations do not converge to the precisely-desired constraint value. An alternative to evaluating these derivatives is given by the perturbative cranking approxima- tion [102, 109, 110], in which the derivatives are approximated by matrix elements of the deformation operator. The perturbative approximation leads to large differences in the fis- sion pathway, as well as the tunneling half-life, as compared to the nonperturbative case [109]. Nevertheless, we use both the perturbative and nonperturbative cranking inertia. So far, the collective coordinates have not been specified. In the classically-allowed region, E > V (q), the dynamics can be described by either the collective Hamiltonian or the TD- HFB equations. Hence, one can determine qα { } from the direct time evolution [20, 111]. Alternatively, one may postulate qα { } using some physical intuition. In simple cases, such as reactions between two α particles, however, there is significant disagreement between the two approaches, and the postulated coordinate is shown to lead to an incorrect inertia [20]. In the classically forbidden region, E < V (q), since the TD-HFB equations do not hold, the coordinates must be postulated. Common choices are described in Sec. 2.3. In later chapters, the choice of collective coordinates will be noted explicitly. 3.1.2 WKB Theory To describe tunneling behavior, we re-quantize the total energy of the system. Doing so results in a multidimensional Schrodinger equation in the collective space. Solving this system for WKB-type wave functions has been formalized in Ref. [112]. However, the tun- 26 neling probability is dominated by isolated one-dimensional trajectories, called instanton trajectories [113, 114]. This is closely related to the path integral formulation, in which the probability amplitude between two points is dominated by classical trajectories [115, 116]. Thus, we search for one-dimensional tunneling wave functions. To re-quantize the Hamiltonian, we introduce a one-dimensional coordinate s, and define the trajectory q(s). The total energy is E(s) = V (q(s)) + 1 2 Meff(s) ˙s2, where Meff(s) dqα ds (cid:88) αβ ≡ Mαβ(q(s)) dqβ ds (3.12) is the effective mass, and the re-quantized Hamiltonian is ˆH(s) = V (q(s)) ℏ2 − 2Meff(s) d2 ds2 . (3.13) We solve this Hamiltonian for WKB-type wave functions [117] ψ(s) = N (s) e− S[s], where S[s] = (cid:90) s 1 ℏ ds′ (cid:112) 2Meff(s)(V (s) E0) − (3.14) is the classical action integral. The normalization factor N (s) is neglected in most treatments. The energy E0 is chosen in this work as the ground-state potential energy, E0 = V (qgs). Far from the classical turning points, where V (s) = E0, this is an approximate eigenstate of the collective Hamiltonian, to lowest order in ℏ. Note that some authors add a zero-point energy correction associated with additional quantum fluctuations [21]. The SF half-life can be estimated as follows: approximate the potential minimum as quadratic, and consider a classical particle bouncing between the barriers of the well. At each contact with the outer point, there is a probability to tunnel through the barrier (and 27 the same probability to tunnel back to the ground state). The half-life t1/2 is then determined by the frequency of assaults on this barrier [118]: t1/2 = ln 2 nP , where P = 1 1 + exp[2S[Lmin]] (3.15) is the penetration probability along the least-action path Lmin and n is the number of assaults on the barrier per unit time. While n may be computed quantum mechanically [118, 119] we take it to be 1/n = 10− 21 s. It remains to determine the trajectory Lmin along which this wave function is valid. Lmin is precisely the trajectory that minimizes the classical action S[L], and as such will be called the least action path (LAP). Before discussing how to compute it, I will present an illustrative example of the PES, together with the computed LAPs. 3.1.3 An Illustrative Example Figure 3.1 shows the PES for 256Fm computed using the D1S EDF. This example considers two collective coordinates, (Q20, Q30), see Sec. 2.3.2. Axial symmetry is enforced, so all moments Qµλ are zero for λ = 0. All other moments are left unconstrained, and hence are minimized over. First, notice the ground state configuration. The density profile is shown; as can be seen, the ground state is compact and reflection-symmetric. The HFB energy is normalized to V (qgs) = 0. Next, notice the white curve. It separates the classically forbidden (V (q) > 0) and classically allowed (V (q) < 0) regions. We refer to this as the outer turning line (OTL). In higher dimensions, it is a hypersurface; nevertheless I will refer to this as the OTL. The computed LAPs are shown in red and green. They connect the ground state to the OTL. The connection points, marked as dots, are called the exit points. The red (green) 28 ̸ Figure 3.1: An example PES computed using the D1S EDF. Taken from Ref. [13]. See the text for a detailed description. path corresponds to symmetric (asymmetric) fission. Density profiles are shown for both exit points; as can be seen, they are reflection-symmetric (reflection-asymmetric). Notice also that the nucleus has not yet split into two fragments; this typically occurs at a much larger elongation than is shown on the figure, at scission points. The identification of scission points has been discussed considerably in recent years [21]. 3.2 Determining the Least Action Path There are a number of different approaches to computing the LAP. Refs. [120–122] assumed a parametric form for the LAP, and minimized over the variational parameters. The eikonal equation has been solved directly, see Refs. [113, 123, 124]. This approach is closely related to solving the Euler-Lagrange (EL) equations for the trajectory q(s), see Ref. [13]. These approaches are, however, prone to numerical difficulties [125, 126]. Grid-based methods have 29 also been widely used, see [109, 127, 128]. And, recently, we have implemented an iterative minimization scheme based on the discretization of the LAP, called the nudged elastic band method [13]. Note that a common approach is to search for minimum energy pathways (MEPs), which follow the minimum energy valley as the nucleus elongates. This is equivalent to the LAP when the collective inertia is assumed constant, hence it is also called the static fission pathway. This has previously been used to study nuclei for which multiple fission valleys appear, see [91] and Ch. 4. As discussed in Ref. [13], it is only under restrictive conditions that the LAP is also an MEP, although for nuclei with well-defined fission valleys, the two pathways tend to be similar. In this section I will first review two common grid-based approaches, then discuss the nudged elastic band approach. This content is drawn largely from my publication, Ref. [13]. 3.2.1 Grid-Based Approaches The PES is, strictly speaking, defined as a discretized grid of calculations, with further evaluations defined via some interpolation scheme. As such, one option for computing the LAP is to compute the path through the gridpoints that minimizes a discrete approximation of the action functional between an initial point, qin, and a final point, qfin. These paths can be defined recursively, following a general strategy known as dynamic programming. When presenting these algorithms, I will assume a two-dimensional uniformly-spaced grid, for which a point is labeled qij = (xi, yj), 1 i ≤ ≤ Nx, 1 j ≤ ≤ Ny. Extensions to higher dimensions and non-uniform grids are conceptually straightforward. One variant of dynamic programming has been applied to LAP calculations in Refs. [109, 129]. This implementation will be referred to as the dynamic programming method (DPM). DPM finds paths that traverse diagonally from a given cell: from cell qij, cells qi+1,j can 30 Figure 3.2: An example of pathways that may be found using the dynamic programming method (in black) and the NEB method (red). The common starting point qin is marked in green. The endpoint qfin is marked in black (red) for DPM (NEB). be reached for all 1 j ≤ ≤ Ny. Example pathways are shown in Fig. 3.2. This may be understood physically if x corresponds to the elongation of the nucleus. It seems reasonable that a fissioning nucleus should only ever increase in its elongation, and hence the LAP should increase monotonically in x. So, paths that backtrack or stay constant in x need not be considered. Still, there are approximately N Nx y such paths, which cannot all be enumerated. Instead, the LAP from qin to qfin is constructed iteratively. For a cell qij, there are Ny possible paths, each passing through a cell at xi − 1. The LAP from qin to qij is selected and stored in memory. This is repeated for every cell with x = xi, for a total of Ny possible paths. Once qfin is reached, there are only NxNy paths, and the LAP is selected from these. The DPM algorithm is detailed in Algorithm 2 of the supplemental material of Ref. [13]. 31 The algorithm is closely related to Dijkstra’s algorithm [130]. The principal difference is the allowed grid paths. While DPM assumes paths that are monotonic in x, it does allow for nonlocal jumps in y, such as from y1 to yNy . Different grid connections (or “stencils”) may be considered. An alternate stencil is nearest-neighbors: paths to qij can go through any point qi ± 1,j ± 1 that is not already part of the path. This stencil is also implemented in Ref. [13]. As will be seen, this stencil tends to predict greater action values than when using DPM, corresponding to the physical intuition mentioned above. 3.2.2 The Nudged Elastic Band The grid-based methods previously described have a number of limitations. First, they scale poorly with both the number of collective coordinates and the number of gridpoints. The former is problematic because, as discussed in Ch. 2, there are many collective coordinates important for fission. The latter is problematic because improving the precision of the lifetime calculations by increasing the density of gridpoints (say, by interpolation) has an exponential computational cost. Second, grid-based algorithms cannot easily take advantage of previously-known information. For instance, the MEP mentioned above is typically close to the LAP, which could help speed up LAP calculations. Third, grid-based algorithms only compute the global LAP. Many nuclei have competing fission modes (see Ch. 4), each of which corresponds to a local minimum of the action integral. Grid-based calculations are unable to describe this phenomenon. The nudged elastic band (NEB) algorithm addresses each of these challenges. It was first implemented to describe MEPs in molecular transitions between reactant and product states [131–133], and refined in Refs. [134, 135]. The NEB algorithm was subsequently mod- ified to describe LAPs in atomic tunneling [136]. We have implemented it in an open-source Python package labeled PyNEB (see https://pyneb.dev/) for use in fission calculations, as 32 part of the publication of Ref. [13]. The NEB algorithm discretizes the path L into a number of points, called images. The location of the i’th image will be denoted as qi. The position of each image is updated according to Newton’s second law, as ¨qi = F net i , as if it were a point mass. The net force i = F net(qi) is split into two terms: F net i = F k F net i + g⊥i . The spring force F k i is defined as F k i = k( qi+1 − | qi| − | qi − qi − )τi, 1| (3.16) (3.17) where τi is the unit vector tangent to the line connecting images i 1 and i + 1 [135]. The − spring force on the boundary images is given in Ref. [13]. The spring force keeps the images from either bunching up or spreading out, hence the images form a band. The parameter k controls the strength of this force term; we take it to be the same for every image. The force term gi is taken either as the gradient of the PES itself, gi = V (qi), or the gradient −∇ of the action functional, gi = ∂S ∂qi − . The former choice corresponds to MEPs, the latter to LAPs. The contribution to the net force is taken as the component of gi perpendicular to the path, g⊥i , to help prevent the image from sliding along the path [137]. One subtlety is that the images at the endpoints need not be fixed. This is especially useful for multidimensional PESs, as the exit point can be determined via the iterative process itself. An example pathway using NEB is shown in Fig. 3.2: note that the endpoint qfin is not fixed to a gridpoint. Typically, only the outer image is left free. Then, an 33 additional harmonic force term is added to pull the image to the outer turning line [13]. Note that this force term also has parameter controlling its strength. The particular choice of force update is important for the convergence rate of this algo- rithm. Originally [131–135], a simple Verlet update was used [138]. While stable, this update algorithm can require many iterations before convergence is achieved, especially for nearly- flat PESs. Later, an inertia-based scheme, combined with an adaptive timestep, called the Fast Inertial Relaxation Engine (FIRE), was introduced [139, 140]. Both schemes are im- plemented in PyNEB, and indeed the FIRE algorithm requires an order-of-magnitude fewer iterations than the Verlet optimizer. Note that these optimizers do introduce additional hy- perparameters. The optimization is largely stable for reasonable choice of hyperparameters, and future studies will provide guidance on the optimal choice. 3.3 Initial Applications We have benchmarked PyNEB by comparison with both analytic surfaces and realistic PESs. A qualitative comparison is obtained simply by looking at the pathways and noting whether they appear similar. A quantitative estimate is obtained by examining the action values along the pathways. For the latter, the grid-based algorithms are an especially good point of comparison, as they provably minimize the action for their choice of stencil. The analytic surfaces are described in detail in Ref. [13]. For these calculations, the collective inertia is neglected. The action values are replicated in Table 3.1. Also shown are the MEP results, labeled NEB-MEP, and the numerical solutions to the EL equations (determined by the shooting method), labeled EL. The DPM tends to determine paths with a lower action than both Dijkstra’s algorithm and the EL equations, the latter due in part to numerical challenges with solving the EL equations. The sole counterexample occurs 34 Table 3.1: Action integrals for the 6-Camel-Back (CB-S and CB-A) and M¨uller-Brown (MB) surfaces (defined in Ref. [13]). The integrals have been calculated using a linear spline interpolation evaluated at 500 points along each trajectory. NEB-MEP NEB-LAP DPM 5.524 6.405 22.909 5.522 6.793 28.491 5.518 6.404 22.875 EL 5.536 6.407 22.871 Dijkstra 5.563 6.886 23.427 CB-S CB-A MB when the LAP requires a non-monotonic trajectory in the M¨uller-Brown surface, which is not expected to occur in fission calculations as discussed above. Note also that the MEP action values are consistently higher than the LAP action values, indicating the poorness of the static pathway for these particular surfaces. Finally, for all cases, the NEB algorithm determines LAPs with a comparable, but slightly lower, action value as compared to DPM. This indicates that the NEB algorithm indeed is able to accurately determine LAPs. Next, two actinide nuclei are considered. Figure 3.3 shows 232U in the 2D collective space (Q20, Q30) using the SkM∗ EDF. As can be seen, the LAPs are all in excellent agreement with each other. Additionally, the LAPs are quite similar to the MEP, due to the well- defined fission valley, in contrast to many of the analytic surfaces. Note also that the exit point also agrees well between the different approaches; the importance of this will be seen in Ch. 4. And, the action integrals along all pathways, including the MEP, agree to within approximately 1% (see Ref. [13] for the numerical values). Once again, this demonstrates that the NEB algorithm is able to determine LAPs comparable to previously-used algorithms. Finally, 240Pu is considered in the 3D collective space (Q20, Q30, λ2), recall Sec. 2.3.3. The PES is shown in Fig. 3.4, starting from the fission isomer. Note that the pathways found using Dijkstra’s algorithm are not shown here; otherwise, the pathways are marked identically to Fig. 3.3. As can be seen, while the LAPs agree well with constant inertia, they disagree quite noticeably when the cranking inertia is considered, primarily in the 35 Figure 3.3: The PES for 232U computed using the SkM∗ EDF. The OTL is shown in white. Solid (dotted) lines mark the LAPs and MEP obtained with the constant (nonperturbative) inertia tensor. The blue, orange, purple and black curves represent the LAPs calculated using the NEB, DPM, EL, and Dijkstra methods, respectively. The green curve is the MEP calculated using NEB. Taken from Ref. [13]. dynamical pairing degree of freedom. Curiously, despite the large difference between the pathways found using DPM and NEB, the action integrals are quite similar. This suggests that NEB is falling into a local minimum, which was confirmed by initializing NEB with the DPM pathway. While a possible suggestion was given to use DPM to initialize NEB, this seems unnecessary in general, given that the action integrals agree at the sub-percent level. Finally, I note that, while the exit points disagree considerably in the λ2 degree of freedom, they agree reasonably well in the (Q20, Q30) degree of freedom. Thus, as will be discussed in the next chapter, the primary fission fragment yields are expected to agree well. 36 0100200300400Q20(b)01020304050Q30(b3/2)336669912121515181821212424272703691215182124 Figure 3.4: The PES for 240Pu computed using the SkM∗ EDF, normalized to the fission isomer energy. The outer turning surface is shown as a blue mesh. The coloring of the paths is the same as in Fig. 3.3. The 2D cross-section at λ2 = 0 is also shown. Taken from Ref. [13]. 3.4 Conclusion In this chapter, I have described the tunneling process in spontaneous fission, beginning from the reduction to collective coordinates, and ending with the practical calculation of tunneling pathways. As has been established, the nudged elastic band algorithm is indeed the appropriate tool for computing tunneling paths. Since the initial publication, I have optimized PyNEB for runtime performance, achieving more than 100-times speedup over the initial implementation, which has greatly aided the calculations done in later chapters. A study of the optimal hyperparameters of the algorithm is in progress. Additionally, PyNEB is currently in use by a number of other research groups, such as the recent publication [141]. 37 AAAB8HicbVDLSgMxFL1TX7W+qi7dBIvgqswUUZdFNy4r2Ie0Q8lkMm1okhmSjFCGfoUbF4q49XPc+Tem01lo64HA4Zxzyb0nSDjTxnW/ndLa+sbmVnm7srO7t39QPTzq6DhVhLZJzGPVC7CmnEnaNsxw2ksUxSLgtBtMbud+94kqzWL5YKYJ9QUeSRYxgo2VHgfcRkM8bAyrNbfu5kCrxCtIDQq0htWvQRiTVFBpCMda9z03MX6GlWGE01llkGqaYDLBI9q3VGJBtZ/lC8/QmVVCFMXKPmlQrv6eyLDQeioCmxTYjPWyNxf/8/qpia79jMkkNVSSxUdRypGJ0fx6FDJFieFTSzBRzO6KyBgrTIztqGJL8JZPXiWdRt27rHv3F7XmTVFHGU7gFM7Bgytowh20oA0EBDzDK7w5ynlx3p2PRbTkFDPH8AfO5w9qS5Ao2AAAB8nicbVBNSwMxEJ31s9avqkcvwSJUkLJbRD0WvXhswX7AdinZNNuGZpMlyQpl6c/w4kERr/4ab/4b03YP2vpg4PHeDDPzwoQzbVz321lb39jc2i7sFHf39g8OS0fHbS1TRWiLSC5VN8SaciZoyzDDaTdRFMchp51wfD/zO09UaSbFo5kkNIjxULCIEWys5Df7Wc2d9i4r4UW/VHar7hxolXg5KUOORr/01RtIksZUGMKx1r7nJibIsDKMcDot9lJNE0zGeEh9SwWOqQ6y+clTdG6VAYqksiUMmqu/JzIcaz2JQ9sZYzPSy95M/M/zUxPdBhkTSWqoIItFUcqRkWj2PxowRYnhE0swUczeisgIK0yMTaloQ/CWX14l7VrVu656zaty/S6PowCncAYV8OAG6vAADWgBAQnP8ApvjnFenHfnY9G65uQzJ/AHzucPqm+QMw==Q20(b)AAAB+nicbVDLSsNAFJ3UV62vVJduBotQQWrSirosunHZgn1AG8NkOmmHTiZhZqKUmE9x40IRt36JO//GaZuFth64cDjnXu69x4sYlcqyvo3cyura+kZ+s7C1vbO7Zxb32zKMBSYtHLJQdD0kCaOctBRVjHQjQVDgMdLxxjdTv/NAhKQhv1OTiDgBGnLqU4yUllyz2HSTmpX2T8vefVI7q6YnrlmyKtYMcJnYGSmBDA3X/OoPQhwHhCvMkJQ924qUkyChKGYkLfRjSSKEx2hIeppyFBDpJLPTU3islQH0Q6GLKzhTf08kKJByEni6M0BqJBe9qfif14uVf+UklEexIhzPF/kxgyqE0xzggAqCFZtogrCg+laIR0ggrHRaBR2CvfjyMmlXK/ZFxW6el+rXWRx5cAiOQBnY4BLUwS1ogBbA4BE8g1fwZjwZL8a78TFvzRnZzAH4A+PzBwJRkos=Q30(b3/2)100150200250Q20(b)051015202530Q30(b3/2)0.51.01.52.02.52.53.03.03.53.54.04.00.01.22.53.85.06.27.58.810.0 Chapter 4. Multimodal Fission As discussed in the previous chapter, the nudged elastic band algorithm can identify pathways that locally minimize the action integral, and hence can be used to study competing fission modes. This is carried out in my publication, Ref. [14], and is the subject of this chapter. In the constrained DFT framework described above, multiple fission modes can be seen as multiple fission valleys in a PES, see Fig. 3.1. These pathways have been identified numerically using minimum energy trajectories [91, 142, 143]. However, fission pathways are not observable quantities. Experimentally, one way to distinguish multiple modes is the primary fission fragment yields (FY), recall the discussion in Ch. 1. In this chapter, I will discuss how FY can be computed, with a focus on estimating them using the exit point location. I will then discuss multimodal fission in the Fermium chain in two- and three-dimensional collective spaces. 4.1 Primary Fission Fragment Yields There are a number of different models used to compute primary fission fragment yields (FY), which may be divided into two schematic approaches: those based on phenomenological approaches to fission, which I have neglected in this thesis so far, and those rooted in the EDF, primarily in time-dependent frameworks. I will briefly review these different approaches, with the goal of motivating the approach used in this thesis and Ref. [14]. I will then discuss the approach we use. 4.1.1 Phenomenological Approaches Despite the focus on DFT described throughout this work, fission has long been understood from a more phenomenological point of view, beginning from its description as competition between surface energy and the Coulomb force. Such approaches are relatively computation- 38 ally inexpensive, allowing for calculations across regions of the nuclear chart relevant to the r process. For a recent review of these approaches, see Ref. [144]. FY are commonly described by assuming dissipative motion outside or above the effective barrier. This can be understood based on the separation of the nuclear Hamiltonian into collective and intrinsic parts. The collective part of the Hamiltonian corresponds to the PES, while the intrinsic part describes excitations on top of the collective deformation. As the nucleus evolves in time, energy can contribute to either degree of freedom. Considering only the collective energy, then, energy that goes into intrinsic degrees of freedom may be interpreted as being lost from the collective degrees of freedom - hence, dissipation [145]. This motivates the Langevin dynamical approach, in which the following stochastic equations are solved [146]: dpi dt dqi dt = pjpk 2 ∂ ∂qi − 1 M − jk − ∂ ∂qi V − = M − 1 ij pj. 1 ηijM − jk pk + gijΓj(t), (4.1a) (4.1b) Here, p is the conjugate momentum to the collective coordinates q, η is the dissipation tensor, g is the force-strength tensor, and Γ is a time-dependent Gaussian random variable. V and M are the PES and collective inertia described above. Via the fluctuation-dissipation theorem, η and g are related to the nuclear temperature, which itself depends on the level density of the deformed nucleus. A common challenge is the specification of the dissipation tensor: while intuition sug- gests its importance, the precise form is not known from microscopic theory. Nevertheless, such approaches have met with varying degrees of success [144]. Of particular note is a global analysis following this approach, the Brownian shape model(BSM) [147], which will 39 be discussed in relation to the Fm chain later. A recent application in which V and M are computed using DFT is presented in Ref. [148]. The location of the FY peaks is some- what insensitive to both the inertia and the dissipation tensor, although the tails depend noticeably on both. Alternative approaches are not based on dissipative dynamics. Two approaches that will be contrasted with later are instead the general description of fission observables (GEF) [149], which is based on a number of fairly general theorems combined with empirical data, and the scission point yield (SPY) [150], which describes fission properties based solely on the nucleus at scission. While these approaches are reasonably successful at and near stability [144], there are qualitative differences for many nuclei, such as the number of FY peaks and their locations [150]. 4.1.2 Time-Dependent Approaches The past decade has seen rapid advancements in modeling FY using time-dependent frame- works. Modeling of FY can take place entirely within the classically allowed region, allowing for direct time evolution that is otherwise forbidden. These frameworks, while considerably more computationally expensive than the above phenomenological approaches, make fewer assumptions as to the nature of the fission process. For instance, they do not assume adia- baticity, and recent time-dependent DFT (TD-DFT) studies suggest the evolution is indeed non-adiabatic [151]. A first approach to FY consists of directly iterating the TD-DFT equations forward in time, beginning from a highly deformed state outside of the barrier and ending when the fragments are well separated. Then, properties such as the excitation energy may be described for each fragment independently. One known problem is that large fluctuations in observable quantities (such as the par- 40 ticle number dispersion from Sec. 2.3.3) are consistently underestimated. A fix appears with the time-dependent generator coordinate method (TD-GCM), a multi-reference DFT extension, in which the wave function is, for instance, a linear combination of HFB wave functions. Evolution of this trial wave function can also be used to compute FY, now with improved correlations. See Ref. [90] and references within for many references relevant for this discussion, as well as a successful application of the TD-GCM approach. Recent work has also extended this approach to odd-mass nuclei [152]. 4.1.3 A Hybrid Approach Direct time-dependent approaches are too computationally expensive for near-term large- scale calculations. Nevertheless, approaches using DFT input are desirable, to avoid many of the assumptions present in phenomenological models. A hybrid approach was developed in Refs. [153, 154], which has been used throughout this thesis. For that reason, I detour to review it here. When examining the nuclear configuration at the exit point, it is convenient to use the nucleon localization functional (NLF) [155, 156], (cid:20) Cq = 1 + (cid:18) τqρq (cid:19)2(cid:21) 1 − 1 4|∇ − ρqτ TF q ρq 2 | , where τ TF q = 3 5 (6π2)2/3ρ 5/3 q (4.2) is the Thomas-Fermi density and is introduced for normalization purposes. The NLF is related to the conditional probability of finding a nucleon at location r′, given that a nucleon has been found at location r. This particular form is normalized to lie between 0 and 1, with 0 (1) representing a low (high) conditional probability [157]. An example is shown in Fig. 4.1. Observe that the nucleus appears to have two well-formed prefragments, connected by 41 Figure 4.1: The localization functional for 258Fm using the UNEDF1HFB EDF at the sym- metric exit point. (a) and (b) are proton and neutron localizations, respectively. Horizontal dashed lines denote the center of the prefragments. a neck of non-zero width. The prefragments are expected to remain well-formed as the nucleus splits. The total number of nucleons assigned to each prefragment is determined by integrating upwards (downwards) in z for the top (bottom) prefragment and doubling the result. In general, this does not add up to the total number of nucleons in the nucleus, and the remaining nucleons are assumed to be in the neck region of the nucleus. These neck nucleons are then distributed between the prefragments according to a microcanonical probability distribution [158]: P (A1, A2) (cid:118) (cid:117) (cid:117) (cid:116) ∝ (cid:32) (A1A2)8 5/3 2 5/3 1 + A )3(A1 + A2)2 (cid:33) a1a2 (a1 + a2)5 (cid:33) E 9/4 r (cid:112) exp[2 (a1 + a2)Er]. (4.3) (A (cid:32) 1 × 1 − (cid:112) 2 (a1 + a2)Er 42 −10−50510r(fm)−15−10−5051015z(fm)(a)(b)0.050.150.250.350.450.550.65 Here, ai = Ai/10 MeV− 1 is the level-density parameter and Er = (Eb1 + Eb2 + EC ) is the − energy of a given fragment configuration [148]. We take the binding energy of each fragment, Eb1 and Eb2, to be the liquid-drop binding energy [159], and the Coulomb energy EC to be that of point charges at the center-of-mass of the fragments. Finally, the mass and charge yields are convolved with Gaussian distributions; for the charge yield, the odd- and even-Z yields are convolved individually to characterize the odd-even staggering. The identification of the prefragments is not entirely a well-posed problem, especially when one of the prefragments appears to be octupole deformed. When the prefragments are not octupole deformed, the center of the prefragment is a local extrema of the NLF along the r = 0 line, and the center is determined by the location of this extremal value [153]. Otherwise, the extremal value does not correspond to the center, and we estimate it following the discussion in the supplemental material of Ref. [14]. There are a number of assumptions built into this approach. Improvements such as microscopic level-density parameters or self-consistent binding energies may be considered. However, a major source of uncertainty is the two-particle uncertainty in the distribution - that is, the number of neck nucleons is generally not an integer. So, the FY is recomputed for the configurations (Ni ± large, and is expected to be larger than the uncertainty from varying Er [153]. Within this 1), i = 1, 2. This source of uncertainty is typically somewhat 1, Zi ± uncertainty, the theoretical FY agree with the experimental values. For these reasons, this FY model was deemed sufficient for this work. Finally, when there are multiple fission modes identified, each trajectory has a probability 43 Pi given by Eq. (3.15). The relative probability ¯Pi of the i-th mode is then ¯Pi = Pi j Pj (cid:80) . (4.4) We assume that individual fission modes are independent decay channels, which is justified in the common case in which a high potential barrier separates the two LAPs. Then, the total FY are determined as Y (A) = (cid:80) i ¯PiYi(A), with Yi(A) the FY along the i-th mode. 4.2 The Fermium Chain The even-even Fm isotopes with 154 N ≤ ≤ 160 are known to transition from asymmetric to symmetric fission as N increases towards N = 164. This transition is related to strong shell effects in the fragments as they approach the doubly magic nucleus 132Sn [160]. Nuclei in this transition region then undergo bimodal fission, which has been investigated in numerous papers [91, 124, 142, 153, 154, 161–174]. For this reason, we chose to investigate this isotopic chain using PyNEB, in Ref. [14]. I note that, in this study, the perturbative cranking approximation to the collective inertia was used, see Sec. 3.1.1. 4.2.1 Two Collective Coordinates First, we consider the nuclei 254,256,258,260Fm with the SkM∗, UNEDF1HFB, and D1S EDFs, in the 2D collective space (Q20, Q30). These PESs are shown in Fig. 4.2. In general, the topologies of PESs agree well between various EDFs. The fission barriers obtained with SkM∗ and D1S are higher than those of UNEDF1HFB. Except for 260Fm with UNEDF1HFB, all the PESs predict a fission isomer (FI) at 100 < Q20 < 160 b and Q30 = 0. Recall that the fission pathway between the ground state and the FI tends to be triaxially deformed, see Sec. 2.3.2. Since triaxiality is ignored in this study, in most cases, one path along Q30 ≈ 0 connects the g.s. and the FI. The reflection-symmetric path continues past 44 Figure 4.2: Potential energy surfaces of 254,256,258,260Fm (in MeV) calculated using UNEDF1HFB (left), SkM∗ (center), and D1S (right) EDFs. The symmetric (dashed lines) and asymmetric (solid and dotted lines) least action paths are drawn from the ground state (filled circle) to the fission isomer (asterisk), to the exit point (open square). The white contour denotes the outer turning line Veff = ∆E = 0. Gray contours are marked at 1 MeV intervals for 0 < Veff < 5 MeV. Taken from Ref. [14]. the FI, and then a bifurcation resulting in a coexistence of symmetric and asymmetric pathways takes place. The PESs tend to be relatively flat in the region beyond the FI, leading to straight path segments. For 258Fm calculated with SkM∗, we observe an additional asymmetric LAP that ends close to the symmetry axis, due to the geometry of the outer turning surface. However, within our hybrid approach, this exit point corresponds to a wide 45 01020UNEDF1HFB5510151520254Fm(a)SkM∗551010151520(e)D1S55101010151520(i)010205510151520256Fm(b)551010151520(f)55101010151520(j)01020551010151520258Fm(c)55101015152020(g)551010151520(k)0801602400102055101520-10260Fm(d)08016024055101015152020(h)080160240551010151520(l)Q20(b)Q30(b3/2) symmetric yield distribution. Figure 4.3: The relative probability of the symmetric mode, ¯Ps, plotted along the Fm isotopic chain. Magenta dots, black s, and green squares correspond to the UNEDF1HFB, SkM∗, and D1S EDFs, respectively. Solid (dashed) lines correspond to 2D (3D) collective spaces. × The relative probability ¯Ps of the symmetric mode is plotted in Fig. 4.3. All EDFs transition from asymmetric-dominant to symmetric-dominant fission path with increasing N . Competition between the modes is predicted in 256,258Fm for UNEDF1HFB, 258Fm for SkM∗, and 256Fm for D1S. This competition is generally consistent with the SkM∗ results of Refs. [88, 91], in which the fission probabilities were computed along minimum energy path- ways. The small differences are most likely due to (i) our inclusion of the multidimensional collective inertia tensor in computing the collective action and (ii) our precise minimization of S(L). Figure 4.4 shows the calculated yields, together with experimental data [175–178]. The experimental data for 254,256Fm show an asymmetric distribution, while for 258Fm, the yields are primarily symmetric, with a small asymmetric shoulder. The error bands in the calculated distributions originate from the two-particle uncertainty, see Sec. 4.1.3. Our 46 254256258260A0.00.20.40.60.81.01.2¯PsUNEDF1HFBSkM∗D1S Figure 4.4: Fission fragment mass (left) and charge (right) yields for 254,256,258,260Fm calcu- lated with UNEDF1HFB (magenta vertical patterns), SkM∗ (black horizontal patterns), and D1S (green patterns). Experimental yields (circles) [175–178] are shown where available. Taken from Ref. [14]. × UNEDF1HFB results are in close agreement with the data for 254Fm. Competition between modes for 256,258Fm is present for each EDF, although the D1S EDF overestimates the symmetric contribution. All EDFs predict overlapping symmetric yields for 260Fm. The transition from asymmetric to symmetric fission is clearly present, albeit at different neu- tron numbers for the EDFs. Despite the overestimation of the symmetric mode for 256Fm, 47 024(a)0510(b)254Fm0510yield(%)(c)01020yield(%)(d)256Fm0510(e)01020(f)258Fm80100120140160180fragmentmass0510(g)3040506070fragmentcharge01020(h)260Fm UNEDF1HFB provides the best description overall. This may be due to the calibration in- cluding the fission isomer energy rather than the fission barrier height, as the former is an observable quantity while the latter is not, see Sec. 2.2. The charge yields using the D1S EDF are largely consistent with the recent TD-GCM calculations of Ref. [171], although the contribution to the symmetric mode for 258Fm is somewhat overestimated in our results. The mass fragment yields agree well with the GEF results [149], and disagree with both the SPY [150] and BSM [147] results. Both references discuss this disagreement. In the former, it is attributed to an underestimate of the contri- bution of two 132Sn fragments. In the latter, a number of model tweaks are suggested, to both the shell-correction approach and the stochastic dynamics used. We have also carried out calculations for 258Rf, 262Sg, and 262Hs, for the UNEDF1HFB and SkM∗ EDFs; the results are presented in Ref. [14]. Similar to the Fm chain, individual fission modes agree well between EDFs, but the relative probability of each mode differs drastically between EDFs. As with the Fm chain, large-scale calculations also make dif- ferent predictions than our approach. The sensitivity of the fragment yields to the EDF parameterization is, therefore, a potentially useful tool for discriminating between various EDFs. 4.2.2 Three Collective Coordinates Next, we include the hexadecapole degree of freedom Q40, see again Sec. 2.3.2. Constraining this moment resolves fission pathways not seen in the (Q20, Q30) collective space, as has been demonstrated in Refs. [70, 91, 142, 179–182]. The D1S EDF is not considered in this section, as it does not predict competing fission modes in the 2D collective space for 258Fm. First, we consider 254Fm, as a nucleus with a single dominant fission mode. With the SkM∗ EDF, an additional elongated fission pathway is identified, consistent with Ref. [91], 48 while for the UNEDF1HFB EDF, no additional pathway appears. Still, the asymmetric pathway is the only pathway that contributes to fission as in the 2D collective space, see Fig. 4.3. The fragment yields are thereby similar between the collective spaces [14]. Similar results are obtained for the UNEDF1HFB EDF. This demonstrates that individual fission modes are somewhat robust to the choice of collective coordinates. Figure 4.5: The least action paths in the 3D collective space for 258Fm using SkM∗ (a) and UNEDF1HFB (b), starting from the FI. The outer turning surface is shown in purple. A 2D PES is shown for constant Q40 = 16 b2 (Q40 = 32 b2) for SkM∗ (UNEDF1HFB). Neutron localization functionals for the identified precission configurations are shown in the insets. Taken from Ref. [14]. Next, we consider 258Fm. The fission pathways, along with the neutron localization at the identified exit points, is shown in Fig. 4.5. Again as in Ref. [91], an additional elongated symmetric fission mode is identified, albeit with relative probability 0. For the ≈ SkM∗ EDF, the geometry of the outer turning surface is similar to the 2D case, leading again to a weakly asymmetric exit point that corresponds (in the hybrid yield framework) to symmetric FY. For the UNEDF1HFB EDF, the asymmetric fission mode is closer to the Q30 = 0 axis, resulting in a more symmetric contribution to the FY. This worsens agreement with the tails of the FY distribution, as shown in Fig. 4.6. The relative probability of the compact symmetric mode is also shown in Fig. 4.3; it corresponds to the symmetric mode 49 identified in 2D. For the SkM∗ EDF, ¯Ps is much larger in 3D than in 2D, resulting in a better agreement with the experimental data. Taken together, this suggests that additional collective coordinates should be considered when multiple fission modes are important, as their addition may induce competition between the modes. Figure 4.6: The fragment mass (left) and charge (right) yields for 258Fm using SkM∗ (top) and UNEDF1HFB (bottom) calculated in the 2D (green solid lines) and 3D (blue dashed lines) collective spaces. Experimental data from [178] is shown as red filled circles. Taken from Ref. [14]. 4.3 Conclusion In this chapter, I have presented results in systems with competing fission modes. Different modes are computable using the nudged elastic band method implemented in Ch. 3, along with the relative probability of each mode. Using an approximation scheme developed in Refs. [153, 154], we have estimated the primary fission fragment yields. When two fission modes compete, these yields strongly depend on the EDF parameterization. They also weakly depend on the number of collective coordinates, although this warrants further study. The UNEDF1HFB EDF is shown to be the most accurate EDF of those studied for the Fm chain. 50 0510yield(%)(a)0102030yield(%)SkM∗(b)70130190fragmentmass0510(c)305070fragmentcharge01020UNEDF1HFB(d) The bulk of my contribution to this work was computing the PESs, LAPs, and fragment yields using the UNEDF1HFB EDF. I have also developed a number of useful software tools for computing PESs and FYs, and have supported my collaborators in their use. Finally, previous experimental [183] and theoretical [184] work has suggested a transition from asymmetric to symmetric fission as temperature increases. Thus, competition between fission modes with increasing temperature is of interest. I am currently using the approach from this chapter to study this competition in a number of different systems. 51 Chapter 5. Neural Network Emulation The approach to fission described in the previous chapters is computationally expensive, primarily due to the HFB calculations that must be carried out. This limits its use in large- scale calculations, which is especially problematic for the rapid neutron capture process (r process). The r process involves thousands of heavy, neutron-rich nuclei that are typically outside of experimental reach. Therefore, fission predictions with quantified theoretical un- certainties are necessary. This, in turn, necessitates emulation techniques. In this chapter I will discuss emulation of the PES and collective inertia tensor across the r process region of the nuclear chart using neural networks (NNs). I will begin with a brief discussion of the r process and the relevance of fission observables. I will then discuss details of the NN used, and present the NN results on both the PES and inertia, as well as (proxies for) observable quantities. The results in this chapter are from my publication, Ref. [15]. 5.1 The Rapid Neutron Capture Process The r process is responsible for the production of roughly half of the abundance of elements heavier than iron. It occurs in explosive stellar environments with high neutron densities. Some types of core-collapse supernovae, as well as neutron-star neutron-star (NS-NS) merg- ers, have long been studied as potential r process sites. Recent observations, especially the NS-NS merger GW170817, have largely confirmed the latter. For recent reviews of the r process, see Refs. [5, 185–188]. In an environment with free neutrons, repeated neutron captures and β− decays trace a pathway through the nuclear chart, resulting in heavier nuclei. For low neutron densities, this path stays close to the valley of stability, leading to the slow neutron capture process (the s process). For high neutron densities, this path pushes outward from stability, into 52 the neutron rich region of the nuclear chart. Eventually, near A 260 [188], fission begins ≈ to compete with neutron capture rates, terminating the r process chain. Given enough free neutrons, the fission fragments capture more neutrons from the environment, repeating the cycle. This is called fission recycling, and is important for explaining actinide abundances [5]. Finally, after the available neutrons are exhausted, fission participates in the decay chain back towards stability, thereby affecting the final abundances. Further, recent studies have shown that fissioning nuclei leave visible signatures in the electromagnetic spectrum released in NS-NS mergers, due to the large amount of energy released in this decay [189, 190]. As such, fission properties such as half-lives (often parameterized by barrier heights) and fragment yields are important inputs for r process reaction network calculations. A number of r process nucleosynthesis studies have been carried out in recent years. Many of them take fission inputs from phenomenological models such as those discussed in Sec. 4.1.1, see Refs. [8, 189, 191, 192]. Due to its computational cost, similar calculations using fission inputs from DFT have been rather limited - for instance, while fission barriers have been taken from DFT, fission fragment yields have typically been taken from statistical models [9, 193]. Fission inputs from DFT calculations are therefore of great interest for r process studies. 5.2 The Neural Network Scheme Machine learning (ML) approaches show promise for speeding up such calculations. ML has been used extensively within nuclear physics, see Ref. [12] for a recent review. Specifically, a number of ML tools have been used to emulate PESs in both quantum chemistry [194– 197] and nuclear physics [198, 199]. While these works have focused on PESs for individual systems, Ref. [200] used a committee of NNs to emulate the PES and collective inertia in 53 the space of axial and triaxial quadrupole deformations. The NN results were then able to predict properties such as the HFB energy to within 500 keV root-mean-square error, along with many of the low-lying states to within a few percent error. This work demonstrated that NNs are capable of emulating PESs. It also demonstrated that learning intermediate terms in a calculation (rather than e.g. just learning the binding energy) may be helpful, and allows for predictions of quantities beyond what the NN was trained on. In Ref. [15], we have trained fully-connected feedforward NNs separately on the potential energy V and the components of the collective inertia M across the r process region of the nuclear chart. We work in the (Q20, Q30) collective space, using the D1S EDF, using the solver HFBAxial. M was computed using the nonperturbative cranking approximation. Each NN takes as input (A, Z, Q20, Q30), specifying the nucleus and deformation, then outputs the value (either V or one of Mµν) at that point. We linearly rescale the NN inputs to lie between zero and one using the same rescaling factor for all data points. Each NN has between 2 and 7 hidden layers, with 200 nodes in the first layer and a decreasing number of nodes in subsequent layers. We use the RELU activation function, and train to minimize the root-mean-square error in the desired quantity. For each variant on the NN depth, we train multiple NNs, forming a committee. The predictions from each committee member are combined in a weighted average, to reduce the error due to stochasticity in the training of any individual NN. As training data, we computed PESs and the collective inertia for 194 nuclei, each on a regular grid of 4 b for 0 Q20 ≤ ≤ 248 b, and 6 b3/2 for 0 Q30 ≤ ≤ 60 b3/2. These nuclei are labeled as either training, combining, or validation. The combining and validation nuclei were sampled from a uniform random distribution, such that no region of the nuclear chart is overrepresented in either dataset; the remaining nuclei form the training set. These 54 different datasets are indicated in Fig. 5.1. For each nucleus, the entire grid is used in the training/combining/validation. The nuclei in the training set are used to train individual NNs, the nuclei in the combining dataset are used to combine predictions from the committee members in a weighted average, and the nuclei in the validation set are used for validation of the NN predictions. The weights for each committee member are chosen to minimize the root-mean-square error on the nuclei in the combining dataset. Most of the nuclei (about 70%) are used for training, with the remaining 30% split equally between the combining and validation datasets. In general, the NN performance is not sensitive to the distribution of training data, provided the NN does not attempt to extrapolate across the nuclear chart. No detailed optimization of the choice of training nuclei was carried out. Recall that the collective inertia depends on the derivatives of ρ and κ with respect to the collective coordinates, see Eq. (3.11). Thus, due to level crossings, M can develop discontinuities and rapid variations, meaning that the tensor components can span many orders of magnitude. To enable the NN to learn the tensor, we consider the eigenvalue decomposition M = U ΣU T . U is the 2 × 2 matrix of eigenvectors and Σ is the diagonal matrix of eigenvalues. Since U is an orthogonal matrix, we parameterize it by the Euler angle θ. The NN is then trained on θ and the log of the eigenvalues. Training on this representation of the tensor is similar to normalizing the network inputs, as both put NN inputs/outputs on a similar scale. Additionally, this forces the tensor predictions to be positive semi-definite. 5.3 Neural Network Quality In this section, I discuss the performance of the NN on both the PES and collective inertia, and on the observable quantities that have been discussed in previous chapters. Through- out, I will refer to quantities computed using DFT (the NN) as the reference data (NN 55 reconstruction). 5.3.1 The Emulated Quantities For a single nucleus, we define the root-mean-square error (RMSE) ∆V (A, Z) in energy over the collective domain considered as ∆V (A, Z)2 = 1 n (cid:88) Q20,Q30 [V DFT(Q20, Q30, A, Z) − V NN(Q20, Q30, A, Z)]2, (5.1) where n = 693 is the number of gridpoints evaluated in the PES. Figure 5.1 shows ∆V (A, Z) across the region of the nuclear chart considered, for the deepest NN (7 hidden layers, with 200-175-150-125-100-75-50 hidden units), with rescaled inputs. As can be seen, for most nuclei, ∆V (A, Z) ≲ 0.5 MeV. Exceptions occur, with most remaining below 1.5 MeV. Figure 5.1: ∆V (A, Z) (in MeV) for the deepest NN. The different shapes indicate which dataset each nucleus belongs to. Taken from Ref. [15]. For some nuclei (e.g. 308Cf, 314Fm, and 318No), relatively poor performance is not surprising: these nuclei are on the outer edge of the region of the nuclear chart considered, so the NN is extrapolating away from the training region. For other nuclei (e.g. 232Th, 280Cm), poor performance is unexpected: these nuclei are surrounded by training nuclei, and so should be emulated fairly well. 56 120140160180200220N90100110ZTrainingCombiningValidationTrainingCombiningValidation0.000.250.500.751.001.251.50 To understand this, Fig. 5.2(a) shows the reference PES for 280Cm, and Fig. 5.2(b) shows the difference between the reference PES and its NN reconstruction. This nucleus is chosen because it has ∆V = 2.15 MeV, the largest of all nuclei in the validation set. The difference is less than 1 MeV across most of the PES, notably in the region relevant for fission (shaded in gray). The energy difference is large elsewhere, with a difference of more than 5 MeV, explaining the ∆V value. We conclude that even for nuclei with larger RMSE, NNs could provide a very reasonable description of the fission path. This aspect will be examined further in Sec. 5.3.2. Figure 5.2: (a): the reference PES for 280Cm in MeV. (b): the difference between the reference and reconstructed PESs, again for 280Cm, in MeV. The reference ground state and the reference LAPs are marked with a × symbol and black lines, respectively, in both panels. The energy range 0.5 12 MeV is shaded in gray in panel (b). Taken from Ref. [15]. V DFT ≤ ≤ To assess the sensitivity of our results with respect to the NN architecture, we repeated our calculations employing different NN sizes and rescaling the inputs. The results are shown in Fig. 3 of Ref. [15]. As is generally expected, the training dataset has a monotoni- cally decreasing RMSE as the NN depth increases, due to the increasing number of tunable parameters in the NN. On the other hand, the RMSE for the combining and validation sets is 57 0100200Q20(b)0102030405060Q30(b3/2)ReferencePES5510101015202525(a)0100200Q20(b)VDFT−VNN(b)−5051015202530−5−4−3−2−1012345 fairly stable with respect to the number of hidden layers of the NN. A general improvement is observed when rescaling the inputs (A, Z, Q20, Q30) to be between 0 and 1. This is because (i) the NN weights are not scale-invariant [201], and (ii) neither is the optimization algo- rithm [202, 203]. We conclude that the NN performance in predicting the PES is relatively stable with respect to the NN architecture; Sec. 5.3.2 will demonstrate that performance on this level is adequate for predicting SF observables. Next, I consider the collective inertia tensor. Figure 5.3 shows the reference inertia components, plotted against the NN reconstructions, for all nuclei considered. The NN used is the 7-layer NN with rescaled inputs, with the number of hidden units as described above. The diagonal components M22 and M33 are predicted fairly well, with distributions in alignment roughly along the diagonal. It is worth noting that the NN tends to underpredict relatively large values and overpredict relatively small values, meaning that the NN is slightly biased towards the mean value of the inertia. The off-diagonal component M23| | is not aligned along the diagonal. This is because this component actually varies across almost 10 orders of magnitude (compare to the 4 orders of magnitude for M22 and M33), and so the NN is biased towards predicting the larger values more accurately. In terms of the angle θ that is actually determined by the NN, it is difficult to predict both small and large angles, and because θ is allowed to be negative, a logarithm transform is not possible. Nevertheless, we obtain a reasonable-looking distribution above M23| | ≳ 10− 4 MeV− 1 b− 5/2, indicating that some learning has indeed taken place. And, the poorly-learned values below 10− 4 MeV− 1 b− 5/2 are truncated at values 10− 6 − 10− 2 MeV− 1 b− 5/2. When changing the depth of the NN, performance is similar. For shallow networks, predictions on the training dataset show a larger bias: the larger (smaller) reference values 58 Figure 5.3: The reference components of M , plotted against the NN reconstructions, for all nuclei considered. The black line is the diagonal, M DFT µν . The blue squares/green cir- cles/orange triangles correspond to the training/combining/validation datasets. Note that, for use on a log scale, we plot the absolute value of M23 (the other components are nonnega- 3. Taken from 1b− tive). M22 is in MeV− Ref. [15]. 5/2, and M33 is in MeV− 2, M23 is in MeV− µν = M NN 1b− 1b− are underestimated (overestimated), to a greater degree than with the deepest network shown in Fig. 5.3. The validation dataset is aligned similar to the deepest network. As the depth of the network is increased, the training data points are aligned closer with the diagonal. This is indicative of the NN tending to overfit on the training data as the number of variational parameters increases. The distribution of M23 values remains approximately the same when increasing NN depth. Overall, the NN performance on the validation dataset is mostly stable when varying the NN depth. The question is whether this performance is sufficient for predicting observable quantities of interest. As with the PES, this question can be directly answered by looking at NN predictions of physical observables. Thus, we consider observable 59 quantities in the next section. 5.3.2 Relation to Observables While encouraging, the results discussed above do not give a perfectly clear estimation of the performance of the NNs. For instance, the NN reconstruction of the PES for 280Cm may be adequate for reproducing SF observables despite the poor RMSE, since the largest deviations occur at deformations that will not be explored by LAPs. Similarly, the NN commonly fails to reproduce the off-diagonal component of the collective inertia, M23, but primarily for small values of M23. We consider two quantities: the lifetime-weighted exit point, as a proxy for the fragment yields, see Ch. 4, and the half-life of the nucleus. For both quantities, we compare three sets of data: the quantity computed using (i) the reconstructed PES and the identity inertia; (ii) using the reference PES and the reconstructed inertia; and (iii) the reconstructed PES and inertia. In this way, we isolate the impact of the PES and inertia emulations, and combine them to assess the overall error of the emulator. In this section, we use the 7-hidden-layer NN with rescaled inputs, with a number of hidden units described above. Based on the relative insensitivity to the depth of the NN, the overall performance is expected to be similar for different NN depths. Figure 5.4 shows the difference in the octupole moment of the lifetime-weighted exit point, for configuration (iii) mentioned above. The octupole moment is chosen because it is critical for explaining multimodality in SF, recall Sec. 2.3.2. The Q30 error is similar for the other configurations, and the quadrupole moment is typically within 1 b for all ± configurations. As can be seen, the Q30 error is mostly 1 b3/2, which we expect results in ± similar FY (within the hybrid method of Sec. 4.1.3). The agreement is mainly due to the accurate PES reconstruction, as previous studies have shown that the exit point location 60 Figure 5.4: The Q30 component of the reconstructed lifetime-weighted exit point, minus the reference Q30 component, in b3/2. These results were computed using configuration (iii), i.e. the NN was used to reconstruct both the PES and the collective inertia. Taken from Ref. [15]. is fairly robust with respect to variations in the collective inertia [72, 109, 153, 204]. This agreement holds even for nuclei whose PES reconstruction has a large error, such as 280Cm, indicating that the qualitative features shown in Fig. 5.2(a) are indeed reconstructed well enough to describe multimodality in SF. Notice that the exit point locations are not reproduced perfectly for some nuclei, espe- cially in the thorium (Z = 90) chain, where the difference can be as much as 5 b3/2. This is not due to the PES reconstruction: Fig. 5.1 shows that the thorium isotopes have RMSE ∆V (Z = 90) ≲ 100 keV, and the exit point reconstruction error of configuration (i) is within 1 b3/2. While a side-by-side comparison of the components of M does not show a systematic deviation between the reference inertia and the NN reconstruction, M is indeed the cause of the error. Random reconstruction error is present for every deformation considered, and it is the accumulation of this random error that causes the discrepancy. While the location of any individual exit point is not sensitive to the collective inertia, the probability of tunneling to a particular point, Eq. (3.15), is exponentially dependent on the action (and therefore M ), small errors can accumulate and switch the dominant exit point from asymmetric to 61 120140160180200220N90100110Z−5−3−1135 symmetric and vice versa. This is especially important for nuclei with a wide fission barrier, as the cumulative error along the path is large. Nevertheless, we observe that both the PES and the collective inertia are emulated well enough to predict exit points that agree with the reference data. For most nuclei, the dominant fission mode is also in agreement. Together, this means that the FY largely agree between the reference data and the NN reconstruction. Figure 5.5: The half-life predicted using the DFT reference data, tsf-DFT half-life computed using the NN reconstruction, tsf-NN . Panel (a) shows configuration (i), in which only the PES is emulated; panel (b) shows configuration (iii), in which both the PES = tsf-NN and the collective inertia are emulated. The black line marks the diagonal: tsf-DFT . 1/2 3, i.e. 3 orders of magnitude above and below the Gray bars are drawn at tsf-DFT 1010 s, to highlight the relevant r-process range. diagonal. Taken from Ref. [15]. Insets show the range 10− , plotted against the 10± 5 1/2 1/2 1/2 1/2 − × Next, I consider the SF half-life. Figure 5.5 shows tsf 1/2 computed using the reference data vs. tsf 1/2 computed using the NN reconstruction, for configurations (i) and (iii) mentioned above (results for configuration (ii) are similar to those of (iii)). As can be seen, the tsf 1/2 predictions agree well, typically within 3 orders of magnitude across the approximately 80 orders of magnitude under consideration. 62 −200204060−200204060log10[tsf-NN1/2(s)](a)−200204060(b)log10[tsf-DFT1/2(s)] Figure 5.5(a) demonstrates that the PES reconstruction is sufficient to predict tsf 1/2 values that agree well with the reference values. As with the FY, this is true even for nuclei with a large ∆V , e.g. 280Cm, once again demonstrating that the PES emulation quality is indeed sufficient to reproduce SF observables. Figure 5.5(b) includes the collective inertia emulation. As can be seen, the reproduced tsf 1/2 values agree less well with the reference values, although the disagreement is still within 3 orders of magnitude for most nuclei. This is not unexpected: as was discussed above, the collective inertia emulation is not accurate enough for all nuclei (although it is sufficient for most). As with the exit points, the disagreement in tsf 1/2 is due to the accumulation of random errors when the fission pathway goes across the fission barrier. Now, rather than changing the dominant fission mode, tsf 1/2 is simply changed from the reference value in a more-or- less random manner. The effect is most prominent for long-lived nuclei, where errors in the collective inertia add up to a fairly large value as the pathway traverses a wider fission barrier. This is demonstrated in Fig. 7 of Ref. [15], which plots the ratio of the half-lives against the barrier width, ∆Q20. As ∆Q20 increases above ∼ 75 b, the difference between reconstructed half-lives and the reference half-lives tends to increase. While it may be desirable in principle to improve the emulation, the nuclei whose tsf 1/2 values are reproduced with a large error are those predicted to be stable to SF, within the (Q20, Q30) collective space. As such, errors in the SF observables have little effect on results that are further dependent on tsf 1/2, such as r process network calculations. 5 Finally, the inset panels in Fig. 5.5 magnify the range 10− 1010 s, to highlight the − relevant r process range. As can be seen, almost all nuclei within this range are reproduced within three orders of magnitude. Therefore, we conclude that NNs are able to reproduce both the PES and the collective inertia well enough that tsf 1/2 is reproduced within 3 orders 63 of magnitude for nuclei for which SF is relevant in the r process region. 5.4 Conclusion In this chapter, I have discussed the emulation of the PES and collective inertia tensor across the r process region of the nuclear chart. As has been shown, neural networks are able to emulate both quantities well, in the deformation region relevant for SF - for instance, the PES is emulated to within ≲ 1 MeV in the barrier region. I have then demonstrated that this accuracy translates to reasonable accuracies in the observables discussed in this thesis, namely the half-life and the primary fission fragment yields. There are a number of desirable improvements with this approach. The first is that the emulated quantities only predict the deformation of the nucleus at the exit point. To actually compute FY, additional HFB calculations are necessary. The second is the amount of training data used. As can be seen in Fig. 5.1, a large fraction of the r process nuclei considered are used to train the NN, necessitating tens of thousands of HFB calculations. It is desirable to reduce this number as much as possible, especially since it scales poorly with the number of collective coordinates used. Third, while the emulator error is indeed adequate for SF observables, it is inadequate for uncertainty quantification efforts. For example, propagating uncertainties from the UNEDF1 calibration leads to statistical uncertainties in the fission barrier on the order of 1-2 MeV [11]. Additional emulator error of ≲ 0.5 MeV thereby results in considerably wider uncertainty bands, reducing the apparent reliability of the model. The next chapter introduces a category of emulators that addresses these issues. 64 Chapter 6. Reduced Order Model Emulation We are interested in computing theoretical uncertainties using Bayesian methods. A Bayesian approach to model calibration explicitly encodes prior-known information about, and results in a posterior distribution over, the model (EDF) parameters [205]. Theoretical uncertainties are then propagated by sampling over the posterior parameter distribution and computing the quantity of interest for each parameter sample. For a review of Bayesian statistical approaches, see e.g. [206]. Bayesian methods involve sampling over model parameters, re- quiring tens- to hundreds-of-thousands of model calculations [10]. To address the resulting computational cost, emulators are required. In this chapter, I will discuss a class of emulation techniques known as reduced order models. I will then present a pair of illustrative examples, and conclude with my recent applications to configuration-space DFT solvers. 6.1 Reduced Order Models Reduced order models (ROMs) have been used widely both inside and outside of nuclear theory. Literature on this subject is commonly called dimensionality reduction [207] or model order reduction [208]. In the nuclear physics context, eigenvector continuation [209, 210] is closely related to the reduced basis method [10, 211], which itself falls under the category of ROMs. The recently-developed parametric matrix models [212] similarly fall under this label. Even neural-network-based approaches are often included in this category. A number of textbooks on the subject have been published [207, 213, 214]. For a review article on their use in nuclear theory, see Ref. [210]. I will not attempt to review ROMs in general, or even their applications in nuclear physics. Instead, I will simply describe the common philosophy of such methods. 65 The methods can be abstractly described as follows: consider a completely generic equa- tion of motion F [z; α] = 0, where z, F [z] RK and α ∈ ∈ RM . z is a vector of independent variables, and α are some model parameters. For the Schr¨odinger equation for the simple harmonic oscillator, for instance, z may be the wave function discretized on a coordinate mesh and α may be the harmonic oscillator width. The goal is to then solve F [z; α] = 0 for a number of different values of α. K is typically large, varying between 100-1000 for a one-dimensional potential (see the example in Sec. 6.2.2). The number of independent variables grows exponentially with the dimensionality of the problem - if 100 mesh points are used in 1 spatial dimension, a naive approach uses 1002, then 1003, for two and three dimensions, respectively (modern solvers commonly exploit pseudospectral methods [215], but the point stands). The solution for a parameter α1 tells us some information about the solution for α2, which tells us additional information about α3, and so on. So, one should be able to solve F [z; αm] = 0 faster as m increases, using some information from solutions obtained for i = 1, . . . , m − 1. A ROM is then a generic label for any approach that uses the information from previous solutions to approximate the next solution. 6.2 Illustrative Examples For illustrative purposes, I will first present the one-dimensional Schr¨odinger equation. I will then present results for a self-consistent problem, still in one dimension. 6.2.1 The Schr¨odinger Equation The time-independent Schr¨odinger equation for a potential Vα(x) can be written as Hα ψα | ⟩ = ϵα ψα | ⟩ , where Hα(x) = ℏ2 − 2m d2 dx2 + Vα(x) (6.1) 66 is the Hamiltonian in position space (generalization to a nonlocal potential is straight- forward). Grid-based solutions discretize the Hamiltonian onto a coordinate mesh x = (x1, . . . , xN ), and diagonalize the resulting matrix; I will refer to this as the coordinate space approach. Another common approach is based on the Ritz variational principle. One expands ψα | ⟩ in a basis ϕi⟩} {| n i=1, as ψα | ⟩ ≈ n (cid:88) i=1 ai(α) . ϕi⟩ | (6.2) After writing the Schr¨odinger equation in this basis, one projects onto for j = 1, . . . , n, ϕj| ⟨ arriving at the system of equations n (cid:88) i=1 ai(α) Hα ϕj| ⟨ ϕi⟩ | = ϵaj(α), j = 1, . . . , n. (6.3) This is an eigenvalue problem for the expansion coefficients a = (a1, . . . , an), and can be solved straightforwardly. I will refer to this as the configuration space approach. Both approaches are arbitrarily precise as N, n increase. The total runtime is dominated by the diagonalization time, which is of order (N 3) and O O (n3) respectively. However, the configuration space approach has an additional numerical cost: the numerical integration to compute the matrix elements Hα ϕj| ⟨ ϕi⟩ | , which scales for a local potential as (an approach using matrix-vector products for the integration actually scales as O (n2N ) O ((nN )2)). Thus, even for efficient numerical integration, the coordinate-space approach is faster than the configuration-space approach, as demonstrated in the example in Sec. 6.2.2. To understand where possible speedups may be present, I show two cases. First, consider Vα(x) = αW (x). The simple harmonic oscillator is the example W (x) = x2. Then, the 67 expansion in the basis results in ϕi⟩ | (T + αW )a(α) = ϵa(α), Tij ≡ (cid:28) ϕi (cid:12) (cid:12) (cid:12) (cid:12) ℏ2 − 2m d2 dx2 (cid:12) (cid:12) (cid:12) (cid:12) ϕj (cid:29) , Wij ≡ ⟨ ϕi| W (x) . ϕj⟩ | (6.4) Note that the matrices T, W can be computed independently of α. This is the key sim- plification: as the parameter α is changed, the diagonalization problem is faster than in coordinate space (provided, as is typical, that n < N ). Explicitly, this precomputation is possible because of the affine parameter dependence in the Hamiltonian. Next, consider a Woods-Saxon type potential, Vα(x) = 1 1 + ex/α . (6.5) This is used prolifically in nuclear physics as a simplified single-nucleon potential [31]. The parameter dependence is non-affine, prohibiting the decomposition above. Empirical inter- polation (EI) [216–218] solves this problem by the decomposition Vα(x) ≈ nei (cid:88) j=1 bj(α)Wj(x). (6.6) With this affine decomposition, the matrices ϕi| ϕk⟩ Wj(x) | ⟨ can be precomputed, and the speedup is retained. To compute Wj(x), we evaluate Vαl(x) for a number of parameters αl} { A l=1 on x, then form the matrix Vij ≡ Vαi(xj). The SVD of V then gives an orthonormal basis informed by the snapshots. Denote the i-th singular value of by Σi. The expansion coefficient V bj(α) is expected to be of order Σi for all α within some range of the sample solutions. This is because this approach projects onto the linear subspace that best describes Vij, and 68 hence data similar to the sample solutions. Note that this range may be fairly large [211]. We choose some threshold εei, and only include basis vectors where Σi/Σ0 ≥ corresponds to an error in the basis truncation of order (εei), and hence an error in the εei. This O final result on the same order. The singular values tend to decay exponentially [219]. So, only a few basis functions are necessary (i.e. nei is small, typically nei ≪ N ). To determine the coefficients bj(α), we enforce equality on some subset of the gridpoints (called the collocation points), Vα(xi) = nei (cid:88) j=1 bj(α)Wj(xi), i 1, . . . N ∈ { . } (6.7) This (in general, overdetermined) system can be approximately solved using the Moore- Penrose inverse. That is, define the matrix Wij ≡ Wj(xi) and the vector vi ≡ Vα(xi). This is then a linear system W b(α) = v, whose approximate solution is b = ( W T 1 )− W W T v. An alternate solution is to evaluate Wj(xi) on exactly nei gridpoints, chosen using the MaxVol algorithm [220, 221] to maximize the information contained in the matrix Wj(xsi), i, j = 1, . . . , nei. Both methods lead to the result bj(α) = n (cid:88) i=1 MiVα(xi), (6.8) with the upper bound being n = N (n = nei) for the overdetermined (MaxVol) case. The advantage of the latter approach is that Vα only has to be evaluated on nei meshpoints, which is especially useful for self-consistent potentials, see Sec. 6.3.3. It remains to choose the wave function basis ϕi⟩} {| . A generic basis, such as the harmonic oscillator, is a common choice. However, in many cases, a basis with fewer elements can be 69 determined. Analogous to determining the expansion functions Wj(x) { , we take the SVD } of the wave functions as α is varied. Note that this introduces another threshold εRBM that must be chosen to truncate the SVD. This approach is what is often referred to as the reduced basis method (RBM): it is simply a basis expansion approach, in a tailored basis. This approach has been applied to self-consistent solvers [10, 211] and widely applied in other areas of nuclear theory; see the recent review [222]. When multiple eigenstates are needed, one may either construct a basis for each wave function separately, or construct a single basis that describes all wave functions simultaneously. The wave function construction will be discussed further in Sec. 6.3. 6.2.2 A Modified Gross-Pitaevskii Equation The Gross-Pitaevskii equation was introduced originally to study superfluidity in bosonic systems [223]. There, the bosons are assumed to have only a contact interaction, plus perhaps an external confining potential V (r). In the Hartree-Fock approximation, which for bosons states that the full wave function Ψ(r1, . . . , rA) = (cid:81)A i=1 ψ(ri), the single-particle state ψ(r) obeys the equation (cid:20) ℏ2 2 + V (r) + α − 2m∇ (cid:21) 2 ψ(r) | | ψ(r) = Eψ(r). (6.9) This involves the self-consistent potential ψ(r) | 2, similar to the EDFs considered above, | motivating its choice as a toy model. To render this model more similar to the Skyrme EDF, we replace the self-consistent potential with ρσ(r), where 2 ψi(r) | | ρ(r) ≡ A (cid:88) i=1 70 (6.10) is the local single-particle density and σ is a (non-integer) tunable parameter. We then solve a coupled system of equations for each wave function ψi(r), i = 1, . . . , A. The use of the density ρ is required to describe fermions rather than bosons. The choice of σ as a free parameter mimics the similar term in the Skyrme EDF, see Sec. 2.2.1. Finally, we specify an external harmonic oscillator potential. Explicitly, in this section we solve the system ˆHGPψi(x) = ϵiψi(x), ˆHGP ≡ − d2 dx2 + kx2 + qρσ, i = 1, . . . , A, (6.11) and I will refer to this equation as the Gross-Pitaevskii (GP) equation. The non-integer power introduces complications. To see this, consider first the case with integer σ, which has been studied in Ref. [211]. Wave functions appear in integer powers, and hence the multinomial theorem guarantees that the basis expansion matrices can be precomputed. Thus, the GP equation can still solved entirely within the reduced space defined by the expansion coefficients. This fails when non-integer σ are considered. I note that Ref. [211] nevertheless obtains a meaningful computational speedup with non-integer σ, albeit with the one-dimensional Skyrme EDF instead of ˆHGP. As with the Woods-Saxon potential, the solution is EI. Specifically, we expand the entire self-consistent potential f (α; x) ≡ qρσ(x) nei (cid:88) j=1 ≈ bj(q, σ)Wj(x), (6.12) with the Wj} { defined via the SVD. An apparent snag is the determination of b(q, σ). With the Woods-Saxon example, the potential is given explicitly, so one may evaluate it at any pair (q, σ; x) so desired. Here, that is not the case. Recall instead the iterative 71 Figure 6.1: A schematic diagram of solving the Gross-Pitaevskii equation using empirical interpolation. The index ℓ is the self-consistent iteration. process described in Sec. 2.1, displayed schematically in Fig. 6.1. So, the potential used to determine b(ℓ+1) may be computed using the wave functions (ℓ) i } . The iterative scheme is ψ { then deemed convergent when max b(ℓ) | − b(ℓ+1) | ≤ εconv for subsequent iterations (ℓ, ℓ + 1), for a specified tolerance εconv. To contrast this approach, I first compare to three different solvers: a harmonic oscillator (HO) basis expansion, a uniform finite-difference grid discretization, and a Chebyshev [215] grid discretization. All approaches use a linear mixing scheme [95]. In the HO expansion, the number of oscillator shells is varied from 20 150. The grid discretization schemes − increase the number of gridpoints used, between 100 1000 and 20 − − 200 for finite difference and Chebyshev grids, respectively. As these values increase, the error on the eigenvalues decreases; conversely, the total runtime increases, as seen in Fig. 6.2(a) for A = 5. As discussed above, the HO expansion option is slower than the finite difference solver, up to a crossover runtime of about 1 second. Additionally, the discretized grid based on Chebyshev meshpoints is faster than the finite difference solver at low accuracy, and scales similar to the HO expansion at high precision. Note that these three solvers use the same tolerance 72 Figure 6.2: The root-mean-square error on the eigenvalues computed using a number of different solution methods. (a) considers the first 5 wave functions; (b) considers the first 20. See the text for details. εconv = 10− 8, and the exact eigenvalue is taken from the highest-precision calculation. Thus, the runtime difference as e.g. the number of HO states increases is due to the increased matrix size, not the total number of iterations. Also shown is the EI implementation with both a HO and a reduced wave function basis (RBM) that describes all A wave functions simultaneously. The EI coefficients are reconstructed in the overdetermined case. The EI (and RBM) tolerance are varied between 3 10− − 10− 7, which corresponds to increasing the number of basis functions. The convergence tolerance is set to εconv = 10 × max (εei, εRBM), as a smaller tolerance is below the truncation 73 10−710−510−310−1(a)HOExpansionFinite-DifferenceChebyshevEIEI+RBM10−310−210−1100101Runtime(s)10−710−510−310−1(b)EigenvalueRMSE error from the expansions. As can be seen, EI in the HO basis is consistently 10 times faster than the Chebyshev approach, and 100 times faster than the HO expansion, for equivalent precision. That EI is only one order of magnitude faster than the Chebyshev grid suggests that pseudospectral methods on tailored grids may provide an additional speedup. The EI + RBM approach demonstrates that an additional factor of 5 may be gained with similar precision. Figure 6.2(b) is the same as Fig. 6.2(a), except for A = 20. Many of the same conclusions hold. The notable exception is that the speedup from the tailored wave function basis is drastically decreased. This is because the wave function basis size must be enlarged to capture higher-lying states. This suggests that a wave function basis becomes less useful as more states are desired; we will return to this point in the next section. Overall, this example demonstrates the effectiveness of the described approach. Next, I will focus on a realistic case with the HFB equations. 6.3 The Axial HFB Case Now, we consider the HFB equations, assuming axial and time-reversal symmetry. A com- mon solution method is to expand in a (transformed) harmonic oscillator basis, as in the solver HFBTHO [33]. Below, I will comment on some of the subtleties that differentiate this from the GP equation, then discuss EI for a representative superheavy nucleus, 254Fm. 6.3.1 Subtleties with the HFB Equations The HFB equations differ considerably from the GP equation discussed in the previous section. The first consideration is the self-consistent potential: in addition to two species of particles (protons and neutrons), there are multiple local densities and currents that contribute to the HFB matrix, and each contributes a term to the self-consistent potential. 74 For each local density O ∈ { ρq, ∆ρq, τq, Jq, . . . q=n,p, } ∇ · ˜ O ∈ { ˜ρq, . . . } q=n,p (6.13) the contribution to the HFB matrix can be written as hOαβ = (cid:90) 0 ∞ dr r (cid:90) ∞ −∞ dz f Oαβ(r, z) δEHFB δ O (r, z) (6.14) for a particular function f Oαβ(r, z), and analogous for the pairing channel h ˜ O. Throughout, the variations (r, z) will be referred to as fields [36]. For the particle densities, the δEHFB δ O fields are related to the central, spin-orbit, and the effective mass [38]. It is, however, convenient to separate the ρq and ∆ρq terms. Then, the full HFB matrix splits into a number of terms: (cid:88) h = hO, ˜h = (cid:88) ˜ O. h O ˜ O (6.15) When split by nucleon species, the sum is over proton/neutron variations (e.g. hn only includes variations with respect to neutron densities). Each of these fields must be recon- structed at every iteration, and the collocation points should be chosen so as to not repeat calculations. In this work, we use the set-wise union of the collocation points for each field. The self-consistent solver then determines the expansion coefficients for each field, b (α) and O (α). Otherwise, the iterative scheme follows that of Sec. 6.2.2, using Broyden mixing [95]. b ˜ O Next, consider the wave functions. An appealing option is to build an emulator for each wave function individually, each requiring few basis states, and hope to achieve a speedup that way. However, this requires tracking wave functions as the EDF parameters 75 are changed. In the GP example above, the number of nodes in the wave function determines which eigenstate is considered (e.g. the ground state or the first excited state) by the well- known Courant nodal domain theorem [224]. Thus, tracking of states is possible. In the HFB case, this theorem no longer holds. Alternatively, one may use clustering techniques to search for similar-looking wave func- tions. Far from (avoided) level crossings, this is indeed possible. Near the crossing, the two wave functions occupy the same subspace and thereby look drastically different from their values far from the crossing. This can be seen in simple two-level systems. The wave functions can be tracked through the crossing using by calculations near the level crossing, but crossings are ubiquitous. This makes the number of required calculations large, reducing the utility of this approach. The only remaining option is to build a reduced basis for all of the wave functions simultaneously. Typical calculations for superheavy nuclei have many more wave functions than nucleons, so this option is far from guaranteed to speed up the calculation. 6.3.2 The Setup In Ref. [11], statistical uncertainties from the UNEDF1LN calibration were propagated to a number of observable quantities using Gaussian Process emulators. As part of that work, 1000 sample EDF parameters were obtained from the posterior distribution. For each set of parameters, we use HFBTHO to compute exact solutions to the HFB equations. We take 100 of these parameters, sampled randomly, as the “training” data. We have checked, and these parameters are evenly distributed throughout the posterior distribution, i.e. we are interpolating in parameter space. We then follow the procedure discussed above to generate EI basis functions, as well as a tailored basis for the single-quasiparticle wave functions. The emulator has been written in Python. 76 Here, I will specifically consider the unconstrained ground state of 254Fm using 25 oscil- lator shells, at the HFB level. Results when varying the number of major harmonic oscillator shells from 20-30, considering LN instead of HFB, computing the unconstrained fission iso- mer instead of the ground state, and studying 236Pu are similar. Figure 6.3: The singular values of the fields for the unconstrained ground state of 254Fm using 25 oscillator shells. The legend displays the field being varied, e.g. ρ corresponds to δE δρ an example threshold εei = 10− . Solid (dashed) lines correspond to protons (neutrons). The dashed black line indicates 5. The singular values of the fields SVD are displayed in Fig. 6.3. As can be seen, the kinetic field δE δτq decays rapidly - for the sample threshold, a mere 10 samples are required. This is due to the constant kinetic energy, Eq. (2.14), whose variation gives a constant ℏ2/(2mq). Conversely, the singular values for the pairing field δE δ ˜ρq decay more slowly, requiring 50-60 sample points for this threshold. They decay more rapidly first, such that a threshold of 10− 4 requires approximately 20 functions. Thus, the threshold chosen sets a minimum on 77 020406080100i10−810−610−410−2100Σi/Σ0ρ∆ρτdivJ˜ρ the total number of data points for which exact calculations must be carried out. Figure 6.4: Panel (a): a sample field δE δρp . Panels (b)-(d): the first three basis field for the variation with respect to ρp. bz and b ⊥ are the oscillator lengths defined in Ref. [36]. For some intuition, Fig. 6.4 shows a sample field δE δρp , along with the first three basis fields. Note that the basis fields are normalized to one. Observe that the basis fields all share the same shape as the sample, and decay to zero for similar (r, z) values. This is true despite the HFB calculations being unconstrained, and is a consequence of the EDF parameters all predicting a rather similar ground state deformation. Contrast with e.g. taking samples all across a PES, where the nuclear shape differs considerably (recall Fig. 3.1). Similarly, the first basis function is essentially the same as the sample field, except for an overall factor. Next, notice that the second basis function is peaked around the nuclear surface. This suggests at first that surface properties vary somewhat with changing EDF parameters. However, the expansion coefficient for this function is rather small, as evidenced by Fig. 6.3. So, the form of the basis function is not due to surface properties of the nucleus. Rather, it is a consequence of generating functions orthogonal to Fig. 6.4(b) (consider, for instance, the spectrum of the harmonic oscillator), and is in fact related to the Courant nodal domain 78 theorem mentioned above. The basis functions for other fields follow the same pattern. 6.3.3 The Payoff Figure 6.5 plots the runtime per iteration against the error in the HFB energy as the two thresholds εei and εRBM are varied. Also shown are the runtime per iteration of HFBTHO and a Python axial HFB solver. The overall runtime is determined by (i) the initial guess and (ii) the convergence tolerance of the iterative procedure. It is not obvious how to obtain equivalent initial guesses between the true solvers and the emulators. Also, the convergence tolerance should in principle be chosen based on the desired precision of the solver. Thus, the overall performance of the emulator may be misrepresented if the total runtime is displayed in place of the runtime per iteration. Figure 6.5: The runtime per iteration for the 25-oscillator-shell unconstrained ground state of 254Fm. EI is empirical interpolation, RBM uses a tailored basis for the quasiparticle wave functions, and MaxVol evaluates the fields on a minimal coordinate grid. First, notice that as the thresholds are increased, the overall error increases, while the 79 10−1100101102RuntimeperIteration(s)10−510−410−310−210−1100EHFBError(MeV)HFBTHOPythonHFBSolverEIEI+RBMEI+MaxVolEI+MaxVol+RBM runtime per iteration remains constant. This is because the runtime per iteration is dom- inated by the HFB matrix diagonalization and the wavefunction reconstruction, neither of which depend on the number of EI basis functions. Extremely precise results are obtain- able, to eV precision - far below both the calibration error of the EDF [40] and the basis truncation error of the solver [199]. Alternatively, keV precision may be obtained with far fewer high-fidelity samples, provided they are distributed adequately across the parameter range (e.g. using a Latin hypercube sampling [221, 225]). This accuracy is attained with an order-of-magnitude speedup over HFBTHO. Next, note the lack of speedup obtained using a tailored basis for the single-quasiparticle wave functions, denoted in Fig. 6.5 with the label “RBM”. This confirms the intuition mentioned in Sec. 6.3.1. Actually, the performance is often worse: for εRBM ≳ 10− 3, high- lying states are not adequately represented in the reduced basis. This leads to errors in EHFB of 10s of MeV. Thus, this approach for reducing the diagonalization cost is insufficient in the best case, and harmful in the worst case. Alternatively, recent work has successfully used machine learning approaches to obtain orbital-free EDFs [226]. Such an approach avoids diagonalizing the HFB matrix entirely, and is therefore an interesting avenue for future research. Finally, it is worth emphasizing an additional utility of this approach: because the den- sities and wave functions are computed exactly, downstream quantities are straightforward to obtain. This includes quantities important for fission calculations, such as the collective inertia (Sec. 3.1.1) and the localization functional (Sec. 4.1.3). A representative example of the latter is shown in Fig. 6.6, for the highest-precision emulator. As can be seen, there is no visible difference between the quantity using the real HFB solver and the emulated result. 80 Figure 6.6: The proton localization functional for the 25-oscillator-shell unconstrained ground state of 254Fm. The left (right) panel uses the densities from the HFB calcula- tion (the emulator). The oscillator widths are bz = 2.242, b = 2.464, as determined in HFBTHO [36]. ⊥ 6.4 Conclusion In this chapter, I have presented the general philosophy behind reduced order models. I have demonstrated how projection-based emulators may be modified for use with self-consistent calculations with non-affine parameter dependence. And, I have demonstrated that this approach leads to order-of-magnitude speedup in axial HFB calculations, with essentially exact accuracy. Explicitly, this approach addresses the concerns raised at the end of Ch. 5. First, the nucleonic densities and currents are computed, avoiding excess HFB runs once the exit point is located. Second, the training data consists of fewer than 100 HFB runs, which is promising despite the more well-constrained parameter region (the EDF parameters, as opposed to (A, Z, Q20, Q30)). And third, the emulator error is arbitrarily precise. Additionally, the 81 −4−2024r/b⊥−4−2024z/bzRealEmulated0.080.160.240.320.400.480.560.640.72 emulator error is straightforward to check, whereas with the NN approach it is impossible. And, the emulator is easy to update as additional training data is generated, as opposed to retraining a NN. There are a number of future directions branching from this topic. As far as emula- tor improvements go, reduced coordinate-space solvers are worth exploring, as mentioned in Sec. 6.2.2. Also worth considering are approximate orbital-free EDFs as discussed in Ref. [226]. As far as applications go, one option is to include (A, Z) and shape parame- ters as additional model parameters, with an eye towards similar large-scale calculations as described in Ch. 5. Another is applications to other EDFs, as were described in Sec. 2.2. And, a third is to integrate this approach into existing solvers, so as to enable quicker error propagation in the tunneling calculations described in Ch. 3. 82 Chapter 7. Conclusions and Outlook In this thesis, we dealt with the implementation and usage of the nudged elastic band method for fission tunneling pathway calculations, and the development of emulators for use in large- scale fission studies. In Ch. 3, I described the tunneling process in spontaneous fission, beginning from the reduction to collective coordinates, and ending with the practical calculation of tunneling pathways. I discussed my development and implementation of the NEB algorithm, which agrees to the percent-level with previously-used algorithms. It also scales well with the number of collective coordinates, allowing for efficient use with three or more coordinates. In Ch. 4, I applied NEB to the Fermium chain to study multimodal fission. In addition to qualitative agreement with previously-identified fission modes, we find that the primary fission fragment yields are quite sensitive to the EDF of choice when multiple competing modes are present, with the UNEDF1HFB EDF agreeing the best with experimental data. The results are only weakly sensitive to the number of collective coordinates used. A number of future directions for this approach are possible. On the development front, the optimal NEB hyperparameters have not been determined. It is also interesting to con- sider fluctuations around the tunneling path, and the impact of additional collective coor- dinates. On the physics front, it is interesting to apply NEB to tunneling in microscopic- macroscopic approaches to fission [227]. I am also also applying this approach to induced fission, where there is evidence suggesting competition between fission modes as the excita- tion energy of the fissioning nucleus is increased [183, 184]. In Ch. 5, I demonstrated that neural networks are able to emulate the PES to within ≲ 1 MeV across the r process region of the nuclear chart. This accuracy translates to reasonable accuracies in the half-life and the primary fission fragment yields. And, in Ch. 6, I have 83 demonstrated that simple reduced order models may be applied to axial HFB calculations to achieve an order-of-magnitude speedup with essentially exact accuracy. The latter approach also requires fewer than 100 HFB calculations as training data, and no additional HFB calculations are necessary to compute other quantities, such as the collective inertia tensor. There are many future directions to explore. In regards to emulator development, reduced coordinate-space solvers are worth exploring, as are approximate orbital-free EDFs [226]. With regards to applications, one option is to include (A, Z) and shape parameters as addi- tional model parameters, with an eye towards large-scale calculations as described in Ch. 5. Another is applications to other EDFs, as described in Sec. 2.2. And, a third is to inte- grate this approach into existing solvers, so as to enable quicker error propagation in PES calculations. Finally, there are a number of recent exciting developments in both fission experiment and theory. On the experimental side, the angular momentum of the primary fragments has recently been measured [228], as has neutrino-induced fission [229], both of which are interesting to explore from a theoretical point of view. And, on the theory side, recent developments on the description of fragment yields [230] and the tunneling process [231] are also exciting. It is indeed an interesting time to study nuclear fission theory. 84 BIBLIOGRAPHY [1] O. Hahn and F. Strassmann. “ ¨Uber den Nachweis und das Verhalten der bei der Bestrahlung des Urans mittels Neutronen entstehenden Erdalkalimetalle”. In: Natur- wissenschaften 27.1 (1939), pp. 11–15. doi: 10.1007/bf01488241. [2] L. Meitner and O. R. Frisch. “Disintegration of uranium by neutrons: a new type of nuclear reaction”. In: Nature 143.3615 (1939), pp. 239–240. [3] U. D. of Energy (USDOE). A New Era of Discovery: The 2023 Long Range Plan for Nuclear Science. Tech. rep. US Department of Energy (USDOE), Washington, DC (United States). Office of Science, Oct. 2023. doi: 10.2172/2280968. [4] M. Bender et al. “Future of nuclear fission theory”. In: J. Phys. G Nucl. Part. Phys. 47.11 (Oct. 2020), p. 113002. doi: 10.1088/1361-6471/abab4f. [5] J. J. Cowan et al. “Origin of the heaviest elements: The rapid neutron-capture pro- cess”. In: Rev. Mod. Phys. 93.1 (1 Feb. 2021), p. 015002. doi: 10.1103/RevModPhys. 93.015002. [6] R. Surman and M. Mumpower. “Masses and lifetimes for r-process nucleosynthesis: FRIB outlook”. In: EPJ Web Conf. 178 (2018), p. 04002. doi: 10 . 1051 / epjconf / 201817804002. [7] L. Neufcourt et al. “Quantified limits of the nuclear landscape”. In: Phys. Rev. C 101 (4 Apr. 2020), p. 044307. doi: 10.1103/PhysRevC.101.044307. [8] N. Vassh et al. “Using excitation-energy dependent fission yields to identify key fission- ing nuclei in r-process nucleosynthesis”. In: J. Phys. G 46.6 (June 2019), p. 065202. doi: 10.1088/1361-6471/ab0bea. [9] S. A. Giuliani et al. “Fission and the r-process nucleosynthesis of translead nuclei in neutron star mergers”. In: Phys. Rev. C 102 (4 Oct. 2020), p. 045804. doi: 10.1103/ PhysRevC.102.045804. [10] P. Giuliani et al. “Bayes goes fast: Uncertainty quantification for a covariant energy density functional emulated by the reduced basis method”. In: Front. Phys. 10 (2023). doi: 10.3389/fphy.2022.1054524. [11] J. D. McDonnell et al. “Uncertainty Quantification for Nuclear Density Functional Theory and Information Content of New Measurements”. In: Phys. Rev. Lett. 114 (12 Mar. 2015), p. 122501. doi: 10.1103/PhysRevLett.114.122501. 85 [12] A. Boehnlein et al. “Colloquium: Machine learning in nuclear physics”. In: Rev. Mod. Phys. 94 (3 Sept. 2022), p. 031003. doi: 10.1103/RevModPhys.94.031003. [13] E. Flynn et al. “Nudged elastic band approach to nuclear fission pathways”. In: Phys. Rev. C 105.5 (5 May 2022), p. 054302. doi: 10.1103/PhysRevC.105.054302. [14] D. Lay et al. “Multimodal fission from self-consistent calculations”. In: Phys. Rev. C 109.4 (2024), p. 044306. doi: 10.1103/PhysRevC.109.044306. [15] D. Lay et al. “Neural network emulation of spontaneous fission”. In: Phys. Rev. C 109.4 (2024), p. 044305. doi: 10.1103/PhysRevC.109.044305. [16] R. G. Parr and W. Yang. Density-functional theory of atoms and molecules. Oxford University Press, 1989. [17] D. S. Sholl and J. A. Steckel. Density functional theory: a practical introduction. John Wiley & Sons, 2022. [18] P. Hohenberg and W. Kohn. “Inhomogeneous Electron Gas”. In: Phys. Rev. 136 (3b Nov. 1964), B864–b871. doi: 10.1103/PhysRev.136.B864. [19] W. Kohn and L. J. Sham. “Self-Consistent Equations Including Exchange and Cor- relation Effects”. In: Phys. Rev. 140 (4a Nov. 1965), A1133–a1138. doi: 10 . 1103 / PhysRev.140.A1133. [20] N. Schunck, ed. Energy Density Functional Methods for Atomic Nuclei. 2053-2563. IOP Publishing, 2019. doi: 10.1088/2053-2563/aae0ed. [21] N. Schunck and L. M. Robledo. “Microscopic theory of nuclear fission: a review”. In: Rep. Prog. Phys. 79.11 (Oct. 2016), p. 116301. doi: 10.1088/0034-4885/79/11/116301. [22] P. Ring and P. Schuck. The nuclear many-body problem. New York: Springer-Verlag, 1980. [23] J.-P. Blaizot and G. Ripka. Quantum theory of finite systems. eng. Cambridge, Mass: MIT Press, 1986. [24] J. Dobaczewski, H. Flocard, and J. Treiner. “Hartree-Fock-Bogolyubov description of nuclei near the neutron-drip line”. In: Nuclear Physics A 422.1 (1984), pp. 103–139. doi: 10.1016/0375-9474(84)90433-0. [25] J. Dobaczewski et al. “Mean-field description of ground-state properties of drip-line nuclei: Pairing and continuum effects”. In: Physical Review C 53.6 (1996), p. 2809. 86 [26] J. Dobaczewski, W. Nazarewicz, and M. V. Stoitsov. “Nuclear ground-state properties from mean-field calculations”. In: Eur. Phys. J. A 15.1 (2002), pp. 21–26. doi: 10. 1140/epja/i2001-10218-8. [27] L. M. Robledo, T. R. Rodr´ıguez, and R. R. Rodr´ıguez-Guzm´an. “Mean field and beyond description of nuclear structure with the Gogny force: a review”. In: J. Phys. G: Nucl. Part. Phys. 46.1 (Jan. 2019), p. 013001. doi: 10.1088/1361-6471/aadebd. [28] M. Bender, P.-H. Heenen, and P.-G. Reinhard. “Self-consistent mean-field models for nuclear structure”. In: Rev. Mod. Phys. 75.1 (1 Jan. 2003), pp. 121–180. doi: 10.1103/RevModPhys.75.121. [29] B. A. Brown and B. Wildenthal. “Status of the nuclear shell model”. In: Annual Review of Nuclear and Particle Sciences 38.1 (1988), pp. 29–66. [30] E. Caurier et al. “The shell model as a unified view of nuclear structure”. In: Rev. Mod. Phys. 77 (2 June 2005), pp. 427–488. doi: 10.1103/RevModPhys.77.427. [31] A. N. Bohr and B. R. Mottelson. Nuclear Structure (in 2 volumes). World Scientific Publishing Company, 1998. [32] D. J. Dean and M. Hjorth-Jensen. “Pairing in nuclear systems: from neutron stars to finite nuclei”. In: Rev. Mod. Phys. 75 (2 Apr. 2003), pp. 607–656. doi: 10.1103/ RevModPhys.75.607. [33] R. N. Perez et al. “Axially deformed solution of the Skyrme–Hartree–Fock–Bogolyubov equations using the transformed harmonic oscillator basis (III) hfbtho (v3.00): A new version of the program”. In: Computer Physics Communications 220 (2017), pp. 363– 375. doi: 10.1016/j.cpc.2017.06.022. [34] N. Schunck et al. “Solution of the Skyrme–Hartree–Fock–Bogolyubov equations in the Cartesian deformed harmonic-oscillator basis.: (VII) hfodd (v2.49t): A new version of the program”. In: Computer Physics Communications 183.1 (2012), pp. 166–192. doi: 10.1016/j.cpc.2011.08.013. [35] L. M. Robledo and G. F. Bertsch. “Application of the gradient method to Hartree- Fock-Bogoliubov theory”. In: Phys. Rev. C 84.1 (July 2011), p. 014312. doi: 10.1103/ PhysRevC.84.014312. [36] M. Stoitsov et al. “Axially deformed solution of the Skyrme–Hartree–Fock–Bogolyubov equations using the transformed harmonic oscillator basis. The program HFBTHO (v1.66p)”. In: Computer Physics Communications 167.1 (2005), pp. 43–63. doi: 10. 1016/j.cpc.2005.01.001. 87 [37] D. Vautherin and D. M. Brink. “Hartree-Fock Calculations with Skyrme’s Interaction. I. Spherical Nuclei”. In: Phys. Rev. C 5 (3 Mar. 1972), pp. 626–647. doi: 10.1103/ PhysRevC.5.626. [38] D. Vautherin. “Hartree-Fock calculations with Skyrme’s interaction. II. Axially de- formed nuclei”. In: Physical Review C 7.1 (1973), p. 296. [39] M. Kortelainen et al. “Nuclear energy density optimization”. In: Phys. Rev. C 82 (2 Aug. 2010), p. 024313. doi: 10.1103/PhysRevC.82.024313. [40] M. Kortelainen et al. “Nuclear energy density optimization: Large deformations”. In: Phys. Rev. C 85 (2 Feb. 2012), p. 024304. doi: 10.1103/PhysRevC.85.024304. [41] M. Kortelainen et al. “Nuclear energy density optimization: Shell structure”. In: Phys. Rev. C 89 (5 May 2014), p. 054314. doi: 10.1103/PhysRevC.89.054314. [42] J. C. Slater. “A Simplification of the Hartree-Fock Method”. In: Phys. Rev. 81 (3 Feb. 1951), pp. 385–390. doi: 10.1103/PhysRev.81.385. [43] N. Schunck et al. “Solution of the Skyrme-Hartree–Fock–Bogolyubovequations in the Cartesian deformed harmonic-oscillator basis. (VIII) hfodd (v2.73y): A new version of the program”. In: Computer Physics Communications 216 (2017), pp. 145–174. doi: 10.1016/j.cpc.2017.03.007. [44] P. Marevi´c et al. “Axially-deformed solution of the Skyrme-Hartree-Fock-Bogoliubov equations using the transformed harmonic oscillator basis (IV) HFBTHO (v4. 0): A new version of the program”. In: Computer physics communications 276 (2022), p. 108367. [45] T. Naito et al. “Coulomb exchange functional with generalized gradient approxima- tion for self-consistent Skyrme Hartree-Fock calculations”. In: Phys. Rev. C 99 (2 Feb. 2019), p. 024309. doi: 10.1103/PhysRevC.99.024309. [46] C. Titin-Schnaider and P. Quentin. “Coulomb exchange contribution in nuclear Hartree- Fock calculations”. In: Physics Letters B 49.5 (1974), pp. 397–400. doi: 10.1016/0370- 2693(74)90617-0. [47] H.-Q. Gu et al. “Slater approximation for Coulomb exchange effects in nuclear co- variant density functional theory”. In: Phys. Rev. C 87 (4 Apr. 2013), p. 041301. doi: 10.1103/PhysRevC.87.041301. [48] T. Skyrme. “The effective nuclear potential”. In: Nuclear Physics 9.4 (1958), pp. 615– 634. doi: 10.1016/0029-5582(58)90345-6. 88 [49] D. Gogny. “Simple separable expansions for calculating matrix elements of two-body local interactions with harmonic oscillator functions”. In: Nuclear Physics A 237.3 (1975), pp. 399–418. doi: 10.1016/0375-9474(75)90407-8. [50] J. Decharg´e and D. Gogny. “Hartree-Fock-Bogolyubov calculations with the D1 ef- fective interaction on spherical nuclei”. In: Phys. Rev. C 21 (4 Apr. 1980), pp. 1568– 1593. doi: 10.1103/PhysRevC.21.1568. [51] S. Fayans. “Towards a universal nuclear density functional”. In: Journal of Experi- mental and Theoretical Physics Letters 68 (1998), pp. 169–174. [52] P.-G. Reinhard and W. Nazarewicz. “Toward a global description of nuclear charge radii: Exploring the Fayans energy density functional”. In: Phys. Rev. C 95 (6 June 2017), p. 064328. doi: 10.1103/PhysRevC.95.064328. [53] M. Baldo et al. “Accurate nuclear masses from a three parameter Kohn-Sham DFT approach (BCPM)”. In: AIP Conference Proceedings 1491.1 (Oct. 2012), pp. 218–221. doi: 10.1063/1.4764242. [54] A. Bulgac et al. “Minimal nuclear energy density functional”. In: Phys. Rev. C 97 (4 Apr. 2018), p. 044313. doi: 10.1103/PhysRevC.97.044313. [55] J. Drut, R. Furnstahl, and L. Platter. “Toward ab initio density functional theory for nuclei”. In: Progress in Particle and Nuclear Physics 64.1 (2010), pp. 120–168. doi: 10.1016/j.ppnp.2009.09.001. [56] T. Duguet et al. “Rooting the EDF method into the ab initio framework: PGCM- PT formalism based on MR-IMSRG pre-processed Hamiltonians”. In: The European Physical Journal A 59.1 (2023), p. 13. [57] F. Marino et al. “Nuclear energy density functionals grounded in ab initio calcula- tions”. In: Phys. Rev. C 104 (2 Aug. 2021), p. 024315. doi: 10.1103/PhysRevC.104. 024315. [58] L. Zurek et al. “Optimized nuclear energy density functionals including long-range pion contributions”. In: Phys. Rev. C 109 (1 Jan. 2024), p. 014319. doi: 10.1103/ PhysRevC.109.014319. [59] J. Dobaczewski. “Ab initio derivation of model energy density functionals”. In: Jour- nal of Physics G: Nuclear and Particle Physics 43.4 (Feb. 2016), 04lt01. doi: 10 . 1088/0954-3899/43/4/04lt01. [60] P. Ring. “Relativistic mean field theory in finite nuclei”. In: Progress in Particle and Nuclear Physics 37 (1996), pp. 193–263. doi: https://doi.org/10.1016/0146- 6410(96)00054-3. 89 [61] T. Nikˇsi´c, D. Vretenar, and P. Ring. “Relativistic nuclear energy density functionals: Mean-field and beyond”. In: Progress in Particle and Nuclear Physics 66.3 (2011), pp. 519–548. doi: https://doi.org/10.1016/j.ppnp.2011.01.055. [62] E. Perlinska et al. “Local density approximation for proton-neutron pairing corre- lations: Formalism”. In: Phys. Rev. C 69 (1 Jan. 2004), p. 014316. doi: 10.1103/ PhysRevC.69.014316. [63] H. K¨ohler. “Skyrme force and the mass formula”. In: Nuclear Physics A 258.2 (1976), pp. 301–316. [64] J. Bartel et al. “Towards a better parametrisation of Skyrme-like effective forces: A critical study of the SkM force”. In: Nucl. Phys. A 386.1 (1982), pp. 79–100. doi: 10.1016/0375-9474(82)90403-1. [65] P. J. Borycki et al. “Pairing renormalization and regularization within the local den- sity approximation”. In: Phys. Rev. C 73 (4 Apr. 2006), p. 044319. doi: 10.1103/ PhysRevC.73.044319. [66] J. Stone and P.-G. Reinhard. “The Skyrme interaction in finite nuclei and nuclear matter”. In: Progress in Particle and Nuclear Physics 58.2 (2007), pp. 587–657. doi: 10.1016/j.ppnp.2006.07.001. [67] H. Krivine, J. Treiner, and O. Bohigas. “Derivation of a fluid-dynamical lagrangian and electric giant resonances”. In: Nuclear Physics A 336.2 (1980), pp. 155–184. doi: 10.1016/0375-9474(80)90618-1. [68] D. Gogny. Self-consistent pairing calculations. Ed. by G. Ripka and M. Porneuf. 1975. [69] J. Decharge, M. Girod, and D. Gogny. “Self consistent calculations and quadrupole moments of even Sm isotopes”. In: Physics Letters B 55.4 (1975), pp. 361–364. doi: 10.1016/0370-2693(75)90359-7. [70] J. Berger, M. Girod, and D. Gogny. “Microscopic analysis of collective dynamics in low energy fission”. In: Nuclear Physics A 428 (1984), pp. 23–36. doi: 10.1016/0375- 9474(84)90240-9. [71] J. Berger, M. Girod, and D. Gogny. “Time-dependent quantum collective dynam- ics applied to nuclear fission”. In: Computer Physics Communications 63.1 (1991), pp. 365–374. doi: https://doi.org/10.1016/0010-4655(91)90263-K. [72] J. Sadhukhan et al. “Pairing-induced speedup of nuclear spontaneous fission”. In: Phys. Rev. C 90 (6 Dec. 2014), 061304(r). doi: 10.1103/PhysRevC.90.061304. 90 [73] M. Bender, T. Duguet, and D. Lacroix. “Particle-number restoration within the en- ergy density functional formalism”. In: Phys. Rev. C 79 (4 Apr. 2009), p. 044319. doi: 10.1103/PhysRevC.79.044319. [74] L. M. Robledo. “Particle Number Restoration: Its Implementation And Impact In Nuclear Structure Calculations”. In: International Journal of Modern Physics E 16.02 (2007), pp. 337–351. doi: 10.1142/s0218301307005776. [75] G. Hupin, D. Lacroix, and M. Bender. “Formulation of functional theory for pairing with particle number restoration”. In: Phys. Rev. C 84 (1 July 2011), p. 014309. doi: 10.1103/PhysRevC.84.014309. [76] M. V. Stoitsov et al. “Variation after particle-number projection for the Hartree-Fock- Bogoliubov method with the Skyrme energy density functional”. In: Phys. Rev. C 76 (1 July 2007), p. 014308. doi: 10.1103/PhysRevC.76.014308. [77] J. A. Sheikh et al. “Symmetry restoration in mean-field approaches”. In: Journal of Physics G: Nuclear and Particle Physics 48.12 (Nov. 2021), p. 123001. doi: 10.1088/ 1361-6471/ac288a. [78] H. Hergert. “A Guided Tour of ab initio Nuclear Many-Body Theory”. In: Frontiers in Physics 8 (2020). doi: 10.3389/fphy.2020.00379. [79] H. J. Lipkin. “Collective motion in many-particle systems: Part 1. The violation of conservation laws”. In: Annals of physics 9.2 (1960), pp. 272–291. [80] Y. Nogami. “Improved Superconductivity Approximation for the Pairing Interaction in Nuclei”. In: Phys. Rev. 134 (2b Apr. 1964), B313–b321. doi: 10.1103/PhysRev. 134.B313. [81] M. V. Stoitsov et al. “Systematic study of deformed nuclei at the drip lines and beyond”. In: Phys. Rev. C 68 (5 Nov. 2003), p. 054312. doi: 10.1103/PhysRevC.68. 054312. [82] P.-G. Reinhard et al. “Lipkin-Nogami pairing scheme in self-consistent nuclear struc- ture calculations”. In: Phys. Rev. C 53 (6 June 1996), pp. 2776–2785. doi: 10.1103/ PhysRevC.53.2776. [83] N. Schunck et al. “Error analysis in nuclear density functional theory”. In: J. Phys. G 42.3 (Feb. 2015), p. 034024. doi: 10.1088/0954-3899/42/3/034024. [84] J. Dobaczewski et al. HFODD (v2.40h) User’s Guide. 2009. url: https://arxiv.org/ abs/0909.3626. 91 [85] R. W. Hasse and W. D. Myers. Geometrical relationships of macroscopic nuclear physics. Springer Science & Business Media, 1988. [86] S. Larsson, I. Ragnarsson, and S. Nilsson. “Fission barriers and the inclusion of axial asymmetry”. In: Physics Letters B 38.5 (1972), pp. 269–273. doi: 10 . 1016 / 0370 - 2693(72)90243-2. [87] A. Staszczak, J. Dobaczewski, and W. Nazarewicz. “Skyrme–hartree–fock Calcu- lations Of Fission Barriers Of The Heaviest And Superheavy Nuclei”. In: Inter- national Journal of Modern Physics E 14.03 (2005), pp. 395–402. doi: 10 . 1142 / s0218301305003181. [88] A. Staszczak, A. Baran, and W. Nazarewicz. “Breaking Of Axial And Reflection Symmetries In Spontaneous Fission Of Fermium Isotopes”. In: Int. J. Mod. Phys. E 20.02 (2011), pp. 552–556. doi: 10.1142/s0218301311017995. [89] Y. Cao et al. “Landscape of pear-shaped even-even nuclei”. In: Phys. Rev. C 102 (2 Aug. 2020), p. 024311. doi: 10.1103/PhysRevC.102.024311. [90] N. Schunck and D. Regnier. “Theory of nuclear fission”. In: Prog. Part. Nucl. Phys. 125 (July 2022), p. 103963. doi: 10.1016/j.ppnp.2022.103963. [91] A. Staszczak et al. “Microscopic description of complex nuclear decay: Multimodal fission”. In: Phys. Rev. C 80 (1 July 2009), p. 014309. doi: 10.1103/PhysRevC.80. 014309. [92] B. Giraud, J. Le Tourneux, and S. Wong. “Constrained Hartree-Fock calculations for a quadrupole generator coordinate”. In: Physics Letters B 32.1 (1970), pp. 23–25. doi: 10.1016/0370-2693(70)90325-4. [93] W. Younes and D. Gogny. “Microscopic calculation of 240Pu scission with a finite- range effective force”. In: Phys. Rev. C 80 (5 Nov. 2009), p. 054313. doi: 10.1103/ PhysRevC.80.054313. [94] M. Stoitsov et al. “Axially deformed solution of the Skyrme-Hartree–Fock–Bogoliubov equations using the transformed harmonic oscillator basis (II) hfbtho v2. 00d: A new version of the program”. In: Computer Physics Communications 184.6 (2013), pp. 1592–1604. [95] A. Baran et al. “Broyden’s method in nuclear structure calculations”. In: Physical Review C–Nuclear Physics 78.1 (2008), p. 014318. [96] D. M. Brink and R. A. Broglia. Nuclear superfluidity: pairing in finite systems. Cam- bridge University Press, 2005. 92 [97] A. Staszczak et al. “Coupling of the pairing vibrations with the fission mode”. In: Physics Letters B 161.4 (1985), pp. 227–230. doi: 10.1016/0370-2693(85)90750-6. [98] N. L. Vaquero, T. R. Rodr´ıguez, and J. L. Egido. “On the impact of large amplitude pairing fluctuations on nuclear spectra”. In: Phys. Lett. B 704 (2011), pp. 520–526. doi: 10.1016/j.physletb.2011.09.073. [99] N. L. Vaquero, J. L. Egido, and T. R. Rodr´ıguez. “Large-amplitude pairing fluctua- tions in atomic nuclei”. In: Phys. Rev. C 88 (2013), p. 064311. doi: 10.1103/physrevc. 88.064311. [100] E. Runge and E. K. U. Gross. “Density-Functional Theory for Time-Dependent Systems”. In: Phys. Rev. Lett. 52 (12 Mar. 1984), pp. 997–1000. doi: 10 . 1103 / PhysRevLett.52.997. [101] M. Baranger and M. V´en´eroni. “An adiabatic time-dependent Hartree-Fock theory of collective motion in finite systems”. In: Annals of Physics 114.1 (1978), pp. 123–200. doi: 10.1016/0003-4916(78)90265-8. [102] A. Baran et al. “Quadrupole collective inertia in nuclear fission: Cranking approxi- mation”. In: Phys. Rev. C 84.5 (5 Nov. 2011), p. 054321. doi: 10.1103/PhysRevC. 84.054321. [103] J. J. Sakurai and J. Napolitano. Modern quantum mechanics; 2nd ed. San Francisco, CA: Addison-Wesley, 2011. url: https://cds.cern.ch/record/1341875. [104] S. A. Giuliani and L. M. Robledo. “Non-perturbative collective inertias for fission: A comparative study”. In: Phys. Lett. B 787 (Dec. 2018), pp. 134–140. doi: 10.1016/j. physletb.2018.10.045. [105] T. Nakatsukasa, T. Inakura, and K. Yabana. “Finite amplitude method for the so- lution of the random-phase approximation”. In: Phys. Rev. C 76 (2 Aug. 2007), p. 024318. doi: 10.1103/PhysRevC.76.024318. [106] N. Hinohara. “Collective inertia of the Nambu-Goldstone mode from linear response theory”. In: Phys. Rev. C 92 (3 Sept. 2015), p. 034321. doi: 10.1103/PhysRevC.92. 034321. [107] M. J. Giannoni and P. Quentin. “Mass parameters in the adiabatic time-dependent Hartree-Fock approximation. II. Results for the isoscalar quadrupole mode”. In: Phys. Rev. C 21 (5 May 1980), pp. 2076–2093. doi: 10.1103/PhysRevC.21.2076. [108] E. Yuldashbaeva et al. “Mass parameters for large amplitude collective motion: A perturbative microscopic approach”. In: Phys. Lett. B 461.1-2 (Aug. 1999), pp. 1–8. doi: 10.1016/s0370-2693(99)00836-9. 93 [109] J. Sadhukhan et al. “Spontaneous fission lifetimes from the minimization of self- consistent collective action”. In: Phys. Rev. C 88 (6 Dec. 2013), p. 064314. doi: 10.1103/PhysRevC.88.064314. [110] A. Staszczak, S. Pi(cid:32)lat, and K. Pomorski. “Influence of the pairing vibrations on spontaneous fission probability”. In: Nucl. Phys. A 504.3 (1989), pp. 589–604. doi: 10.1016/0375-9474(89)90559-9. [111] T. Nakatsukasa et al. “Time-dependent density-functional description of nuclear dy- namics”. In: Rev. Mod. Phys. 88 (4 Nov. 2016), p. 045004. doi: 10.1103/RevModPhys. 88.045004. [112] Z. H. Huang et al. “Wentzel-Kramers-Brillouin method in multidimensional tunnel- ing”. In: Phys. Rev. A 41 (1 Jan. 1990), pp. 32–41. doi: 10.1103/PhysRevA.41.32. [113] A. Schmid. “Quasiclassical wave function in multidimensional quantum decay prob- lems”. In: Annals of Physics 170.2 (1986), pp. 333–369. doi: 10.1016/0003-4916(86) 90096-5. [114] K. Takatsuka, H. Ushiyama, and A. Inoue-Ushiyama. “Tunneling paths in multi- dimensional semiclassical dynamics”. In: Physics Reports 322.5 (1999), pp. 347–417. doi: 10.1016/s0370-1573(99)00036-8. [115] R. P. Feynman, A. R. Hibbs, and D. F. Styer. Quantum mechanics and path integrals. Courier Corporation, 2010. [116] S. Coleman. “Quantum tunneling and negative eigenvalues”. In: Nuclear Physics B 298.1 (1988), pp. 178–186. doi: 10.1016/0550-3213(88)90308-2. [117] C. M. Bender and S. A. Orszag. Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory. Springer Science & Business Media, 2013. [118] Y. Zhu and J. C. Pei. “Thermal fission rates with temperature dependent fission barriers”. In: Phys. Rev. C 94 (2 Aug. 2016), p. 024329. doi: 10.1103/PhysRevC.94. 024329. [119] N. G. Kelkar and H. M. Casta˜neda. “Critical view of WKB decay widths”. In: Phys. Rev. C 76 (6 Dec. 2007), p. 064605. doi: 10.1103/PhysRevC.76.064605. [120] T. Ledergerber and H.-C. Pauli. “On the dynamics of fission: The role of reflection asymmetry in the nuclear shape”. In: Nucl. Phys. A 207.1 (1973), pp. 1–32. doi: 10.1016/0375-9474(73)90022-5. 94 [121] A. Baran. “Some dynamical aspects of the fission process”. In: Phys. Lett. B 76.1 (1978), pp. 8–10. doi: 10.1016/0370-2693(78)90085-0. [122] K. Okada et al. “Cassini-oval description of the multidimensional potential energy surface for 236U: Role of octupole deformation and calculation of the most probable fission path”. In: Phys. Rev. C 107 (3 Mar. 2023), p. 034608. doi: 10.1103/PhysRevC. 107.034608. [123] A. Iwamoto. “Multi-dimensional tunneling and nuclear fission process”. In: Z. Phys. A 349.3 (1994), p. 265. doi: 10.1007/bf01288972. [124] G. Scamps, C. Simenel, and D. Lacroix. “Superfluid dynamics of 258Fm fission”. In: Phys. Rev. C 92 (1 July 2015), 011602(r). doi: 10.1103/PhysRevC.92.011602. [125] J. Skalski. “Nuclear fission with mean-field instantons”. In: Phys. Rev. C 77 (6 June 2008), p. 064610. doi: 10.1103/PhysRevC.77.064610. [126] M. Baranger, M. Strayer, and J.-S. Wu. “Nuclear collective motion with full nonlinear- ity”. In: Phys. Rev. C 67 (1 Jan. 2003), p. 014318. doi: 10.1103/PhysRevC.67.014318. [127] J. Zhao et al. “Multidimensionally constrained relativistic Hartree-Bogoliubov study of spontaneous nuclear fission”. In: Phys. Rev. C 92 (6 Dec. 2015), p. 064315. doi: 10.1103/PhysRevC.92.064315. [128] F. Mercier et al. “Microscopic Description of 2α Decay in 212Po and 224Ra Isotopes”. In: Phys. Rev. Lett. 127 (1 July 2021), p. 012501. doi: 10.1103/PhysRevLett.127. 012501. [129] A. Baran et al. “A dynamic analysis of spontaneous-fission half-lives”. In: Nuclear Physics A 361.1 (1981), pp. 83–101. doi: 10.1016/0375-9474(81)90471-1. [130] E. W. Dijkstra. “A note on two problems in connexion with graphs”. In: Numer. Math. 1.1 (1959), pp. 269–271. doi: 10.1007/bf01386390. [131] G. Mills and H. J´onsson. “Quantum and thermal effects in H2 dissociative adsorption: Evaluation of free energy barriers in multidimensional quantum systems”. In: Phys. Rev. Lett. 72 (7 Feb. 1994), pp. 1124–1127. doi: 10.1103/PhysRevLett.72.1124. [132] G. Mills, H. J´onsson, and G. K. Schenter. “Reversible work transition state the- ory: application to dissociative adsorption of hydrogen”. In: Surf. Sci. 324.2 (1995), pp. 305–337. doi: 10.1016/0039-6028(94)00731-4. [133] H. J´onsson, G. Mills, and K. W. Jacobsen. “Nudged elastic band method for find- ing minimum energy paths of transitions”. In: Classical and Quantum Dynamics in Condensed Phase Simulations, pp. 385–404. doi: 10.1142/9789812839664 0016. 95 [134] G. Henkelman, B. P. Uberuaga, and H. J´onsson. “A climbing image nudged elastic band method for finding saddle points and minimum energy paths”. In: J. Chem. Phys. 113.22 (2000), pp. 9901–9904. doi: 10.1063/1.1329672. [135] G. Henkelman and H. J´onsson. “Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points”. In: J. Chem. Phys. 113.22 (2000), pp. 9978–9985. doi: 10.1063/1.1323224. [136] V. ´Asgeirsson, A. Arnaldsson, and H. J´onsson. “Efficient evaluation of atom tunneling combined with electronic structure calculations”. In: J. Chem. Phys. 148.10 (2018), p. 102334. doi: 10.1063/1.5007180. [137] B. J. Berne, G. Ciccotti, and D. F. Coker. Classical and quantum dynamics in con- densed phase simulations: Proceedings of the International School of Physics. World Scientific, 1998. [138] L. Verlet. “Computer ”Experiments” on Classical Fluids. I. Thermodynamical Prop- erties of Lennard-Jones Molecules”. In: Phys. Rev. 159 (1 July 1967), pp. 98–103. doi: 10.1103/PhysRev.159.98. [139] E. Bitzek et al. “Structural Relaxation Made Simple”. In: Phys. Rev. Lett. 97 (17 Oct. 2006), p. 170201. doi: 10.1103/PhysRevLett.97.170201. [140] J. Gu´enol´e et al. “Assessment and optimization of the fast inertial relaxation engine (fire) for energy minimization in atomistic simulations and its implementation in lammps”. In: Computational Materials Science 175 (2020), p. 109584. doi: 10.1016/ j.commatsci.2020.109584. [141] G. Grams et al. “Skyrme–Hartree–Fock–Bogoliubov mass models on a 3D mesh: IV. Improved description of the isospin dependence of pairing”. In: The European Physical Journal A 61.2 (2025), pp. 1–13. [142] M. Warda et al. “Self-consistent calculations of fission barriers in the Fm region”. In: Phys. Rev. C 66.1 (1 July 2002), p. 014310. doi: 10.1103/PhysRevC.66.014310. [143] N. Dubray, H. Goutte, and J.-P. Delaroche. “Structure properties of 226Th and 256,258,260Fm fission fragments: Mean-field analysis with the Gogny force”. In: Phys. Rev. C 77 (1 Jan. 2008), p. 014310. doi: 10.1103/PhysRevC.77.014310. [144] K. H. Schmidt and B. Jurado. “Review on the progress in nuclear fission: experimen- tal methods and theoretical descriptions”. In: Rep. Prog. Phys. 81.10 (Sept. 2018), p. 106301. doi: 10.1088/1361-6633/aacfa7. 96 [145] N. Schunck, D. Duke, and H. Carr. “Description of induced nuclear fission with Skyrme energy functionals. II. Finite temperature effects”. In: Phys. Rev. C 91 (3 Mar. 2015), p. 034327. doi: 10.1103/PhysRevC.91.034327. [146] Y. Abe et al. “On stochastic approaches of nuclear dynamics”. In: Physics reports 275.2-3 (1996), pp. 49–196. [147] M. R. Mumpower et al. “Primary fission fragment mass yields across the chart of nuclides”. In: Phys. Rev. C 101 (5 May 2020), p. 054607. doi: 10.1103/PhysRevC. 101.054607. [148] J. Sadhukhan, W. Nazarewicz, and N. Schunck. “Microscopic modeling of mass and charge distributions in the spontaneous fission of 240Pu”. In: Phys. Rev. C 93 (1 Jan. 2016), 011304(r). doi: 10.1103/PhysRevC.93.011304. [149] K.-H. Schmidt et al. “General description of fission observables: GEF model code”. In: Nucl. Data Sheets 131 (2016). Special Issue on Nuclear Reaction Data, pp. 107– 221. doi: 10.1016/j.nds.2015.12.009. [150] J.-F. Lemaˆıtre et al. “Fission fragment distributions and their impact on the r- process nucleosynthesis in neutron star mergers”. In: Phys. Rev. C 103 (2 Feb. 2021), p. 025806. doi: 10.1103/PhysRevC.103.025806. [151] A. Bulgac, S. Jin, and I. Stetcu. “Nuclear Fission Dynamics: Past, Present, Needs, and Future”. In: Frontiers in Physics 8 (2020). doi: 10.3389/fphy.2020.00063. [152] N. Schunck et al. “Microscopic calculation of fission product yields for odd-mass nuclei”. In: Physical Review C 107.4 (2023), p. 044312. [153] J. Sadhukhan et al. “Efficient method for estimation of fission fragment yields of r-process nuclei”. In: Phys. Rev. C 101 (6 June 2020), p. 065803. doi: 10 . 1103 / PhysRevC.101.065803. [154] J. Sadhukhan, S. A. Giuliani, and W. Nazarewicz. “Theoretical description of fission yields: Toward a fast and efficient global model”. In: Phys. Rev. C 105 (1 Jan. 2022), p. 014619. doi: 10.1103/PhysRevC.105.014619. [155] P.-G. Reinhard et al. “Localization in light nuclei”. In: Phys. Rev. C 83 (3 Mar. 2011), p. 034312. doi: 10.1103/PhysRevC.83.034312. [156] C. L. Zhang, B. Schuetrumpf, and W. Nazarewicz. “Nucleon localization and fragment formation in nuclear fission”. In: Phys. Rev. C 94 (6 Dec. 2016), p. 064323. doi: 10.1103/PhysRevC.94.064323. 97 [157] T. Li et al. “Nucleon localization function in rotating nuclei”. In: Phys. Rev. C 102 (4 Oct. 2020), p. 044305. doi: 10.1103/PhysRevC.102.044305. [158] P. Fong. “Asymmetric fission”. In: Physical Review 89.1 (1953), p. 332. [159] W. D. Myers and W. J. Swiatecki. “Nuclear masses and deformations”. In: Nuclear Physics 81.1 (1966), pp. 1–60. url: https://escholarship.org/uc/item/65v5x6mx. [160] A. Sandulescu. “A new radioactivity”. In: Journal of Physics G: Nuclear and Particle Physics 15.5 (May 1989), p. 529. doi: 10.1088/0954-3899/15/5/008. [161] K. Depta et al. “Bimodal Fission in 258Fm”. In: Mod. Phys. Lett. A 01.06 (1986), pp. 377–381. doi: 10.1142/s0217732386000464. [162] P. M¨oller, J. Nix, and W. Swiatecki. “Calculated fission properties of the heaviest elements”. In: Nucl. Phys. A 469.1 (1987), pp. 1–50. doi: 10.1016/0375- 9474(87) 90083-2. [163] P. M¨oller, J. Nix, and W. Swiatecki. “New developments in the calculation of heavy- element fission barriers”. In: Nucl. Phys. A 492.3 (1989), pp. 349–387. doi: 10.1016/ 0375-9474(89)90403-x. [164] P. M¨oller et al. “Five-dimensional potential-energy surfaces and coexisting fission modes in heavy nuclei”. In: AIP Conf. Proc. 561.1 (Apr. 2001), pp. 455–468. doi: 10.1063/1.1372821. [165] M. Warda et al. “Fission paths in Fm region calculated with the Gogny forces”. In: Phys. At. Nucl. 66 (2003), pp. 1178–1181. doi: 10.1134/1.1586434. [166] M. Warda et al. “Microscopic Structure of The Bimodal Fission of 258Fm”. In: Int. J. of Mod. Phys. E 13.01 (2004), pp. 169–174. doi: 10.1142/s0218301304001904. [167] T. Asano et al. “Dynamical calculation of multi-modal nuclear fission of fermium nuclei”. In: J. Nucl. Radiochem. Sci. 5.1 (2004), pp. 1–5. doi: 10.14494/jnrs2000.5.1. [168] L. Bonneau. “Fission modes of 256Fm and 258Fm in a microscopic approach”. In: Phys. Rev. C 74 (1 July 2006), p. 014301. doi: 10.1103/PhysRevC.74.014301. [169] C. Simenel and A. S. Umar. “Formation and dynamics of fission fragments”. In: Phys. Rev. C 89 (3 Mar. 2014), 031601(r). doi: 10.1103/PhysRevC.89.031601. [170] Y. Miyamoto et al. “Origin of the dramatic change of fission mode in fermium isotopes investigated using Langevin equations”. In: Phys. Rev. C 99 (5 May 2019), 051601(r). doi: 10.1103/PhysRevC.99.051601. 98 [171] D. Regnier, N. Dubray, and N. Schunck. “From asymmetric to symmetric fission in the fermium isotopes within the time-dependent generator-coordinate-method for- malism”. In: Phys. Rev. C 99.2 (Feb. 2019), p. 024611. doi: 10.1103/PhysRevC.99. 024611. [172] M. Albertsson et al. “Super-short fission mode in fermium isotopes”. In: Phys. Rev. C 104 (6 Dec. 2021), p. 064616. doi: 10.1103/PhysRevC.104.064616. [173] K. Pomorski et al. “Fission Fragment Mass And Kinetic Energy Yields Of Fermium Isotopes”. In: Acta Phys. Pol. B 54 (2023), 9–a2. doi: 10.5506/APhysPolB.54.9-A2. [174] R. Bernard, C. Simenel, and G. Blanchon. “Hartree–Fock–Bogoliubov study of quan- tum shell effects on the path to fission in 180Hg, 236U and 256Fm”. In: Eur. Phys. J. A 59.3 (2023), p. 51. [175] R. M. Harbour et al. “Mass and Nuclear Charge Distributions from the Spontaneous Fission of 254Fm”. In: Phys. Rev. C 8 (1973), pp. 1488–1493. doi: 10.1103/PhysRevC. 8.1488. [176] J. E. Gindler et al. “Distribution of mass, kinetic energy, and neutron yield in the spontaneous fission of 254Fm”. In: Phys. Rev. C 16 (4 Oct. 1977), pp. 1483–1492. doi: 10.1103/PhysRevC.16.1483. [177] K. F. Flynn et al. “Distribution of Mass in the Spontaneous Fission of 256Fm”. In: Phys. Rev. C 5 (5 May 1972), pp. 1725–1729. doi: 10.1103/PhysRevC.5.1725. [178] D. C. Hoffman et al. “12.3-min 256Cf and 43-min 258Md and systematics of the spontaneous fission properties of heavy nuclides”. In: Phys. Rev. C 21 (3 Mar. 1980), pp. 972–981. doi: 10.1103/PhysRevC.21.972. [179] J. Berger, M. Girod, and D. Gogny. “Constrained hartree-fock and beyond”. In: Nucl. Phys. A 502 (1989), pp. 85–104. doi: 10.1016/0375-9474(89)90656-8. [180] I. Tsekhanovich et al. “Observation of the competing fission modes in 178Pt”. In: Phys. Lett. B 790 (2019), p. 583. doi: 10.1016/j.physletb.2019.02.006. [181] A. Zdeb, M. Warda, and L. M. Robledo. “Description of the multidimensional potential- energy surface in fission of 252Cf and 258No”. In: Phys. Rev. C 104 (1 July 2021), p. 014610. doi: 10.1103/PhysRevC.104.014610. [182] J. Chi et al. “Role of adecapole deformation in fission potential energy surfaces of 240Pu”. In: Nucl. Phys. A 1032 (2023), p. 122626. doi: 10.1016/j.nuclphysa.2023. 122626. [183] C. Wagemans. “The nuclear fission process”. In: (1991). 99 [184] J. A. Sheikh, W. Nazarewicz, and J. C. Pei. “Systematic study of fission barriers of excited superheavy nuclei”. In: Phys. Rev. C 80 (1 July 2009), p. 011302. doi: 10.1103/PhysRevC.80.011302. [185] M. Arnould, S. Goriely, and K. Takahashi. “The r-process of stellar nucleosynthesis: Astrophysics and nuclear physics achievements and mysteries”. In: Physics Reports 450.4 (2007), pp. 97–213. doi: 10.1016/j.physrep.2007.06.002. [186] F.-K. Thielemann et al. “What are the astrophysical sites for the r-process and the production of heavy elements?” In: Progress in Particle and Nuclear Physics 66.2 (2011). Particle and Nuclear Astrophysics, pp. 346–353. doi: 10.1016/j.ppnp.2011. 01.032. [187] T. Kajino et al. “Current status of r-process nucleosynthesis”. In: Prog. Part. Nucl. Phys. 107 (July 2019), pp. 109–166. doi: 10.1016/j.ppnp.2019.02.008. [188] C. J. Horowitz et al. “r-process nucleosynthesis: connecting rare-isotope beam facilities with the cosmos”. In: J. Phys. G 46.8 (Aug. 2019), p. 083001. doi: 10.1088/1361- 6471/ab0849. [189] Y. L. Zhu et al. “Modeling Kilonova Light Curves: Dependence on Nuclear Inputs”. In: The Astrophysical Journal 906.2 (Jan. 2021), p. 94. doi: 10.3847/1538-4357/abc69e. [190] M.-R. Wu et al. “Fingerprints of Heavy-Element Nucleosynthesis in the Late-Time Lightcurves of Kilonovae”. In: Phys. Rev. Lett. 122.6 (6 Feb. 2019), p. 062701. doi: 10.1103/PhysRevLett.122.062701. [191] J. d. J. Mendoza-Temis et al. “Nuclear robustness of the r process in neutron-star mergers”. In: Phys. Rev. C 92 (5 Nov. 2015), p. 055805. doi: 10.1103/PhysRevC.92. 055805. [192] L. F. Roberts et al. “Electromagnetic Transients Powered By Nuclear Decay In The Tidal Tails Of Coalescing Compact Binaries”. In: The Astrophysical Journal Letters 736.1 (July 2011), p. L21. doi: 10.1088/2041-8205/736/1/l21. [193] S. Goriely. “The fundamental role of fission during r-process nucleosynthesis in neu- tron star mergers”. In: Eur. Phys. J. A 51.2 (Feb. 2015), p. 22. doi: 10.1140/epja/ i2015-15022-3. [194] S. Manzhos and T. J. Carrington. “Neural Network Potential Energy Surfaces for Small Molecules and Reactions”. In: Chem. Rev. 121.16 (2021), pp. 10187–10217. doi: 10.1021/acs.chemrev.0c00665. 100 [195] R. Nagai, R. Akashi, and O. Sugino. “Completing density functional theory by ma- chine learning hidden messages from molecules”. In: npj Comput. Mater. 6.1 (2020). doi: 10.1038/s41524-020-0310-0. [196] B. Jiang, J. Li, and H. Guo. “High-Fidelity Potential Energy Surfaces for Gas-Phase and Gas–Surface Scattering Processes from Machine Learning”. In: J. Phys. Chem. Lett. 11.13 (2020), pp. 5120–5131. doi: 10.1021/acs.jpclett.0c00989. [197] P. O. Dral et al. “Hierarchical machine learning of potential energy surfaces”. In: J. Chem. Phys. 152.20 (2020), p. 204110. doi: 10.1063/5.0006498. [198] S. Akkoyun et al. “Consistent empirical physical formulas for potential energy curves 66Ti isotopes by using neural networks”. In: Phys. Part. Nucl. 10.6 (Nov. 2013), of 38 pp. 528–534. doi: 10.1134/s1547477113060022. − [199] M. Verriere et al. “Building surrogate models of nuclear density functional theory with Gaussian processes and autoencoders”. In: Front. in Phys. 10 (2022). doi: 10. 3389/fphy.2022.1028370. [200] R.-D. Lasseri et al. “Taming Nuclear Complexity with a Committee of Multilayer Neural Networks”. In: Phys. Rev. Lett. 124.16 (16 Apr. 2020), p. 162502. doi: 10. 1103/PhysRevLett.124.162502. [201] X. Glorot and Y. Bengio. “Understanding the difficulty of training deep feedforward neural networks”. In: Proceedings of the Thirteenth International Conference on Ar- tificial Intelligence and Statistics. Ed. by Y. W. Teh and M. Titterington. Vol. 9. Pro- ceedings of Machine Learning Research. Chia Laguna Resort, Sardinia, Italy: Pmlr, 13–15 May 2010, pp. 249–256. url: https://proceedings.mlr.press/v9/glorot10a.html. [202] H. B. Curry. “The method of steepest descent for non-linear minimization problems”. In: Q. of Appl. Math. 2 (1944), pp. 258–261. url: https://api.semanticscholar.org/ CorpusID:125304075. [203] A. Rigler, J. Irvine, and T. Vogl. “Rescaling of variables in back propagation learning”. In: Neural Netw. 4.2 (1991), pp. 225–229. doi: 10.1016/0893-6080(91)90006-q. [204] Z. Matheson et al. “Cluster radioactivity of 294Og”. In: Phys. Rev. C 99.4 (Apr. 2019), 041304(r). doi: 10.1103/PhysRevC.99.041304. [205] D. R. Phillips et al. “Get on the BAND Wagon: a Bayesian framework for quantifying model uncertainties in nuclear dynamics”. In: Journal of Physics G: Nuclear and Particle Physics 48.7 (May 2021), p. 072001. doi: 10.1088/1361-6471/abf1df. [206] A. Gelman et al. Bayesian data analysis. Chapman and Hall/CRC, 1995. 101 [207] S. L. Brunton and J. N. Kutz. Data-Driven Science and Engineering: Machine Learn- ing, Dynamical Systems, and Control. Cambridge University Press, 2019. [208] W. H. Schilders, H. A. Van der Vorst, and J. Rommes. Model order reduction: theory, research aspects and applications. Vol. 13. Springer, 2008. [209] D. Frame et al. “Eigenvector continuation with subspace learning”. In: Physical review letters 121.3 (2018), p. 032501. [210] T. Duguet et al. “Colloquium: Eigenvector continuation and projection-based emula- tors”. In: Rev. Mod. Phys. 96 (3 Aug. 2024), p. 031002. doi: 10.1103/RevModPhys. 96.031002. [211] E. Bonilla et al. “Training and projecting: A reduced basis method emulator for many-body physics”. In: Phys. Rev. C 106 (5 Nov. 2022), p. 054322. doi: 10.1103/ PhysRevC.106.054322. [212] P. Cook et al. “Parametric Matrix Models”. In: (2025). url: https://arxiv.org/abs/ 2401.11694. [213] J. Wang. Geometric structure of high-dimensional data. Springer, 2012. [214] A. Quarteroni, A. Manzoni, and F. Negri. Reduced basis methods for partial differential equations: an introduction. Vol. 92. Springer, 2015. [215] J. P. Boyd. Chebyshev and Fourier spectral methods. Courier Corporation, 2001. [216] M. Barrault et al. “An ‘empirical interpolation’ method: application to efficient reduced- basis discretization of partial differential equations”. In: Comptes Rendus Mathema- tique 339.9 (2004), pp. 667–672. doi: 10.1016/j.crma.2004.08.006. [217] S. Chaturantabut and D. C. Sorensen. “Nonlinear Model Reduction via Discrete Empirical Interpolation”. In: SIAM Journal on Scientific Computing 32.5 (2010), pp. 2737–2764. doi: 10.1137/090766498. [218] Y. Maday, O. Mula, and G. Turinici. “Convergence analysis of the Generalized Em- pirical Interpolation Method”. In: SIAM Journal on Numerical Analysis 54.3 (2016), pp. 1713–1731. doi: 10.1137/140978843. [219] K. Godbey et al. Dimensionality reduction in nuclear physics. url: https://dr.ascsn. net/landing.html. [220] V. Olshevsky. Matrix methods: theory, algorithms and applications. World Scientific, 2010. 102 [221] D. Odell et al. “ROSE: A reduced-order scattering emulator for optical models”. In: Phys. Rev. C 109 (4 Apr. 2024), p. 044612. doi: 10.1103/PhysRevC.109.044612. [222] C. Drischler et al. “BUQEYE guide to projection-based emulators in nuclear physics”. In: Frontiers in Physics 10 (2023). doi: 10.3389/fphy.2022.1092931. [223] L. P. Pitaevskii. “Vortex lines in an imperfect Bose gas”. In: Sov. Phys. JETP 13.2 (1961), pp. 451–454. [224] R. Courant and D. Hilbert. Methods of mathematical physics. CUP Archive, 1985. [225] M. D. McKay, R. J. Beckman, and W. J. Conover. “A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code”. In: Technometrics 21.2 (1979), pp. 239–245. url: http : / / www . jstor . org / stable/1268522 (visited on 02/17/2025). [226] N. Hizawa, K. Hagino, and K. Yoshida. “Analysis of a Skyrme energy density func- tional with deep learning”. In: Phys. Rev. C 108 (3 Sept. 2023), p. 034311. doi: 10.1103/PhysRevC.108.034311. [227] P. M¨oller et al. “Heavy-element fission barriers”. In: Phys. Rev. C 79 (6 June 2009), p. 064304. doi: 10.1103/PhysRevC.79.064304. [228] J. Wilson et al. “Angular momentum generation in nuclear fission”. In: Nature 590.7847 (2021), pp. 566–570. [229] T. Johnson. “The First Indication of Neutrino-Induced Nuclear Fission”. English. Copyright - Database copyright ProQuest LLC; ProQuest does not claim copyright in the individual underlying works; Last updated - 2025-01-17. PhD thesis. 2024, p. 397. url: https : / / www . proquest . com / dissertations - theses / first - indication - neutrino-induced-nuclear-fission/docview/3156296008/se-2. [230] M. Verriere, N. Schunck, and D. Regnier. “Microscopic calculation of fission product yields with particle-number projection”. In: Phys. Rev. C 103 (5 May 2021), p. 054602. doi: 10.1103/PhysRevC.103.054602. [231] E. Flynn. “Nuclear Fission and Tunneling Phenomena”. PhD thesis. Michigan State University, in prep. 103