APPLICATION OF DENSITY FUNCTIONAL THEORY IN NUCLEAR STRUCTURE By Tong Li 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 2022 ABSTRACT APPLICATION OF DENSITY FUNCTIONAL THEORY IN NUCLEAR STRUCTURE By Tong Li The nuclear density functional theory (DFT) is a microscopic self-consistent framework suitable for describing heavy nuclei and performing large-scale studies. In this dissertation I discuss my research works on the development and application of the Skyrme nuclear-DFT framework, covering a broad range of topics, including the nucleon localization in rotating systems, the origin of reflection-asymmetric deformations, the parameter calibration for beta decays, and the development of a new coordinate-space DFT solver. The nucleon localization function (NLF), discussed in the first part, is a useful tool for the visualization of structure information. It has been utilized to characterize clustering and shell structure. How the NLF pattern evolves in rotating systems, how it visualizes internal nuclear structure, and how it is connected with single-particle (s.p.) orbits are discussed in this dissertation. The second part deals with nuclei having reflection-asymmetric shapes, which are important candidates for the search of permanent electric dipole moments. In this dissertation, the origin of pear-like deformation is investigated through both the multipole expansion of the energy density functional and the spectrum of canonical s.p. states. The- oretical predictions of beta-decay rates are discussed next; they are important for r-process simulations that involves nuclei whose experimental beta-decay data are unknown. To pro- vide reliable predictions with quantified uncertainties, the χ2 optimization is performed to constrain parameters that significantly affect beta-decay transitions in proton-neutron finite- amplitude-method calculations. Besides a well calibrated functional, a reliable and efficient DFT solver is also crucial. The Hartree-Fock-Bogoliubov (HFB) method in the coordinate space is preferred for deformed and weakly bound nuclei, as solvers based on basis expansions often have difficulty correctly describing continuum effects. A new HFB solver based on the canonical-basis HFB formalism in the three-dimensional coordinate space is developed in this dissertation. It is a well parallelized solver and has been carefully benchmarked against other established HFB solvers. ACKNOWLEDGMENTS I would like to first thank my advisor, Prof. Witold Nazarewicz. He is nice and always enthusiastic about research, which inspires me a lot during my PhD study. He always does his best to help me solve troublesome issues in both research and life. I am also grateful for H. Metin Aktulga, Shanker Balasubramaniam, Alexandra Gade, and Heiko Hergert, who serve as members of my PhD guidance committee. I would also thank my advisor for my undergraduate thesis, Prof. Furong Xu at Peking University. He brought me into the field of theoretical and computational nuclear physics and helped me a lot on my career development. This dissertation could not be finished without my research collaborators. I would like to thank Shrijita Bhattacharya, Mengzhi Chen, Jacek Dobaczewski, Jon Engel, Samuel Giuliani, Vojtech Kejzlar, Markus Kortelainen, Tapabrata Maiti, Evan Ney, Paul-Gerhard Reinhard, Bastian Schütrumpf, Mookyong Son, Stefan Wild, and Chunli Zhang. It has been pleasant to collaborate with all of them. I want to express special thanks to my closest collaborators, i.e., Chunli Zhang, Mengzhi Chen, Evan Ney and Mookyong Son; discussions with them are really helpful for my research. I would love to thank Kyle Godbey and Pablo Giuliani for their comments and suggestions on the draft of this dissertation. I also sincerely appreciate help offered by other members of my research group, including Sylvester Agbemava, Yuchen Cao, Eric Flynn, Kevin Fossez, Yannen Jaganathen, Rahul Jain, Daniel Lay, Xingze Mao, Zachary Matheson, Eric Olson, Simin Wang, and Joshua Wylie. Besides, help from my office mate Pierre Nzabahimana is also appreciated. Apart from work, I also enjoy the time with Xi Dong, Jie Fu, Chenglong He, Zhen Li, Hao Lin, Avik Sarkar, Chunqiang Xu, and Heda Zhang; I am grateful for their company as my friends. iv I want to thank professors and instructors I encountered in my PhD coursework, such as Martin Berz, Scott Bogner, Alex Brown, Morten Hjorth-Jensen, Mark Iwen, Michael Murillo, Filomena Nunes, Brian O’Shea, Scott Pratt, Andrea Shindler, Mohsen Zayernouri, and Vladimir Zelevinsky. I would like to extend my thanks to supporting staff at Department of Physics & Astronomy and NSCL/FRIB, especially Kim Crosslan, Elizabeth Deliyski, and Gillian Olsen. Supports from my family and friends in China are also appreciated, and it is always comforting to know that someone far away is caring about me. There are many other people who helped me during my PhD study at Michigan State University, and it is a pity that I cannot list all of them in the acknowledgment. Last but not least, I would like to thank everybody who works diligently to keep ev- erything working smoothly during the COVID-19 pandemic. It has been a hard time for everyone, and research could not continue without everyone’s effort. Computational resources for projects discussed in this dissertation were provided by the Institute for Cyber-Enabled Research at Michigan State University, the Laboratory Comput- ing Resource Center (LCRC) and the General Computing Environment (GCE) at Argonne National Laboratory. v TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x LIST OF ALGORITHMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Nuclear density functional theory (DFT) . . . . . . . . . . . . . . . . . . . . 1 1.2 Deformation, rotation and nucleonic localization . . . . . . . . . . . . . . . . 2 1.3 Model calibration for beta-decay calculations . . . . . . . . . . . . . . . . . . 3 1.4 Reliable and efficient DFT solver in coordinate space . . . . . . . . . . . . . 5 1.5 Organization of this dissertation . . . . . . . . . . . . . . . . . . . . . . . . . 6 Chapter 2 Theoretical framework . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 Skyrme energy density functional . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.1 Local densities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.2 Skyrme EDF formalism . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Hartree-Fock-Bogoliubov method . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2.1 HFB equations in the quasiparticle basis . . . . . . . . . . . . . . . . 15 2.2.2 HFB equations in the canonical basis . . . . . . . . . . . . . . . . . . 18 2.2.3 Constrained calculations . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.2.4 Numerical solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.3 Charge-changing finite amplitude method . . . . . . . . . . . . . . . . . . . . 27 2.3.1 FAM equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.3.2 FAM response function . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.3.3 Beta-decay rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.3.4 Gamow-Teller resonance . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.3.5 Numerical PNFAM solver . . . . . . . . . . . . . . . . . . . . . . . . 37 Chapter 3 Nucleon localization function in rotating nuclei . . . . . . . . . . . . . 39 3.1 Nucleon localization function . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.1.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.1.2 Alternative interpretation . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2 Cranked calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2.1 Cranked Hartree-Fock calculation . . . . . . . . . . . . . . . . . . . . 43 3.2.2 Cranked harmonic-oscillator (CHO) calculation . . . . . . . . . . . . 45 3.3 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.3.1 Time-odd local densities . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.3.2 Simplified nucleon localization function . . . . . . . . . . . . . . . . . 49 3.3.3 Dependence on the choice of spin-quantization axis . . . . . . . . . . 52 3.3.4 Angular-momentum alignment in the CHO model . . . . . . . . . . . 52 vi 3.3.5 Angular momentum alignment in the CHF calculation . . . . . . . . 57 3.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 Chapter 4 Origin of reflection-asymmetric shapes . . . . . . . . . . . . . . . . . . 64 4.1 Multipole expansion of the EDF . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.2 Energy expansion and reflection asymmetry . . . . . . . . . . . . . . . . . . 66 4.3 Relation to single-particle spectra . . . . . . . . . . . . . . . . . . . . . . . . 68 4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 Chapter 5 Model calibration for beta-decay studies . . . . . . . . . . . . . . . . . 73 5.1 Model parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.2 Fit observables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.3 Least-squares fit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.4 Numerical details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.5 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.5.1 Optimizations involving one observable type . . . . . . . . . . . . . . 83 5.5.2 Weight determination . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.5.3 Optimizations involving two observable types . . . . . . . . . . . . . 86 5.5.4 Number of effective model parameters . . . . . . . . . . . . . . . . . 87 5.5.5 Sensitivity analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Chapter 6 New HFB solver HFBFFT in three-dimensional coordinate space . . . 93 6.1 Numerical framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.1.1 Numerical realization on a grid . . . . . . . . . . . . . . . . . . . . . 93 6.1.2 Iteration scheme in coordinate space . . . . . . . . . . . . . . . . . . 96 6.1.3 Sub-iteration scheme in the configuration space . . . . . . . . . . . . 97 6.1.4 Convergence criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 6.1.5 Strategies for pairing . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.2 Benchmarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.2.1 Effect of Hermiticity restoration . . . . . . . . . . . . . . . . . . . . . 103 6.2.2 Benchmarks without pairing . . . . . . . . . . . . . . . . . . . . . . . 104 6.2.3 Benchmarks with pairing . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 Chapter 7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 vii LIST OF TABLES Table 5.1: Four optimization schemes with different free and fixed parameter sets. . 76 Table 5.2: GTR energies and their experimental errors (in MeV) taken from Refs. [197, 199, 200, 201] for the four nuclei selected in this work. . . . . . . . 76 Table 5.3: Beta-minus-decay half lives and their experimental errors (in second) taken from Ref. [202] for the 25 nuclei selected in this work. Half lives are listed in ascending order. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Table 5.4: Optimal parameter values obtained from fits with GTR energies only. Stan- dard deviations are given in brackets. The first letter of the fit indicates the scheme (see Table 5.1). Schemes B and D are absent as gA does not affect the GTR energy. Root-mean-square errors (RMSEs) are also presented for comparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 Table 5.5: Similar to Table 5.4 but with beta-decay half-life data only. Optimization schemes C and D with g1′ as a free parameter are absent because of severe numerical instability. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 Table 5.6: Optimal parameter values obtained from fits with weights wGTR = 0.3 MeV and wβ = 1 (relative weight γ = 0.3 MeV). Standard deviations are given in brackets, and the χ2 values are also presented for comparison. . 86 Table 5.7: Correlation matrix (5.7) obtained from Fit A in Table 5.6. . . . . . . . . . 87 Table 5.8: Similar to Table 5.7 but for Fit B given in Table 5.6. Correlations larger than 0.6 in absolute value are marked in italics. . . . . . . . . . . . . . . . 87 Table 5.9: Correlation matrix (5.7) obtained from Fit C in Table 5.6. . . . . . . . . . 88 Table 5.10: Similar to Table 5.9 but for Fit D given in Table 5.6. . . . . . . . . . . . 88 Table 6.1: Parameters adopted in different solvers for benchmark calculations: HF- BTHO – the number of HO shells NHO and the axial deformation of the HO basis β2 ; Sky1D, Sky2D and HFBFFT – the number of points Ni and grid spacing δi in the direction of i; Common parameter – the pairing cutoff energy Ecut . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 Table 6.2: Total energies Etot (in MeV) and ∆S (in MeV) obtained from HFBFFT without and with Hermiticity restoration. Digits that do not coincide are marked in bold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 viii Table 6.3: Total energies Etot of 132 Sn and 208 Pb obtained from HFBTHO, HFBFFT, Sky1D and Sky2D. Various energy components are also listed: Subscript of densities denotes contributions from corresponding ED terms in Eq. (2.17), “Coul” the Coulomb energy, and “kin” the kinetic energy. Digits that do not coincide with HFBFFT are marked in bold. . . . . . . . . . . . . . . . 104 Table 6.4: Results of spherical 120 Sn obtained from HFBTHO, HFBFFT, Sky1D and Sky2D, including total energies Etot , some energy components, average pairing gaps ∆ and total root-mean-square (rms) radii rrms (in fm). All energies are in MeV. Pairing strengths Vn and Vp adopted in these calcula- tions are equal and listed in the last row. Digits that do not coincide with HFBFFT are marked in bold and numbers used for the pairing renormal- ization are in italics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Table 6.5: Results of axially deformed 102 Zr and 110 Zr obtained from HFBTHO, HF- BFFT and Sky2D. Besides quantities shown in Table 6.4, quadrupole mo- ments Q20 (in fm2 ) are also listed here. Pairing strengths Vn and Vp adopted in these calculations are listed in the last two rows. Digits that do not coincide with HFBFFT are marked in bold and numbers used for the pairing renormalization are in italics. . . . . . . . . . . . . . . . . . . . . . 107 Table 6.6: Similar to Table 6.5, but for the ground state and fission isomer of 240 Pu. 107 ix LIST OF FIGURES Figure 3.1: Single-particle neutron (a) and proton (b) Routhians as functions of ω, ob- tained from the CHF+SkM* calculations for the SD yrast band of 152 Dy. The (π, ry ) combinations are denoted by solid lines (+, +i), dotted lines (+, −i), dot-dashed lines (−, +i), and dashed lines (−, −i). Single-particle Routhians of the lowest neutron N = 7, proton N = 6, and proton [541]1/2 levels are marked by thicker lines. Large SD gaps at Z = 66 and N = 86 can be clearly seen in this figure. . . . . . . . . . . . . . . . . . 45 Figure 3.2: Single-particle Routhians of the SD CHO model belonging to supershells Nshell = 6 and 7. The CHO quantum numbers [n1 , n2 , n3 ] are given in brackets. Positive-parity and negative-parity orbits are marked by solid and dashed lines, q respectively. The rotational frequency ω is expressed in units of ω0 = 3 ωz ω⊥ 2 while the Routhians E in units of ℏω . Each level is z doubly degenerate due to the two possible spin orientations. The crossing between the lowest N = 7 Routhian [0,0,7] and the [3,0,0] Routhian at ω/ω0 ⪅ 0.2 is denoted by the arrow. . . . . . . . . . . . . . . . . . . . . 46 Figure 3.3: Current density j in the x-z (y = 0) plane obtained from the CHO model, as a function of the rotational frequency ω (in units of ω0 ). The magnitude |j| (in fm−4 ) is shown by color and line thickness. . . . . . . . . . . . . 48 Figure 3.4: Current density j in the x-z (y = 0) plane for neutrons (top) and protons (bottom) in the SD yrast band of 152 Dy, as a function of the rotational frequency ω (in MeV/ℏ). The magnitude |j| (in fm−4 ) is shown by color and line thickness. The FAM result is presented in the leftmost column with a different color range. . . . . . . . . . . . . . . . . . . . . . . . . . 49 Figure 3.5: Spin density s (top) and spin-kinetic density T (bottom) in the x-y (z = 0) plane for neutrons in the SD yrast band of 152 Dy, as functions of the rotational frequency ω (in MeV/ℏ). The magnitudes, |s| (in fm−3 ) and |T | (in fm−5 ), are shown by color and line thickness. The FAM results are presented in the leftmost column with a different color range. . . . . 50 Figure 3.6: Spin-current tensor density J · ey in the x-z (y = 0) plane for neutrons in the SD yrast band of 152 Dy, as a function of the rotational frequency ω (in MeV/ℏ). Its magnitude (in fm−4 ) is shown by color and line thickness. 50 Figure 3.7: C (top), C τ (middle), and their differences (bottom) in the x-z (y = 0) plane obtained from the CHO model, as functions of the rotational fre- quency ω (in units of ω0 ). . . . . . . . . . . . . . . . . . . . . . . . . . . 51 x Figure 3.8: C (left), C τ (middle), and their difference (right) in the x-z (y = 0) plane for neutrons with σ̆y = −1 (y-simplex ry = −1) in the SD yrast band of 152 Dy at rotational frequency ℏω = 0.9 MeV. . . . . . . . . . . . . . . . 52 Figure 3.9: Cqsµ (3.7) and Cqs τ µ (3.14) in the x-z (y = 0) plane for three spin quanti- zation directions µ = x, y, z, obtained from CHF calculations for the SD yrast band of 152 Dy at rotational frequency ℏω = 0.5 MeV. The symbols ↑ and ↓ represent sµ = +1 and −1, respectively. . . . . . . . . . . . . . 53 Figure 3.10: Similar to Fig. 3.9 but shown in the y-z (x = 0) plane. . . . . . . . . . . 54 Figure 3.11: C τ (top), τ (in fm−5 , middle) and τ TF (in fm−5 , bottom) in the x-z (y = 0) plane, obtained from the CHO model. The leftmost column shows the reference plots at ω = 0 while the other columns show the rotational dependence relative to the ω = 0 reference as a function of the rotational frequency ω (in units of ω0 ). . . . . . . . . . . . . . . . . . . . . . . . . 55 Figure 3.12: C τ (thick solid line), τ (solid line), and τ TF (dashed line) in the non- rotating HO model along the z axis (x = y = 0). (a) 3D SD HO potential with 60 particles. The density profile of the [0,0,6] orbit is marked by a dotted line. (b) 1D case. HO orbits with quantum number N ≤ 6 are occupied. The density profile of the N = 6 orbit is marked by a dotted line; here τ TF = π 2 ρ3 /3. Some quantities are scaled for better visualization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Figure 3.13: Variations in the kinetic density τ due to p-h excitations (at ω = 0) from supershell Nshell = 6 to Nshell = 7 shown in Fig. 3.2. These excitations are induced by the cranking term in the CHO model. The rightmost panel displays the average of all the p-h contributions. . . . . . . . . . . . . . 57 Figure 3.14: The NLF Cqσ̆ τ in the x-z (y = 0) plane as a function of ω (in MeV/ℏ), y obtained from CHF calculation for the SD yrast band of 152 Dy. The symbols ↑ and ↓ represent σ̆y = +1 and −1 (y-simplex ry = +i and −i), respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 Figure 3.15: Similar to Fig. 3.14 but for ∆Cqσ̆τ . The reference value of C τ at ω = 0 is y shown in the leftmost column of Fig. 3.14. . . . . . . . . . . . . . . . . . 59 Figure 3.16: Similar to Fig. 3.15, but for ∆τqσ̆y (in fm−5 ). The reference value of τqσ̆y at ω = 0 is shown in the leftmost column. . . . . . . . . . . . . . . . . . 60 Figure 3.17: Similar to Fig. 3.16, but for ∆τqσ̆TF (in fm−5 ). . . . . . . . . . . . . . . . 61 y xi Figure 3.18: ∆τ (in fm−5 ) of neutrons (top) and protons (bottom) in the x-z (y = 0) plane for different parity-signature combinations (π, ry ) in the SD yrast band of 152 Dy at rotational frequency ℏω = 0.1 MeV. . . . . . . . . . . . 62 Figure 3.19: Contributions to the variation ∆τ (in fm−5 ) in the x-z (y = 0) plane for different parity-signature blocks from individual s.p. Routhians in 152 Dy at ℏω = 0.1 MeV: the four positive-parity neutron levels [651]1/2, [642]5/2, [413]5/2 and [411]1/2 with signature ry = +i (a) and −i (b) that lie below the N = 86 shell gap in Fig. 3.1 (see Fig. 1 of Ref. [152] for the asymptotic quantum numbers [N nz Λ]Ω of s.p. levels in the SD 152 Dy); the N = 7 neutron intruder states 71 (c) and 72 (d); the N = 6 proton intruder states 62 + 64 (e) and 61 + 63 (f); and the [541]1/22 (g) and [541]1/21 (h) proton states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 Figure 4.1: The octupole deformation energy ∆E of 224 Ra (a), its multipole expan- sion ∆E[L] (b), and the components of the octupole part ∆E[3] in isospin and proton-neutron representations (c), as functions of the octupole de- formation β3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Figure 4.2: The octupole deformation energy ∆E (a), its multipole components ∆E[L] (b), and the cumulative sum ∆E[0−L] (c), along the isotopic chain of Ra. They are calculated at a fixed octupole deformation β3 = 0.05. Isotopes with negative total deformabilities (octupole-deformed ground states) are marked by shading. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 Figure 4.3: Similar to Fig. 4.2 but along the isotopic chain of Yb. . . . . . . . . . . 67 Figure 4.4: Canonical single-particle energies for neutrons (top) and protons (bottom) in 224 Ra as functions of the quadrupole deformation β2 . The reflection symmetry is imposed (β3 = 0). Solid (dashed) lines represent levels with positive (negative) parity. The equilibrium quadrupole deformation of 224 Ra is given by a vertical dotted line, and the Fermi energies of 224 Ra are marked by dash-dotted lines. Quadrupole deformations and Fermi energies of even-even Ra isotopes with neutron number N = 130 ∼ 144 are denoted by circles. All the Fermi energies are shifted with respect to the positions of spherical neutron 2g9/2 and proton 1h9/2 shells. The line color indicates the angular-momentum projection Ω. . . . . . . . . . . . 71 Figure 4.5: Similar to Fig. 4.4 but for 176 Yb. . . . . . . . . . . . . . . . . . . . . . 72 Figure 5.1: Values of w̃ and r for fits with a series of relative weights γ = wGTR /wβ . The optimization scheme A in Table 5.1 is used and all the data in Tables 5.2 and 5.3 are included in the fits. Actual calculations are performed with wβ = 1 and wGTR = γ. . . . . . . . . . . . . . . . . . . . . . . . . 84 xii Figure 5.2: Relation between the RMSEs of two types of data (GTR energies EGTR and logarithms of beta-decay half lives lg T1/2 ), obtained from the fits shown in Fig. 5.1. Relative weights γ (in MeV) are represented by the colors of square markers. . . . . . . . . . . . . . . . . . . . . . . . . . . 85 Figure 5.3: Individual weighted residuals obtained from the four fits presented in Ta- ble 5.6. The labeling of data points is consistent with that in Tables 5.2 and 5.3. The vertical dashed line separates the two data types, and the horizontal dashed line indicates zero. . . . . . . . . . . . . . . . . . . . 86 Figure 5.4: Squared singular values (in descending order) obtained from Eq. (5.12), the SVD of the normalized Jacobian matrix J, ˇ for the four fits presented in Table 5.6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 Figure 5.5: Cumulative fractions Sm (5.13) for the four fits presented in Table 5.6. The gray horizontal line indicates the 0.99 threshold. . . . . . . . . . . . 89 Figure 5.6: Squared amplitudes of the first (left) and second (right) principal com- ponents (column vectors of V̌ ) that correspond to the first and second largest singular values, for the four fits presented in Table 5.6. They are obtained from Eq. (5.12), the SVD of the normalized Jacobian matrix J. ˇ 90 Figure 5.7: Normalized sensitivity matrices (5.11) in absolute value for the fits pre- sented in Table 5.6, with weights wGTR = 0.3 MeV and wβ = 1 (relative weight γ = 0.3 MeV). The labeling of data points is consistent with that in Tables 5.2 and 5.3. The vertical dashed line separates the points of GTR energies and the logarithms of beta-decay half lives. . . . . . . . . 91 xiii LIST OF ALGORITHMS m Algorithm 6.1: Compute the 1D differentiation ψ (m) ≡ ddxmψ , m = 1, 2, 3, · · · . . . . . 95 h i Algorithm 6.2: Compute the 1D position-varying differentiation χ ≡ dx d B(x) dψ . . 96 dx Algorithm 6.3: Damped gradient iteration scheme of HFBFFT. . . . . . . . . . . . . 98 xiv Chapter 1 Introduction 1.1 Nuclear density functional theory (DFT) Being a many-body quantum system, the atomic nucleus is difficult to describe with theoret- ical and computational models. The nuclear density functional theory (DFT) [1, 2] provides a microscopic self-consistent mean-field model suitable for studying medium- to heavy-mass nuclei and for performing global surveys over the whole nuclear landscape. In a mean-field model, every nucleon moves approximately in a mean potential generated by other nucle- ons, and the concept of the mean field is justified by shell effects observed experimentally (e.g., the existence of magic numbers); see Ref. [3] for a comprehensive discussion. This dissertation discusses my research on the applications of the nuclear DFT. In theoretical and computational studies of nuclear structure, it is of great importance to choose appropriate degrees of freedom and approximations to achieve good descriptions with computational expenses under control. From this perspective, the nuclear DFT model lies between the shell model [4, 5] and microscopic-macroscopic (mic-mac) method [6, 7, 8, 9, 10]. Compared with the DFT model, the mic-mac method is computationally inexpensive but not reliable enough for extrapolations beyond experimentally accessible region, while the ab initio approach [11] and shell model can provide beyond-mean-field descriptions but are too computationally expensive for heavy systems and large-scale surveys. Hence, the nuclear DFT is usually the appropriate choice for such studies. Within the nuclear-DFT framework, the energy of a system is given by the energy density 1 functional (EDF), and there are three functional forms that are widely used – Skyrme [12, 13, 14, 15], Gogny [16, 17, 18] and relativistic [19, 20] EDFs. Inside each category, there also exist variations and different parameterizations. This dissertation focuses on the application of the Skyrme DFT; the underlying theoretical framework is presented in Chapter 2. 1.2 Deformation, rotation and nucleonic localization Nuclear deformation comes from spontaneous symmetry breaking related to the nuclear Jahn-Teller effect [21, 22, 23, 24], a concept originally proposed in molecular physics to explain the geometric distortion of non-linear molecules. Although the vast majority of isotopes exhibit spherical, prolate or oblate ground-state shapes, there is abundant experimental evidence supporting the existence of stable pear- like deformations which break the reflection symmetry in the intrinsic frame [25, 26]. As even-even nuclei with reflection-asymmetric shapes usually have low-energy negative-parity excitations related to octupole collective modes, they are also referred to as “octupole de- formed.” Atoms with pear-shaped nuclei have also been shown to be good laboratories for searching the permanent electric dipole moment [27, 28, 29], which is strongly correlated with the nuclear octupole moment [30]. Within the nuclear-DFT framework, a global survey for reflection-asymmetric ground states of even-even nuclei has recently been conducted in Ref. [31]. Following this work, we investigate the microscopic origin of reflection-asymmetric shapes from two perspectives: the multipole expansion of the EDF and single-particle (s.p.) spectra. An in-depth discussion on this topic is given in Chapter 4. A number of nuclear collective modes, such as rotations and vibrations, are related to nuclear deformation. The observation of rotational bands provides abundant information about the nuclear ground-state shape and underlying shell structures, as well as the interplay 2 between collective and s.p. degrees of freedom [32, 33, 34, 35]. This interplay results from the fact that one cannot distinguish slow and fast wave-function components in nuclear systems as is usually done in molecules. Perfect nuclear rotors are rare and the adiabatic approximation is usually inapplicable; thus, the fully self-consistent and microscopic nuclear- DFT model is a good choice for the descriptions of nuclear rotation. Recently, the nucleon localization function (NLF) has been employed in the nuclear- DFT studies for better visualization of clusters and shell effects. It was originally introduced in electronic-DFT studies to characterize shell structure in atoms and chemical bonds in molecules [36, 37, 38, 39]. In nuclear physics, the NLF has been proved to be a useful tool for the identification and description of (i) clusters in light nuclei [40, 41, 42] and heavy-ion collisions [43], (ii) fragment formation along the fission pathway [44, 45, 46, 47, 48, 49], (iii) nuclear pasta phases in the inner crust of neutron stars [41], and (vi) shell structures of electrons and nucleons in Oganesson [50]. Compared with the particle distribution that is quite smooth in the nuclear interior, the NLF clearly displays the internal structure of a nucleus through its characteristic oscillating or cluster patterns. In Chapter 3 I discuss how the NLF characterizes the nuclear response to rotation and show why the NLF constitutes a powerful visualization tool for nuclear structure studies. Based on my work, the concept of the Pauli kinetic energy is brought up in Ref. [51], which quantifies the effect of the Pauli exclusion principle in the heavy-ion collision. 1.3 Model calibration for beta-decay calculations How heavy elements are created in the universe is a fascinating open question [52]. Various astrophysical processes contribute to the synthesis of these elements, and both slow and rapid neutron-capture processes (s and r processes) are believed to be predominant [53, 54, 3 55, 56, 57]. The study of these processes requires ample astronomical observations as well as data on the properties of nuclei involved. Theoretical predictions of beta-decay rates are crucial nuclear inputs for the simulations of some astrophysical processes, especially the r process that traverses the neutron-rich region that experiments cannot currently access. The competition between neutron captures, photodissociation, and beta decays determines how the r process proceeds in different astronomical environments, and the final abundances of stable nuclei are strongly affected by the beta-decay rates of their progenitors [7, 58, 59, 60, 61, 62]. Therefore, the reliability of r-process simulations depends heavily on the quality of beta-decay predictions. Within the nuclear-DFT framework, the finite amplitude method (FAM) is an efficient approach for the solution of the (quasiparticle) random phase approximation (QRPA) [63, 64]. Thanks to recent developments of the Skyrme proton-neutron FAM (PNFAM), it is now feasible to calculate charge-changing transitions in deformed nuclei and to conduct large-scale calculations for beta-decay rates [62, 65, 66, 67]. In the PNFAM, the time-odd isovector Skyrme couplings, isoscalar pairing strength, and effective axial-vector coupling have a strong impact on beta-decay transitions, but they are not constrained by ground- state properties of even-even nuclei and should thus be calibrated based on experimental data related to beta decays. In this work we carry out the χ2 optimization for the model calibration, following the procedure discussed in Refs. [68, 69]. After the model calibration we will be able to provide reliable predictions with quantified uncertainties for r-process simulations. Details pertaining to this topic are presented in Chapter 5. 4 1.4 Reliable and efficient DFT solver in coordinate space Exotic nuclei, which are far from the valley of stability and hence weakly bound, provide rich information that helps us improve nuclear models; they also play important roles in a number of astrophysical processes. Within the nuclear-DFT framework, the ground state of an even-even nucleus is solved by the Hartree-Fock (HF) + Bardeen-Cooper-Schrieffer (BCS) or Hartree-Fock-Bogoliubov (HFB) methods. The former method is simpler but can be problematic in exotic nuclei, and the full HFB scheme that treats nuclear pairing in a fully self-consistent way should instead be employed for such systems [70, 71, 72, 73, 74]. The HFB calculations of exotic nuclei, however, can be computationally expensive, as a large model space is needed for weakly bound s.p. states; these s.p. states have broad density distributions and are strongly affected by the continuum. The problem becomes even more severe for symmetry-unrestricted calculations. There have been a number of HF+BCS and HFB solvers developed for nuclear-DFT studies, which can be divided into two categories: Some are formulated in the coordinate- space representation while others are based on basis expansions of wave functions. Table 2 in Ref. [75] gives a summary of these solvers. For exotic nuclei, the coordinate-space HFB solver is a preferred choice as basis-based solvers often have difficulty producing correct asymptotic behaviors of weakly bound s.p. orbits. We have developed a reliable and efficient Skyrme- HFB solver HFBFFT in the three-dimensional Cartesian coordinate representation with no spatial symmetry imposed. It is named after the fast-Fourier-transform (FFT) technique used for numerical differentiation. The new solver is based on the canonical-basis HFB formalism proposed in Refs. [76, 77] and finds the HFB solution via the damped gradient method [76, 78, 79, 80, 81]. Compared with the quasiparticle basis, the canonical basis 5 does not have an intractably huge level density as canonical states are spatially localized; the computational cost is thus under control. Furthermore, HFBFFT is well optimized and highly parallelized; it is also benchmarked against other established HFB solvers. Numerical details and benchmark results can be found in Chapter 6. 1.5 Organization of this dissertation This dissertation is organized as follows. Chapter 2 discusses the theoretical formulation of the Skyrme-DFT approach and covers the numerical methods employed in following chapters. The study of the NLF in rotating systems can be found in Chapter 3, which is followed by the investigation about the origin of reflection-asymmetric shapes in Chapter 4. In Chapter 5 I discuss the model calibration for beta-decay calculations, and Chapter 6 describes the new solver HFBFFT. Finally, the conclusions are given in Chapter 7. 6 Chapter 2 Theoretical framework The foundation of the DFT is the Hohenberg-Kohn (H-K) theorem [82] and Kohn-Sham (K-S) equation [83]. Here only the main points of the formalism are discussed, while a comprehensive discussion can be found in Ref. [84]. The H-K theorem states that for a many-body quantum system in an external field vext (x), there exists a unique energy functional Ev [ρ] of the density ρ(x) such that Z X Ev [ρ] = F [ρ] + dx vext (x)ρ(x), (2.1) which is minimized if and only if ρ(x) is the ground-state density. Here x = (r, s, τ ) includes spatial, spin, and isospin coordinates. The F [ρ] part is independent of the external field and thus universal. Although the H-K theorem does not tell how F [ρ] should be constructed, it allows us to search for or to guess an approximate functional. The K-S theory suggests a practical approach to minimize the energy functional: It constructs an auxiliary system of non-interacting particles in the K-S potential, whose density is identical to the true ground state of interacting particles in the external field vext (x). The many-body wave function of the auxiliary system is a Slater determinant, simpler to handle than the true wave function which is a superposition of Slater determinants. The DFT framework was originally developed in the context of electronic systems. In a nuclear system, however, a Slater determinant given by the K-S equation is not enough to account for pairing correlations, so the HFB method has been utilized, in which the nuclear 7 system is described by a product state consisting of non-interacting quasiparticles. In the following I outline the theoretical framework of the Skyrme DFT, which has been widely employed for the global studies of various nuclear properties, such as ground-state energies, deformations, and low-lying excitations [1, 2, 85, 86]. 2.1 Skyrme energy density functional In this section we follow the notation employed in Refs. [1, 15] and define various local densities and the Skyrme EDF. 2.1.1 Local densities Starting from a many-body wave function |Ψ⟩ and a complete s.p. basis {ψ1 , ψ2 , · · · }, one can define the one-body density matrix ρ and pairing tensor κ as † ρij = ⟨Ψ|âj âi |Ψ⟩, κij = ⟨Ψ|âj âi |Ψ⟩, (2.2) † where âi and âi create and annihilate, respectively, a nucleon in the s.p. state ψi . In the coordinate space, the non-local density matrix and pairing tensor can be similarly defined as † ρ rsτ, r ′ s′ τ ′ = ⟨Ψ|â ′ ′ ′ ârsτ |Ψ⟩, κ rsτ, r ′ s′ τ ′ = ⟨Ψ|âr′ s′ τ ′ ârsτ |Ψ⟩, (2.3)   rsτ † where ârsτ and ârsτ create and annihilate, respectively, a nucleon at point r with spin s = ± 12 and isospin τ = ± 21 . 8 In terms of spin and isospin components, the density matrix can be decomposed as  1 ρ rsτ, r ′ s′ τ ′ = ρ00 r, r ′ δss′ + s00 r, r ′ · σ ss′ δτ τ ′    4 +1 (2.4) 1 X h i ρ1t3 r, r ′ δss′ + s1t3 r, r ′ · σ ss′ ⃗ττ τ ′ t    + 4 3 t3 =−1 where σ ss′ = s|σ|s′ and ⃗ττ τ ′ = τ |⃗τ |τ ′ are matrix elements of Pauli matrices in spin and   isospin spaces, respectively. Then we have ρ00 r, r ′ = ρ rsτ, r ′ sτ ,  X (2.5a)  sτ ′ ρ rsτ, r ′ sτ ′ X (2.5b)   ρ1t3 r, r = ⃗ττ ′ τ t , 3 sτ τ ′ s00 r, r ′ = ρ rsτ, r ′ s′ τ σ s′ s ,  X (2.5c)  ss′ τ s1t3 r, r ′ = ρ rsτ, r ′ s′ τ ′ σ s′ s ⃗ττ ′ τ t . X (2.5d)    3 ss′ τ τ ′ Starting from the non-local density matrix, one can define local densities and currents as follows. 1. Time-even densities: ρtt3 (r) = ρtt3 (r, r), (2.6a) τtt3 (r) = ∇ · ∇′ ρtt3 r, r ′ (2.6b)  , r=r ′ i ∇′ − ∇ ⊗ stt3 r, r ′ (2.6c)   Jtt3 (r) = . 2 r=r ′ These are the particle density, kinetic density, and spin-current tensor, respectively. They do not change sign under the time-reversal operation. 9 2. Time-odd densities: stt3 (r) = stt3 (r, r), (2.7a) i ∇′ − ∇ ρtt3 r, r ′ (2.7b)   j tt3 (r) = , 2 r=r ′ T tt3 (r) = ∇ · ∇′ stt3 r, r ′ (2.7c)  ′ , r=r 1 ∇ ⊗ ∇′ + ∇′ ⊗ ∇ · stt3 r, r ′ (2.7d)   F tt3 (r) = . 2 r=r ′ These are the spin density, current, spin-kinetic density, and tensor-kinetic density, respectively. They change sign under the time-reversal operation. For a stationary state of an even-even nucleus, all the time-odd densities vanish due to the time-reversal symmetry. As for the spin-current tensor J, only its anti-symmetric part is usually considered in the Skyrme EDF, which can be written as the spin-orbit current J tt3 = ijk ϵijk ei ej · J · ek , where ei is the unit vector in the i direction and ϵijk is the P  Levi-Civita symbol. Usually the proton-neutron mixing is not considered, which means τ = τ ′ in Eq. (2.3) and t3 = 0 in Eqs. (2.4) ∼ (2.7). In this case we can define density matrices separately for neutrons and protons as ρq (rs, r ′ s′ ) = ρ(rsτq , r ′ s′ τq ), (2.8) where q ∈ {n, p} stands for neutrons or protons and τq = ± 21 for q = n or p. Then the local densities of neutrons and protons are given by ϱ00 = ϱn + ϱp , ϱ10 = ϱn − ϱp , (2.9) 10 where ϱ = ρ, τ, J, J , s, j, T , F . As for the pairing channel, instead of the pairing tensor κ, it is more convenient to formulate the EDF in terms of the pairing density matrix ρ̃ rsτ, r ′ s′ τ ′ = −2s′ κ(rsτ ; r′ , −s′ , τ ′ ), ρ̆ rsτ, r ′ s′ τ ′ = 4s′ τ ′ κ(rsτ ; r′ , −s′ , −τ ′ ), (2.10)   where ρ̃ is usually adopted when there is no proton-neutron mixing, while ρ̆ is more conve- nient to use in the isospin representation. One can decompose ρ̆ rsτ, r ′ s′ τ ′ in the same  way as Eq. (2.4) and then define various local quantities following Eqs. (2.5) ∼ (2.7). But in this dissertation we only care about the time-even pairing density ρ̆ rsτ, rsτ ′ X (2.11)   ρ̆1t3 (r) = ⃗ττ ′ τ t , 3 sτ τ ′ and time-odd pairing spin density ρ̆ rsτ, rs′ τ σ s′ s . X (2.12)  s̆00 (r) = ss′ τ In the case of no proton-neutron mixing, it is more convenient to use the “tilde” density ρ̃ rsτ, r ′ s′ τ ′ and define local neutron and proton pairing densities separately:  ρ̃q (rsτq , r ′ sτq ), q ∈ {n, p}. X ρ̃q (r) = (2.13) ss′ While the particle density ρq (r) gives the probability of finding a nucleon at position r, the pairing density ρ̃q (r) describes the enhancement of the probability of finding a pair of nucleons with opposite spins due to the pairing correlation. 11 2.1.2 Skyrme EDF formalism Within the Skyrme EDF framework, the total energy of a nucleus can be written as an integral over the whole space: Z Z (2.14)   E= dr H(r) = dr Hkin (r) + HSk (r) + Hpair (r) + HCoul (r) + Hcm (r) , where H(r) is the energy density (ED) that consists of various terms that are discussed in the following. The kinetic ED is X ℏ2 Hkin (r) = τq (r). (2.15) 2mq q∈{n,p} 2 Some parameterizations assume mn = mp = m and the kinetic ED then becomes 2m ℏ τ . 00 The Skyrme interaction ED (in the particle-hole channel) is 1 X t X (even) (odd) HSk (r) = Htt (r) + Htt (r), (2.16) 3 3 t=0 t3 =−t where (even) ρ ∆ρ ρ∇J Htt (r) = Ct ρ2tt + Ct ρtt3 ∇2 ρtt3 + Ctτ ρtt3 τtt3 + CtJ J2tt + Ct ρtt3 ∇ · J tt3 (2.17) 3 3 3 is bilinear in time-even local densities, and (odd) j Htt (r) =Cts s2tt + Ct∆s stt3 · ∇2 stt3 + Ct j 2tt + CtT stt3 · T tt3 3 3 3  2 (2.18) s∇j F ∇s + Ct stt3 · ∇ × j tt3 + Ct stt3 · F tt3 + Ct ∇ · stt3 12 is bilinear in time-odd local densities. All the coupling constants C can be density-dependent, ρ but in most Skyrme parameterizations only Ct and Cts depend on the isoscalar particle density ρ00 as ρ00 γ   Ct [ρ00 ] = Ct [0] + (Ct [ρc ] − Ct [0]) , (2.19) ρc where ρc ≈ 0.16 fm−3 is the nuclear saturation density. The Skyrme interaction ED can be either derived from the Hartree-Fock calculation with an effective zero-range momentum- dependent two-body nuclear force proposed by Skyrme [87, 88], or from the density-matrix expansion without reference to an effective force [89, 90]. The first option results in strong dependencies among coupling constants C, while the latter offers more degrees of freedom. Terms involving CtJ , CtT , CtF and Ct∇s in Eqs. (2.17) and (2.18) are known as tensor terms, as they are related to the local two-body tensor interaction [87, 88, 91, 15]. There are a number of publications discussing the effects of tensor terms; see, e.g., Refs. [92, 93, 94] for systematic discussions. The pairing ED (in the particle-particle channel) is usually based on a density-dependent δ force. When the isospin symmetry is preserved, it can be written as   1 2   1 ρ (r)  1 − 00 V0 |s̆00 |2 + V1 X Hpair (r) = ρ̆1t3  , (2.20) 8 ρref t3 =−1 where V0 and V1 are isoscalar and isovector pairing strengths, respectively. When the proton- neutron mixing is absent, we can break the isospin symmetry by assigning different pairing strengths for neutrons and protons:   1 ρ00 (r) |ρ̃q |2 . X Hpair (r) = Vq 1 − (2.21) 4 ρref q∈{n,p} 13 One can adjust the reference density ρref for different types of density dependencies [95, 96]. A pure contact interaction (volume pairing) corresponds to ρref → ∞; a value around the nuclear-matter saturation density ρref = ρc leads to pairing around the nuclear surface (surface pairing); any value in between delivers a mixed-pairing prescription. A widely used value is ρref = 0.32 fm−3 = 2ρc . The Coulomb term is  1 e2 ρp (r)ρp r ′ 3e2 3 3 4 Z  HCoul (r) = dr − ρp3 (r), (2.22) 2 |r − r ′ | 4 π where the second term (exchange term) is given by the Slater approximation. As for the center-of-mass (c.m.) correction Hcm (r), there are several recipes and only those employed in this dissertation are discussed. In Skyrme parameterizations SkM* [97] and SLy4 [98], the c.m. correction is introduced by renormalizing the mass of nucleons:   1 → 1 1 − 1 , where A is the number of nucleons in a nucleus. This mass renormalization m m A is derived from an approximate c.m. correction [86, 99] * A + (diag) 1 p̂2k X Ecm = Ecm = . (2.23) 2m k=1 In UNEDF1-HFB [100], a parameterization designed for the study of fission, however, no c.m. correction is added, because the c.m. correction is not additive in the particle number and causes problems in fission calculations [101]. There are a number of model parameters in the Skyrme EDF, which should be adjusted to reproduce selected experimental data. Thus, there exist various Skyrme parameterizations that are fitted to different data sets with probably different assumptions. The fit data and 14 underlying assumptions are based on the selection of nuclear properties that researchers want to describe well. For example, the energies of fission isomers are included in the fits of UNEDF1 [101] and UNEDF1-HFB [100] so that the nuclear fission can be well described by these two parameterizations. 2.2 Hartree-Fock-Bogoliubov method 2.2.1 HFB equations in the quasiparticle basis Within the nuclear-DFT framework, the ground state of an even-even nucleus is obtained through the variational principle. We choose the product state as the trial wave function and minimize the total energy under various constraints, which leads to the self-consistent HFB equation [1, 2, 3]. A product HFB state |Φ⟩ is a vacuum with respect to quasiparticles: β̂k |Φ⟩ = 0 for all k = 1, 2, 3, · · · , (2.24) where the creation and annihilation operators of quasiparticles are defined by the Bogoliubov transformation: † X †  X  β̂k = Ulk âl + Vlk âl , β̂k = Ulk∗ â + V ∗ ↠. (2.25) l lk l l l This transformation can be written in a compact manner as        † †  β̂  †  â U V   â   =W  =  . (2.26) β̂ † â † V T U T ↠15 The matrix W must be unitary to ensure that β and β † obey the canonical anticommutation relations for fermions. Combining Eqs. (2.25) and (2.2), we can express the density matrix and pairing tensor of the quasiparticle vacuum |Φ⟩ as ρ = V ∗V T , κ = V ∗U T . (2.27) In general, the product state |Φ⟩ does not conserve neutron or proton numbers, so La- grange multipliers need to be introduced to constrain average particle numbers. The total HFB Routhian is then written as (q) D E R [ρ, κ, κ∗ ] E [ρ, κ, κ∗ ] − X = ϵF N̂q , (2.28) q∈{n,p} where N̂q is the particle number operator for neutrons (q = n) or protons (q = p), and (q) the Lagrange multipliers ϵF are also known as chemical potentials or Fermi energies. The minimization of R with regard to ρ and κ yields the HFB equation:        Uk  h − ϵF ∆  Uk  Uk  H =    = Ek   , (2.29) Vk −∆∗ −h∗ + ϵF Vk Vk U  where h is called the HF mean field while ∆ is the pairing field. The eigenvector Vk k is the k-th column of W and the corresponding eigenvalue Ek gives the quasiparticle energy. Uk and Vk are often referred to as the upper and lower components of the quasiparticle wave function, respectively. The matrix elements of H are given by δE [ρ, κ, κ∗ ] δE [ρ, κ, κ∗ ] (q ) hij = = h∗ji , ∆ij = ∗ = −∆ji , (ϵF )ij = δij ϵF i , (2.30) δρji δκij 16 where qi ∈ {n, p} is the isospin of the s.p. state i. One should note that the expression of ϵF will not have the above simple form when there is proton-neutron mixing [15]. Explicit expressions of these matrix elements can be found in a number of publications, e.g., Refs. [1, 15, 86]. Generally speaking, the HFB Hamiltonian depends on the solution U and V , so the HFB equation is non-linear and must be solved in a self-consistent iterative approach (e.g., iterative diagonalization or gradient descent). The HFB equation (2.29) is derived from the variation with respect to the density matrix ρ and pairing tensor κ, but one can also replace κ with the pairing density ρ̃ or ρ̆ and obtain a similar equation that obviously yields the same physics [15, 70]. When pairing collapses in a closed-shell system, κ = 0, ∆ = 0, and Eq. (2.29) reduces to the HF equation that only involves the diagonalization of the HF mean field h. Meanwhile, the product state |Φ⟩ becomes a Slater determinant that has conserved neutron and proton numbers. Therefore, the HF approach can be treated as a special case of the HFB method. Based on the HF method, the pairing correlation can be included via the BCS approximation, and the HF+BCS approach is less computationally expensive than the full HFB method. However, the HF+BCS method can produce an unphysical particle gas surrounding the nucleus and thus a full HFB scheme is preferred [72]. It should be noted that in both HF+BCS and HFB frameworks, the pairing-rearrangement term brought by the density dependence of the pairing functional is also included in the HF mean field h. The HFB equation is numerically solved in a truncated model space, and it is well known that pairing functionals based on the δ interaction, such as those defined in Eqs. (2.20, 2.21), diverge as the model space increases, so a cutoff has to be imposed on the pairing-active space [70, 72, 102]. This cutoff can be done in the quasiparticle space: The contribution of one quasiparticle state to all the densities (including the pairing density or pairing tensor) is 17 multiplied by a function of the quasiparticle energy. A commonly used soft cutoff function is [103, 104] 1 w(e) = , (2.31) e−Ecut  1 + exp ∆Ecut where Ecut is the cutoff energy, ∆Ecut provides smearing around E0 , and e is the reference (equivalent) s.p. energy (q ) ek = (1 − 2Pk ) Ek + ϵF k , (2.32) where Ek is the quasiparticle energy and Pk represents the norm of the lower component Vk . This reference spectrum is similar to the spectrum of the canonical basis; the definition of the canonical basis is to be presented in the next section. When ∆Ecut → 0, w(e) becomes a hard cutoff at energy Ecut . 2.2.2 HFB equations in the canonical basis Besides the formulation based on quasiparticles, the HFB theory can also be formulated in the s.p. basis of canonical states (or natural orbitals), in which the one-body density matrix ρ is diagonal. The connection between the quasiparticle and canonical states is discussed in Refs. [72, 105]. On the one hand, once all the quasiparticle states are obtained, the canonical basis can be computed by diagonalizing the density matrix ρ, and the eigenvalues of ρ represent the occupations of corresponding canonical states. On the other hand, one can directly solve the HFB problem in the canonical basis without any reference to the quasiparticle states. The canonical-basis HFB formalism adopted in my work was first discussed in Ref. [76] and further developed in [77]. The canonical basis is given by a set of orthonormal s.p. wave functions ψα with occu- 18 pation amplitudes vα : {ψα , vα , α = 1, ..., Ω} , vα ∈ [0, 1], (2.33) where vα2 and ψα are the eigenvalues and eigenvectors of the density matrix ρ, respectively, and Ω is the size of the active s.p. space. The non-occupation amplitude is defined as uα = 1 − vα2 . In the canonical basis, the HFB product state takes the BCS-like form p Y  † † |Φ⟩ = uα + vα âα âα |0⟩, (2.34) α>0 † where |0⟩ is the vacuum state, âα generates a particle in the state ψα , and α is the conjugate partner of α that corresponds to the same eigenvalue of ρ. In this section we focus on the stationary state of an even-even nucleus, so the partner α is assumed to be the time-reversed state of α. The density matrix in the coordinate space can be expressed in terms of canonical wave functions: ρ rsτ, r ′ s′ τ ′ =  X 2 vα ψα (r, s, τ )ψα∗ (r ′ , s′ , τ ′ ), (2.35) α from which one can calculate all the local densities involved in the particle-hole channel. Meanwhile, the pairing density matrix is ρ̃ rsτ, r ′ s′ τ ′ = uα vα −2s′ ψα (r ′ , −s′ , τ ′ )ψα (r, s, τ ).  X (2.36)  α Without proton-neutron mixing, the local pairing density can be written as X X ρ̃q (r) = uα vα (−2s) ψα (r, −s)ψα (r, s), (2.37) α∈q s 19 which can be further simplified when the conjugate partner α is the time-reversed state of α: uα vα |ψα (r, s)|2 . XX ρ̃q (r) = (2.38) α∈q s One should note that pairing cutoffs defined in the canonical and quasiparticle bases are not fully equivalent, and one way to achieve equivalence is to employ a pairing functional that does not require a cutoff (e.g., the momentum-dependent pairing functional proposed in Ref. [77]); but the effect brought by the difference should be minor when the cutoff is high enough. The total HFB Routhian to minimize in the canonical-basis representation is h i R [ψ, ψ ∗ , v] =E ρ [ψ, ψ ∗ , v] , ρ̃ [ψ, ψ ∗ , v] , ρ̃† [ψ, ψ ∗ , v] X (q) X 2 X  (2.39) − ϵF vα − λαβ ⟨ψβ |ψα ⟩ − δαβ , q∈{n,p} α∈q αβ (q) where ϵF is the Fermi energy of neutrons (q = n) or protons (q = p), and λ is the matrix of Lagrange multipliers that guarantee the orthonormality of canonical states. The HFB energy E here is originally formulated as a functional of the density matrix ρ, pairing density matrix ρ̃ and its Hermitian conjugate ρ̃† (see Sec. 2.1), which then becomes a functional of canonical wave functions ψ, their complex conjugates ψ ∗ and occupation amplitudes v through Eqs. (2.35, 2.36). Due to ⟨ψβ |ψα ⟩ = ⟨ψα |ψβ ⟩∗ , the quantities we need to constrain for orthonormality constitute a Hermitian matrix, and thus the matrix λ should also be Hermitian so that the number of independent Lagrange multipliers in it coincides with the number of independent constraints. 20 The variation of R with respect to ψα∗ yields the mean-field equation δR δE X X 0= = − ψ λ β βα = Ĥ α ψα − ψβ λβα , (2.40) δψα∗ δψα∗ β β where ˆ Ĥα = vα2 ĥ + uα vα h̃, (2.41a) 1 λβα = ⟨ψβ |Ĥα + Ĥβ |ψα ⟩. (2.41b) 2 There are three points worth discussing about the mean-field equation. First, δψ δE can be ∗ α reduced to δE δE δρ δE δ ρ̃ ˆ , ∗ = ∗ + ∗ = vα2 ĥψα + uα vα h̃ψ α (2.42) δψα δρ δψα δ ρ̃ δψα which delivers the explicit expressions of the HF Hamiltonian ĥ and the pairing field h̃ ˆ in the same way as Eq. (2.30). Second, the generalized Hamiltonian Ĥα is state-dependent and hence the full matrix λ has to be taken into account to preserve orthonormality. In contrast, the HF calculation only requires Lagrange multipliers for normalization: ĥψα − εα ψα = 0, εα = ⟨ψα |ĥ|ψα ⟩, (2.43) where all the s.p. states experience the same Hamiltonian ĥ, and εα is the s.p. HF energy. Third, the Hermiticity of the matrix λ is enforced by explicit symmetrization in Eq. (2.41b). The variation of R with respect to vα yields the gap equation: vα2   (q ) h i 0= 2vα hαα − ϵF α + uα − h̃αα , (2.44) uα 21 ˆ where hαα = ⟨ψα |ĥ|ψα ⟩ is the canonical s.p. energy, h̃αα = ⟨ψα |h̃|ψ α ⟩ is the state-dependent pairing gap, and qα ∈ {n, p} stands for the isospin of the state α. The gap equation can be solved in a closed form:   v (q ) hαα − ϵF α u  vα    u1 1 = u ∓ rh , (2.45) u  t2 2 (qα ) 2 i  u   α hαα − ϵF + h̃2αα (q) where the Fermi energy ϵF should be adjusted to fulfill the particle-number constraints vα2 = N, vα2 = Z. X X (2.46) α∈n α∈p One can note that only the diagonal elements of the HF Hamiltonian ĥ and the pairing ˆ in the canonical basis enter the gap equation; therefore, no information about the field h̃ non-diagonal elements is needed to determine the occupation amplitudes. It should also be noted that the HF+BCS approach solves a gap equation in the same form of Eq. (2.44), but it uses matrix elements computed in the HF basis instead. The mean-field equation (2.40, 2.41) and gap equation (2.44) together constitute the self-consistent HFB equations in the canonical basis, which can be numerically solved by the gradient-descent algorithm. The same cutoff scheme as discussed at the end of Sec. 2.2.1 can also be employed in the canonical space: The contribution of one canonical state to all the densities is multiplied by a function of the canonical s.p. energy εα ≡ hαα , and the cutoff function can still take the form of Eq. (2.31). 22 2.2.3 Constrained calculations In previous sections, the HFB Routhian is minimized under particle-number constraints. Other quantities, such as multipole moments, angular momenta. and particle-number fluc- tuations, can also be constrained in a similar manner [3]. The HFB energy obtained under these constraints becomes a function of the expectation values of these observables, and this function defines the potential-energy curve (one constraint) or surface (multiple constraints) of the system. The constrained HFB calculation can be employed for many problems. For instance, deformation-constrained calculations can help us locate the global minimum that corresponds to the ground state, while an unconstrained calculation can be easily stuck in a local minimum. The potential-energy surface obtained from constrained calculations is also a useful tool to study large-amplitude collective motions. Nuclear deformations are usually extracted from the multipole moments, which are de- fined in the spherical coordinate system as Z qλµ = dr ρ(r)rλ Yλµ (Ω), (2.47) where ρ(r) can be the neutron, proton or total (isoscalar) particle density, which gives the neutron, proton or total multipole moment, respectively. Components with µ ̸= 0 vanish when the axial symmetry is preserved. The total isoscalar dipole moment (λ = 1) is related to the shift of the center of mass, which must be constrained at zero when the parity symmetry is violated so that the center of mass is fixed at the origin. The quadrupole and octupole moments can also be expressed in the Cartesian coordinate system: Z   Z h  i Q20 = dr ρ(r) 2z 2 − x2 − y2 , Q30 = 2 2 dr ρ(r) z 2z − 3x − 3y 2 . (2.48) 23 The dimensionless deformation parameters are then r ! r ! 16π 3 16π 3 β2 = Q20 / AR02 , β3 = Q30 / AR03 , (2.49) 5 4π 7 4π where R0 = 1.2A1/3 fm is the semi-empirical expression for the nuclear radius. It is more convenient to use dimensionless β values for the comparison of deformations among different nuclei. There are several methods for constraining multipole moments. One approach is adding a quadratic penalty term to the HFB Routhian 2 R′ = R − cλµ qλµ − q̄λµ , (2.50) where q̄λµ is the desired value of the multipole moment and the Lagrange multiplier cλµ should be large enough to push the minimum to the point of qλµ = q̄λµ . Another choice is adding a linear constraint [106] R′ = R − cλµ qλµ − q̄λµ , (2.51)  where the Lagrange multiplier cλµ is adjusted based on the quasiparticle random-phase ap- proximation (QRPA). One can also choose the augmented Lagrangian method that combines the quadratic penalty and linear constraint [107]. The constraint on the angular-momentum expectation value appears in the cranking calculation for the description of nuclear rotation, and the same cranking formulation can be derived through the introduction of a rotating intrinsic frame [3, 33, 35, 108]. Here we assume that the system rotates around the y axis; then a linear constraint term is added to 24 the Routhian D E  R′ = R − ω Jˆy − J¯y , (2.52)   and correspondingly the cranking term −ω Jˆy appears in the HF Hamiltonian ĥ. In the description of rotation, we usually specify the value of the Lagrange multiplier ω instead of the angular momentum J¯y , because ω can be interpreted as the angular velocity of the system. 2.2.4 Numerical solvers There are numerous HF and HFB solvers developed by researchers within the nuclear-DFT framework, and Ref. [75] provides a comprehensive table summarizing some widely employed solvers. In this section I briefly discuss the main features of solvers used in this dissertation. It should be noted that all the solvers discussed below assume no proton-neutron mixing. The following two codes solve the HFB equation in the quasiparticle basis (Sec. 2.2.1) via direct diagonalization. • HFBTHO [109, 110, 111, 112] solves the HFB problem with axial and time-reversal symmetries, and one can select whether the parity symmetry is imposed. The HFB Hamiltonian is constructed in the axially symmetric harmonic-oscillator (HO) or trans- formed HO basis. The HO basis is specified by the number of shells NHO and axial deformation β2 , and Gaussian-quadrature formulas are utilized to calculate integrals in the coordinate space. In HFBTHO, linear constraints (2.51) are added for deformation- constrained calculations. It also provides the kickoff mode: The first 10 iterations are carried out under deformation constraints specified in the input, and then the con- straints are released in the following iterations. The kickoff mode helps the program 25 better locate the global minimum without exploring the whole potential-energy surface. • HFODD [113, 114, 115, 116, 117, 118, 119, 120, 121] expands quasiparticle wave func- tions in the three-dimensional (3D) Cartesian deformed HO basis. The oscillator length and numbers of HO quanta in three directions can be varied independently. There is no symmetry restriction, but one can choose to impose the time-reversal symmetry and point-group symmetries like parity or signature. HFODD can perform the crank- ing calculation for the description of nuclear rotation. As for deformation-constrained calculations, one can choose to add quadratic penalties (2.50) or to use the augmented Lagrange method. The following three programs solve the HFB equation in the canonical basis (Sec. 2.2.2) using the gradient-descent method. • HFBFFT [122] solves the HFB equation in the 3D coordinate representation, and the fast Fourier transform (FFT) is employed for numerical differentiation. It assumes the time-reversal symmetry, but no spatial symmetry is imposed. It is based on Sky3D [81, 123, 124], a solver that can perform static HF+BCS and time-dependent HF calculations. The numerical details of HFBFFT are discussed in Chapter 6. • Sky2D and Sky1D [125] also work in the coordinate space, but Sky2D imposes the axial symmetry while Sky1D the spherical symmetry. Sky2D solves the HFB problem in the cylindrical coordinate system while Sky1D in the spherical coordinate system. As for numerical differentiation, Sky2D employs the FFT technique while Sky1D uses the five-point finite difference formula. When the reflection symmetry is imposed in Sky2D, only grid points with z > 0 are taken into account. It is worth mentioning that Sky2D is based on the HF+BCS solver SkyAx [126]. 26 2.3 Charge-changing finite amplitude method 2.3.1 FAM equations The HFB theory discussed in Sec. 2.2 is applicable for the ground-state calculation. As for low-lying excited states and giant resonances, the QRPA is often the tool of choice within the nuclear-DFT framework [1, 2, 3]. The problem is that solving the QRPA equation through direct construction and diagonalization of the QRPA matrix is too computationally demanding, especially for deformed systems with no spherical symmetry. The FAM [63, 64], however, provides an efficient scheme for the solution of the QRPA without explicitly constructing a huge matrix. The FAM is based on the small amplitude limit of the time-dependent HFB (TDHFB) equation. It is assumed that the nuclear system evolves with time but remains a quasiparticle vacuum (HFB product state). The corresponding quasiparticle operators are thus time- dependent: Xn o Xn o † † ∗ (t)â + V ∗ (t)↠. β̂k (t) = Ulk (t)âl + Vlk (t)âl , β̂k (t) = Ulk l lk l (2.53) l l The time evolution of β̂k (t) under a time-dependent external field F̂ (t) is given by the TDHFB equation: ∂ β̂k (t) h i i = Ĥ(t) + F̂ (t), β̂k (t) , (2.54) ∂t where Ĥ(t) is the TDHFB Hamiltonian       â  1 h(t) − ϵF ∆(t)   â      1 Ĥ(t) = † â â H(t) = †  , (2.55) 2 2 â â    â † ∗ −∆ (t) −h (t) + ϵF ∗ ↠27 where normal ordering is assumed, and the matrix H(t) has the same form as the static HFB matrix given in Eq. (2.29). The matrix elements of H(t) are obtained in the same way as Eq. (2.30), but the mean fields h(t) and ∆(t) are time-dependent in the TDHFB equation because the densities (and underlying solutions U and V ) are time-dependent. In the small amplitude limit, the system is perturbed by a weak external field F̂ (t) with a fixed frequency ω:   1 X  20 † 02   X 11 B̂ , F̂ (t) = η F̂ e−iωt + F̂ † eiωt , F̂ = Fkl Âkl + Fkl kl + Fkl kl (2.56) 2 kl kl † † † † where Âkl ≡ β̂k β̂l , B̂kl ≡ β̂k β̂l , and η is small. Because of the canonical anticommutation † † relations {β̂k , β̂l } = {β̂k , β̂l } = 0, F 20 and F 02 are chosen to be anti-symmetric. The oscillation of F̂ (t) induces the oscillations of the density matrix and pairing tensor: n o ρ(t) = ρ0 + ηδρ(t) = ρ0 + η δρ(ω)e−iωt + δρ† (ω)eiωt , (2.57a) n o κ(t) = κ0 + ηδκ(t) = κ0 + η δκ (ω)e (+) −iωt (−) + δκ (ω)e iωt , (2.57b) where ρ0 and κ0 are, respectively, the density matrix and pairing tensor of the stationary HFB state that the system oscillates around. The Hamiltonian is thus also oscillating around the stationary HFB Hamiltonian Ĥ0 : n o Ĥ(t) = Ĥ0 + δ Ĥ(t) = Ĥ0 + η δ Ĥ(ω)e−iωt + δ Ĥ † (ω)eiωt , (2.58a) 1 X n 20 † 02 (ω) o X 11 (ω)B̂ , δ Ĥ(ω) = δHkl (ω)Âkl + δHkl kl + δHkl kl (2.58b) 2 kl kl where δH 20 and δH 02 are anti-symmetric. Moreover, the time-dependent quasiparticle an- 28 nihilation operator can be decomposed in a similar way: n o X †n o β̂k (t) = β̂k + δ β̂k (t) eiEk t , δ β̂k (t) = η β̂l Xlk (ω)e−iωt + Ylk∗ (ω)eiωt (2.59) l where Ek is the quasiparticle energy of the static HFB Hamiltonian Ĥ0 , and δ β̂k (t) is ex- panded only in terms of quasiparticle creation operators so that the anticommutation relation † {β̂k (t), β̂l (t)} = δkl is satisfied up to the linear order of η. The FAM amplitudes X and Y must be anti-symmetric to fulfill the anticommutation relation {β̂k (t), β̂l (t)} = 0. With all the definitions given above, the TDHFB equation up to the linear order of η is ∂δ β̂k (t) h i h i i = Ek δ β̂k (t) + Ĥ0 , δ β̂k (t) + δ Ĥ(t) + F̂ (t), β̂k , (2.60) ∂t which yields the FAM equations (Ek + El − ω) Xkl (ω)+δHkl 20 (ω) = −F 20 , (E + E + ω) Y (ω)+δH 02 (ω) = −F 02 . (2.61) kl k l kl kl kl It should be noted that F 11 and δH 11 are absent in Eq. (2.61) as they do not contribute in the linear expansion. Thanks to the anti-symmetric property of X and Y , only half of their matrix elements must be solved. Furthermore, the frequency ω is, in general, a complex number. Now the task is to calculate the induced fields δH 20 and δH 02 . One can expand them in terms of amplitudes X and Y , and obtain the linear response equation         20  A B  I 0  X(ω) F   −ω    = − , (2.62) ∗ B A ∗ 0 −I Y (ω) F 02 29   where A B is the QRPA matrix and I is the identity matrix. Eq. (2.62) becomes the B ∗ A∗ standard QRPA equation when the right-hand side is set to zero:        A B  X n  I 0  Xn      = Ωn   , (2.63) ∗ B A ∗ Yn 0 −I Yn   Xn where the eigenvalue Ωn is the QRPA energy and the eigenvector Yn contains QRPA amplitudes (see Ref. [127] for the relation between the FAM amplitudes and QRPA eigen- modes). The dimension of the QRPA matrix, which equals the number of two-quasiparticle excitations, is huge, especially for deformed nuclei; so its explicit evaluation is often numer- ically infeasible. In the FAM, however, the induced fields are directly computed with no explicit construction of the QRPA matrix. Based on Eq. (2.59), we can obtain time-dependent quasiparticle wave functions X  Ulk (t) = Ulk e−iEk t + η ∗ X ∗ eiωt + V ∗ Y e−iωt e−iEk t , Vlj jk lj jk (2.64a) j X  Vlk (t) = Vlk e−iEk t + η ∗ X ∗ eiωt + U ∗ Y e−iωt e−iEk t , Ulj (2.64b) jk lj jk j where U gives the quasiparticle wave functions of the stationary HFB solution that the  V system oscillates around. The density matrix ρ(t) and pairing tensor κ(t) in the s.p. basis are then calculated from Uk (t) and Vk (t) via Eq. (2.27). Based on ρ(t) and κ(t) one can construct the TDHFB Hamiltonian matrix (2.55) in the s.p. basis, which is connected with the induced fields δH 20 and δH 02 through the Bogoliubov transformation (2.25). In this procedure we only keep terms up to the first order of η; the value of η should be small enough for linearity, but large enough to avoid reaching the numerical precision limit. Explicit expressions for 30 the procedure are presented in Ref. [64]. One can also note that the time-reversal symmetry is broken in the FAM, so the Skyrme-EDF terms involving time-odd densities must be considered for the evaluation of the induced fields. The general FAM formulation discussed above is applicable for any external transition operator F̂ , but for a specific type of F̂ one can use its features to simplify the numerical procedure. As for beta-decay calculations, the external field is generated by the charge- changing operator that transforms a neutron to a proton or vice versa, and hence the charge- changing FAM is also called the proton-neutron FAM (PNFAM) [65]. In this work we concentrate on the beta-minus decay, which involves one-body transition operators that convert a neutron to a proton. Correspondingly, only proton and neutron quasiparticle states are connected by such operators, and only matrix elements between proton and neutron states are non-zero in the FAM amplitudes X and Y , induced density δρ, induced pairing tensor δκ, and induced fields δH 02 and δH 20 . Therefore, although there is usually no proton-neutron mixing in the stationary HFB state the system oscillates around, such mixing is induced in the PNFAM. Here we only need to solve the proton-neutron part of FAM amplitudes while the neutron-proton part can be obtained via the anti-symmetric property Xνπ = −Xπν , Yνπ = −Yπν . 2.3.2 FAM response function The relation between the FAM amplitudes and transition matrix elements is revealed by connecting the small-amplitude TDHFB framework with the time-dependent perturbation theory [63]. In the first-order approximation with regard to the external perturbation F̂ (t), 31 the system evolves as Z t ′ e−iΩn t dt′ eiΩn t |Φn ⟩⟨Φn |F̂ (t′ )|Φ0 ⟩ X |Ψ(t)⟩ = |Φ0 ⟩ − i n −∞ ! (2.65) X η⟨Φn |F̂ |Φ0 ⟩ −iωt η ∗ ⟨Φn |F̂ † |Φ0 ⟩ iωt = |Φ0 ⟩ + |Φn ⟩ e − e , n ω − Ωn + iϵ ω + Ωn − iϵ where |Φ0 ⟩ is the ground state while |Φn ⟩ represents the n-th excited state whose excitation energy is Ωn . Here the frequency ω in F̂ (t) is replaced by ω ± iε so that the external field vanishes when t → −∞. The expectation value of F̂ † is then ⟨Ψ(t)|F̂ † |Ψ(t)⟩ = ⟨Φ0 |F̂ † |Φ0 ⟩ + ηS(F̂ ; ω)e−iωt + ... (2.66a) ! X |⟨Φn |F̂ |Φ0 ⟩|2 |⟨Φn |F̂ † |Φ0 ⟩|2 S(F̂ ; ω) = − (2.66b) n ω − Ωn + iϵ ω + Ωn − iϵ where S(F ; ω) is the response function. Taking the limit ε → 0, we have the transition strength distribution dB(F̂ ; ω) X 1 ≡ |⟨Φn |F̂ |Φ0 ⟩|2 δ (ω − Ωn ) = − Im S(F̂ ; ω). (2.67) dω n π For a complex frequency ω = ωr + iΓ (Γ > 0), we have ( ) 1 ΓX |⟨Φn |F̂ |Φ0 ⟩|2 |⟨Φn |F̂ † |Φ0 ⟩|2 − Im S(F̂ ; ωr + iΓ) = − , (2.68) π π n (ωr − Ωn )2 + Γ2 (ωr + Ωn )2 + Γ2 which is the strength distribution smeared with a Lorentzian function of width Γ. One can also evaluate the expectation value ⟨Ψ(t)|F̂ † |Ψ(t)⟩ with the TDHFB state in 32 the small amplitude limit, and find the relation 1 X n 20∗ 02∗ o S(F̂ ; ω) = Fkl Xkl (ω) + Fkl Ykl (ω) 2 kl ! (2.69) X |⟨Φn |F̂ |Φ0 ⟩|2 |⟨Φn |F̂ † |Φ0 ⟩|2 = − , n ω − Ωn ω + Ωn which connects the transition matrix elements with the FAM solution. Based on Eqs. (2.62, 2.63), Ref. [127] shows that the excitation energy Ωn equals the QRPA energy and the state |Φn ⟩ is the QRPA state (including the QRPA ground state |Φ0 ⟩). When the QRPA stability condition is met, all the energies Ωn are real-valued and we have S(ω ∗ ) = S ∗ (ω). It is also possible to construct a response function that includes the interference between two distinct transition operators [65]. Starting from ⟨Ψ(t)|Ĝ† |Ψ(t)⟩, we obtain 1 X n 20∗ o χ(F̂ , Ĝ; ω) = Gkl Xkl (F̂ ; ω) + G02∗ Y kl kl (F̂ ; ω) 2 kl ! (2.70) X ⟨Φn |F̂ |Φ0 ⟩⟨Φ0 |Ĝ† |Φn ⟩ ⟨Φ0 |F̂ |Φn ⟩⟨Φn |Ĝ† |Φ0 ⟩ = − , n ω − Ωn ω + Ωn where G20 and G02 are defined in the same manner as F 20 and F 02 in Eq. (2.56). 2.3.3 Beta-decay rate As shown in Ref. [65], the beta-decay rate can be calculated via the contour integration of PNFAM response functions, which eliminates the necessity to extract individual QRPA modes from the response functions. Here we take the Fermi transition as an example to illustrate the procedure. Equations for the rates of allowed (Fermi and Gamow-Teller) and first-forbidden transitions can all be found in Ref. [65]. The Fermi transition rate is proportional to the sum of individual transition strengths 33 Bi (F̂ ) from the ground state of the parent nucleus |Φ0 ⟩ to all the energetically allowed states |Φi ⟩ of the daughter nucleus, weighted by a phase-space factor: ln 2 X λF̂ = f (En ) Bn (F̂ ), (2.71) κ n where κ = (6147.0 ± 2.4) s, F̂ = τ̂− is the Fermi transition operator (isospin lowering operator) for the β − decay, En is the energy released in the transition, and f (En ) is the phase-space factor derived from the final-state lepton kinematics. The explicit expression of f can be found in Refs. [65, 128], and the released energy En is computed via [58] (n) (p) En = ωmax − Ωn = ϵF − ϵF + ∆Mn−H − Ωn , (2.72) where Ωn is the QRPA energy, and ∆Mn−H = 0.78227 MeV is the mass difference between (n) (q) the neutron and hydrogen atom. Here ϵF and ϵF represent neutron and proton Fermi energies of the HFB ground state of the parent nucleus, respectively. Based on Eq. (2.69), the transition strength can be expressed as the residue of the response function, and one can thus perform a contour integration in the complex plane of ω to evaluate the rate: ln 2 X ln 2 X h i λF̂ = g(Ωn )Bn (F̂ ) = g(Ωn ) Res S(F̂ ), Ωn κ n κ n (2.73) ln 2 X h i ln 2 1 I = Res g̃S(F̂ ), Ωn = dωg̃(ω)S(F̂ ; ω), κ n κ 2πi C where g(ω) = f (ωmax − ω), g̃(ω) is the analytic continuation of g(ω), and the contour C encloses all the poles entering the summation (2.71). Our choice is a circular contour that 34 crosses the real axis at origin and ωmax : ωmax   ω(θ) = 1 + eiθ , θ ∈ [−π, π]. (2.74) 2 Thanks to the relation S(ω ∗ ) = S ∗ (ω), we only need to compute the integrand on the upper or lower semicircle. The remaining problem is how the analytic continuation of the phase-space factor f should be done numerically. We employ a polynomial or rational function to interpolate f on the real axis; the interpolating function is analytic in the complex plane and can thus be used in the integrand. Because of Runge’s phenomenon [129], a better choice is Thiele’s interpolation formula [67, 130]. A rational function in the form of a continued fraction is utilized to interpolate f on a grid between 0 and ωmax : x − x1 finterp (x) = f (x1 ) + x−x2 , (2.75) ρ (x1 , x2 ) + x−x3 ρ2 (x1 ,x2 ,x3 )−f (x1 )+ ρ3 (x1 ,x2 ,x3 ,x4 )−ρ(x1 ,x2 )+··· where x1 , x2 , · · · constitute the grid on which f is evaluated and ρ denotes the reciprocal difference: x1 − x2 ρ (x1 , x2 ) = , f (x1 ) − f (x2 ) x1 − x3 ρ (x1 , x2 , x3 ) = + f (x2 ) , ρ (x1 , x2 ) − ρ (x2 , x3 ) (2.76) ··· x1 − xn ρ (x1 , x2 , . . . , xn ) = + ρ (x2 , x3 . . . , xn−1 ) . ρ (x1 , x2 , . . . , xn−1 ) − ρ (x2 , x3 , . . . , xn ) With this interpolation the integrand is smooth, and we can employ the Gauss-Legendre 35 quadrature rule to calculate the integral. The total beta-decay rate λ is obtained by summing the rates of all the possible decay transitions. Currently we include allowed and first-forbidden transitions in our calculations. Once we have the total rate, the half life can be evaluated via T1/2 = lnλ2 . 2.3.4 Gamow-Teller resonance Apart from the decay rate, the Gamow-Teller resonance (GTR) is another phenomenon worth investigating through the PNFAM formalism. Among all the final states that can be reached via GT transitions from the ground state of the parent nucleus, the GTR takes a large portion of the total GT transition strength, and it can be identified as a strong enhancement (peak) in the total GT transition strength distribution 1 X − S(GT, K; ω) (2.77) π K=0,±1 where K is the angular momentum projection of the GT operator. For a spherical parent nucleus, the response function S is independent of K and we only need to consider one K value. To locate the GTR, one can evaluate S(GT, K; ω) along a horizontal line close to the real axis in the ω plane, i.e., ω = ωr + iΓ where Γ is a small positive constant; then the maximum of the smeared total GT strength (2.77) corresponds to the GTR. A direct scan along the line is computationally expensive, because a dense grid of ωr with a tiny smearing width Γ is required for high resolution but the PNFAM does not converge well when Γ approaches zero. Besides, the range of ωr must be large enough to ensure that the GTR is included in the scan. For efficient determination of the GTR, we turn to a sparser ωr grid and a larger Γ when the 36 GTR is well separated from other states, and interpolate the response function between grid points to locate the maximum. The width Γ here should be smaller than the grid spacing to avoid missing peaks. A good choice is to interpolate 1/S with Eq. (2.75), because it produces a rational function with the degree of the numerator being either the same (odd number of grid points) or one greater than (even number of grid points) the degree of the denominator, and 1/S has the latter property according to Eq. (2.69). Complex-conjugate points ω ∗ = ωr − iΓ are also added to the interpolation to ensure S(ω ∗ ) = S ∗ (ω). The GTR energy presented in experimental works is usually the excitation energy with respect to the ground state of the daughter nucleus, while the ωr value corresponding to the GTR peak is the QRPA energy relative to the ground state of the parent nucleus. As Ref. [131] points out, the QRPA energy needs to be converted to (p) (n) Ex = EQRPA − m(Z + 1, A) + m(Z, A) − ∆Mn−H + ϵF − ϵF , (2.78) for direct comparison with the experimental GTR energy. Here m(Z, A) and m(Z + 1, A) are the atomic masses of parent and daughter nuclei, respectively. 2.3.5 Numerical PNFAM solver The development and test of the numerical PNFAM solver has been discussed in Ref. [65], while its applications can be found in Refs. [62, 66, 67]. In PNFAM, the calculation starts from a HFBTHO solution with the axial symmetry imposed. A Python package called PyNFAM [67] has been developed to manage HFBTHO and PNFAM calculations in a parallel manner with MPI, and to calculate observables like GTRs and decay rates based on PNFAM outputs. The parallelization in PyNFAM is straightforward: First, calculations of different 37 nuclei can always be done in parallel; second, for one nucleus, HFBTHO runs with different deformation kickoffs are performed simultaneously, and PNFAM computations at distinct frequencies ω can also run at the same time. In the PNFAM program, the performance has been significantly boosted (≳ 10× speedup) by replacing explicit loops with BLAS [132] matrix-matrix multiplication routines. The Hamiltonian matrix element in the PNFAM is constructed via the Gaussian-quadrature rule: XZ d3 r ψα∗ (r, s)h(r, s, s′ )ψβ (r, s′ ) ≈ wi ψα∗ (r i , s)h(r i , s, s′ )ψβ (r i , s′ ), (2.79) X hαβ = s,s′ s,s′ ,i where ψ(r, s) is the basis wave function at point r with spin s, h(r, s, s′ ) the mean-field matrix element in the coordinate space, and wi the quadrature weight. We can construct (s) √ (s) √ matrices A(s) and B (s) such that Aαi = wi ψα∗ (r i , s), Bjβ = s′ h(r i , s, s′ ) wi ψβ (r i , s′ ). P Then the matrix element hαβ can be expressed in the form of a matrix-matrix multiplication: (s) (s) (s) (s) sA B , which can be efficiently computed by optimized P P P hαβ = s i Aαi Biβ = BLAS libraries. A similar strategy can be used to accelerate the evaluation of local densities. For instance, the particle density is given by ψα (r i , s)ραβ ψβ∗ (r i , s), XX ρ(r i ) = (2.80) s αβ where ραβ is the density matrix in the HO basis. We immediately see that P ∗ β ραβ ψβ (r i , s) is essentially a matrix-matrix multiplication and can be efficiently evaluated. 38 Chapter 3 Nucleon localization function in rotating nuclei This chapter investigates the nucleon localization function (NLF) in rotating systems. Re- sults of this work has been published in Ref. [133]. This chapter is organized as follows. Sec. 3.1 first discusses the definition and interpretation of the NLF, and then the details of cranked calculations are given in Sec. 3.2. Results are presented in Sec. 3.3, with a short summary given in Sec. 3.4. 3.1 Nucleon localization function 3.1.1 Definition Let µ = x, y, z be the spin-quantization axis. Starting from ρq rsµ , r ′ s′µ , the non-local  density of neutrons (q = n) or protons (q = p) defined in Eq. (2.8), one can construct the particle density, kinetic density and current of a specific spin sµ = ± 21 and isospin q as  1 1 ρqsµ (r) = ρq rsµ , rsµ = ρq (r) + σµ sq (r) · eµ , (3.1a) 2 2 1 1 τqsµ (r) = ∇ · ∇′ ρq rsµ , r ′ sµ r=r′ = τq (r) + σµ T q (r) · eµ , (3.1b)  2 2 i 1 1 ∇′ − ∇ ρq rsµ , r ′ sµ (3.1c)   j qsµ (r) = = j q (r) + σµ Jq (r) · eµ , 2 r=r ′ 2 2 where σµ = 2sµ = ±1, eµ is the unit vector in the µ direction, and ρq , sq , τq , T q , j q and Jq are local densities defined in Eq. (2.9). In a rotationally invariant and spin-unpolarized system, densities defined above are independent of the choice of the quantization axis µ, but 39 in a deformed and rotating nucleus one has to explicitly specify µ. For a HF product state |Ψ⟩, the definition of the NLF starts from the probability of finding two nucleons of a given isospin q and spin sµ at spatial locations r and r ′ : † † Pqsµ (r, r ′ ) = ⟨Ψ|arsµ q a ′ ar′ sµ q arsµ q |Ψ⟩ r sµ q (3.2) = ρq (rsµ , rsµ )ρq (r ′ sµ , r ′ sµ ) − |ρq (rsµ , r ′ sµ )|2 . Because of the Pauli exclusion principle, Pqsµ (r, r) = 0. Given that a nucleon with spin sµ and isospin q is located at position r, the conditional probability of finding another nucleon with the same spin and isospin at position r ′ is Pqsµ (r, r ′ ) Rqsµ (r, r ′ ) = . (3.3) ρq (rsµ , rsµ ) To find the conditional probability of like-spin and like-isospin nucleons in the vicinity of each other, we assume that the second nucleon is located within a small spherical shell with radius δ around r; the conditional probability can then be written as: ′ Rqsµ (r, r + δ) = eδ·∇ Rqsµ (r, r ′ ) . (3.4) r=r ′ Let ⟨· · · ⟩ denote an average over the spherical shell. The first non-zero term in the Taylor D E expansion of Rqsµ (r, r + δ) with respect to the radius δ is [134] 1 2 ′2 1 δ ∇ Rqsµ (r, r ′ ) = Dqsµ (r)δ 2 , (3.5) 6 r=r ′ 3 40 where 2 2 1 ∇ρqsµ j qsµ Dqsµ = τqsµ − − , (3.6) 4 ρqsµ ρqsµ is a localization measure that captures the short-range behavior of Rqsµ (r, r ′ ). The deriva- tion for Eq. (3.6) is based on the density-matrix expansion technique [89, 135]; details can be found in Ref. [136]. Following Ref. [36], we define a dimensionless and normalized NLF:  !2 −1 Dqsµ (r) Cqsµ (r) = 1 + TF (r)  , (3.7) τqs µ where τqsTF (r) = 3 6π 2 2/3 ρ5/3 (r) is the Thomas-Fermi kinetic density. The value of C  µ 5 qsµ can vary in the range of [0, 1]; the smaller the conditional probability of finding two like-spin particles near each other, the more localized they are, and the larger the NLF is. According to Eq. (3.1), densities entering the NLF are composed of both time-even and time-odd terms. In a time-reversal invariant system that is spin-saturated or governed by 2 |∇ρ | spin-independent interactions, we have D = 21 τq − 18 ρqq . We also note that the q,sµ =± 21 Jq (r) term does not vanish even when the time-reversal symmetry is conserved [137]. Hence, the current j qsµ (r) does not vanish in the ground-state configuration of an even-even nucleus unless the system is spin-saturated. Although the contribution of the current j qs to the NLF was ignored in previous works [40, 43, 44], it actually almost vanishes inside the nucleus (see Sec. 3.3.2), so we can safely neglect it when using the NLF as a visualization tool. 3.1.2 Alternative interpretation Besides interpreting the NLF as a measure of the conditional probability of like-spin pairs, there is an alternative interpretation which allows the generalization of NLF’s definition. 41 As discussed in Refs. [37, 138], the localization function can also be interpreted in terms of the Pauli exclusion principle. Let us assume that one isolated fermion of given spin sµ and isospin q is located in some region. The wave function of this particle can be written as √ ψqsµ (r) = ρqsµ eiχ(r) , where χ (r) is a position-dependent phase factor whose gradient is related to the current density as j qsµ = ρqsµ ∇χ. The corresponding s.p. kinetic density is the sum of the last two terms in Dqsµ (3.6): 2 2 s.p. 2 1 ∇ρqsµ j qsµ τqsµ = ∇ψqsµ = + , (3.8) 4 ρqsµ ρqsµ where the first term is the von Weizsäcker kinetic density [139]. Therefore, the localization s.p. measure Dqsµ = τqsµ − τqsµ can be interpreted as the excess of kinetic density due to the Pauli exclusion principle. This interpretation of the NLF is more flexible since it does not involve the notion of the conditional probability; with the new interpretation one can straightforwardly generalize the NLF to the case of point-group symmetries. For the rotating systems to be discussed in following sections, parity (P̂ ), y-signature ˆ (R̂y = e−iπJy ), and y-simplex (R̂y = P̂ R̂y ) symmetries are conserved. To display the effects caused by different s.p. orbits, it is convenient to study the NLF of a given y signature ry or y simplex ry . This can be done by expressing local densities in terms of their symmetry- conserving components. For instance, when y simplex is conserved, we have ρq (r) = ρqσ̆y =+1 (r) + ρqσ̆y =−1 (r), (3.9) where σ̆y ≡ ry /i = ±1. In practice, the component ρqσ̆y is computed via summing up the contributions from s.p. levels with y simplex ry = iσ̆y [113, 140, 141]. The kinetic density 42 τq (r) and current j q (r) can be similarly decomposed. With ρqσ̆y , τqσ̆y and j qσ̆y , the NLF for quantum number σ̆y can be defined as 2 2   2 −1 1 ∇ρqσ̆y j qσ̆y Dqσ̆y (r)  Dqσ̆y = τqσ̆y − − , Cqσ̆y (r) = 1 +  TF   , (3.10)  4 ρqσ̆y ρqσ̆y τqσ̆ (r) y TF (r) = 3 6π 2 2/3 ρ5/3 (r). One should note that densities ρ where τqσ̆ qσ̆y , τqσ̆y and j qσ̆y can  y 5 qσ̆y also be decomposed into time-even and time-odd components in the same manner as Eq. (3.1). 3.2 Cranked calculations 3.2.1 Cranked Hartree-Fock calculation The cranked Hartree-Fock (CHF) method is utilized to study the rotational band of su- perdeformed (SD) 152 Dy. Extreme single-particle behavior is seen in this system [142, 143] and the pairing correlation is absent due to large SD gaps at Z = 66 and N = 86 as well as rapid rotation [144, 145, 146]. The collective rotation of 152 Dy has been investigated in many works, e.g., Refs. [94, 147, 148, 149]. We follow the procedure given in Ref. [147] and perform CHF calculations with HFODD [119]. Single-particle wave functions are expanded in a deformed Cartesian HO basis with frequencies ℏωz = 6.246 MeV and ℏω⊥ = 11.200 MeV along the directions parallel and perpendicular to the symmetry axis, respectively, and the total number of basis states is 1013 with HO quanta not exceeding 15 in each direction. The Skyrme parameterization SkM* [97] with its generic time-odd terms [147, 150] is adopted for the CHF calculations. The main idea of the cranked calculation has already been discussed around Eq. (2.52) 43 in Sec. 2.2.3. Here we follow the notation adopted in Sec. 2.2.3 and assume that the system rotates around the y axis. In our CHF calculations, parity, y-signature, and y-simplex symmetries are preserved while time-reversal and axial symmetries are broken; see Refs. [113, 140, 141] for more discussions. One should note that operator R̂y (R̂y ) is time-odd and thus time-reversed s.p. states belong to opposite signature (simplex) eigenvalues. With conserved parity and y signature, the CHF configuration can be labeled in terms of parity-signature blocks [N+,+i , N+,−i , N−,+i , N−,−i ], where Nπry denotes the number of occupied s.p. orbits with parity π ind y-signature ry . The yrast configuration of SD 152 Dy is [22, 22, 21, 21]n ⊗ [16, 16, 17, 17]p [147]. The relative variation of the quadrupole moment Q20 within this configuration is less than 1% in the frequency range ℏω = 0.2 ∼ 0.5 MeV [142], so we constrain Q20 at 42 b to eliminate its possible influence on the NLF. Single-particle Routhians for the SD yrast band of 152 Dy are shown in Fig. 3.1, where SD shell closures are clearly exhibited. In the figure, rotation-aligned s.p. orbits are marked by thicker lines; they are dramatically impacted by rotation and tend to align their angular momenta along the rotational axis (y axis) as the rotational frequency ω rises. Some of them are intruder orbits denoted by their main HO components with large principal quantum numbers N ; they are the lowest N = 7 neutron and N = 6 proton levels [151, 152]. Other rotation-aligned orbits are labeled by the asymptotic quantum numbers (Nilsson quantum numbers) [N nz Λ]Ω of their dominant HO components. On the other hand, some s.p. states around the Fermi level are weakly affected by rotation and thus known as deformation-aligned [32, 153, 154]; they tend to align their angular momenta along the z axis, the symmetry axis when ω = 0. 44 [3, 0, 1] N shell = 7 are identical. Analytical expressions for the s.p. Routhians and wave functions of the CHO model can be found in a number of publications, e.g., Refs. [33, 155, 156]. The CHO s.p. Routhian is 1 1 1 E = ℏΩ+ (n1 + ) + ℏωy (n2 + ) + ℏΩ− (n3 + ) (3.12) 2 2 2 where v u q 2 u 2 2 − ω2 + 8ω 2 ω⊥2 + ω2  t ω⊥ + ωz2 ω⊥ z z Ω± = + ω2 ± , (3.13) 2 2 and n1 , n2 and n3 become HO quantum numbers in x, y and z directions when ω is zero. It should be noted that the consistency relation between the mean-field ellipsoidal deformation and average density distribution [32, 156] is not adopted here. As an analytical counterpart of the SD doubly-magic 152 Dy, the CHO model should also be closed-shell. Figure 3.2 presents the s.p. Routhians of the CHO system with 60 fermions filling SD supershells up to Nshell ≡ 2(n1 + n2 ) + n3 = 6 [32, 154, 157, 158]. As discussed in Ref. [33] and shown in Fig. 3.2, orbits with no CHO quanta along the rotation axis (n2 = 0) and the largest possible value of the difference (n3 −n1 ) carry most of the angular momentum, and thus are most highly aligned; those are the [0, 0, 7] and [0, 0, 6] Routhians shown in the figure. 3.3 Results and discussions 3.3.1 Time-odd local densities In a rotating system, the collective rotational behavior is well characterized by the current density j [149, 159, 160, 161, 162, 163, 164, 165, 166]. Figure 3.3 displays how the current emerges in the CHO model as the rotational frequency ω increases to 0.2ω0 . As ω increases, a 47 ω = 0.05 ω = 0.1 ω = 0.15 ω = 0.2 -2 ×10 10 ~j 1.2 z (fm) 0 0.6 -10 0 -5 0 5 -5 0 5 -5 0 5 -5 0 5 x (fm) Figure 3.3: Current density j in the x-z (y = 0) plane obtained from the CHO model, as a function of the rotational frequency ω (in units of ω0 ). The magnitude |j| (in fm−4 ) is shown by color and line thickness. flow pattern close to a rigid-body rotation gradually develops inside the system. At ω = 0.2ω0 there is dramatic growth in |j|2 , indicating significant angular-momentum alignment; the system becomes more elongated as well. This effect results from the band crossing marked by the arrow in Fig. 3.2, where the prolate s.p. level [0, 0, 7] with large angular momentum becomes occupied while the oblate [3, 0, 0] level becomes empty. Figure 3.4 presents the current distribution of 152 Dy obtained from CHF calculations, with the rotational frequency ℏω ranging from 0 up to 0.8 MeV (angular momentum Jy from 0 to 90ℏ). The leftmost column in the figure displays the result of the benchmark FAM calculation [167], which corresponds to the limit of ω → 0 (see Sec. 2.3 for details). One can find flow patterns close to the rigid-body rotation in both FAM and CHF results. Since the irrotational flow is mainly attributed to pairing correlations [162, 167], it is expected that no significant irrotational current is produced in the CHF calculations without static pairing. Besides the current j, three other time-odd densities are also present in the expression of the NLF: the spin density s, spin-kinetic density T , and spin-current tensor J. Figure 3.5 shows the distributions of s and T at a number of rotational frequencies, while Fig. 3.6 displays J · ey , the projection that enters the NLF when the y axis is selected as the 48 1.0 FAM ω = 0.2 ω = 0.4 ω = 0.6 ω = 0.8 -3 10 ~j n ×10 9.0 0.80 4.5 0.6 -10 z (fm) 0 -3 10 ~j p ×10 0.4 9.0 0.20 4.5 -10 0 0.0 -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 0.0 0.2 0.4 0.6 0.8 1.0 x (fm) Figure 3.4: Current density j in the x-z (y = 0) plane for neutrons (top) and protons (bottom) in the SD yrast band of 152 Dy, as a function of the rotational frequency ω (in MeV/ℏ). The magnitude |j| (in fm−4 ) is shown by color and line thickness. The FAM result is presented in the leftmost column with a different color range. spin-quantization axis. Both s and T are polarized in the y direction, parallel to the aligned angular momentum. As frequency ω increases, the magnitudes of s and T gradually increase while their directions hardly vary. The distributions of J · ey shown in Fig. 3.6 are mainly present at the nuclear surface, similar to the current j. On the other hand, the spin-current tensor does not vanish at ω = 0 and depends weakly on the rotational frequency. 3.3.2 Simplified nucleon localization function An important consequence of the surface characters of j and J·ey is that they only contribute significantly to the NLF at the surface. This observation should be valid in most cases even if the irrotational flow exists (see examples in Refs. [162, 167]). The same feature is also 2 observed in the distribution of squared density gradient ∇ρqsµ since the particle density is quite flat in the nuclear interior. Therefore, we can define a simplified localization function 49 1.0 ω = 0.0 ω = 0.05 ω = 0.1 ω = 0.15 ω = 0.2 C 0.9 10 0.6 0 0.8 0.3 -10 0 Cτ 0.9 10 0.6 0.6 z (fm) 0 0.3 0.4 -10 0 C − Cτ 0.9 10 0.2 0.6 0 0.3 -10 0.0 0 0.0-5 0 5 0.2 -5 0 5 0.4 -5 0 5 0.6 -5 0 5 0.8 -5 0 5 1.0 x (fm) Figure 3.7: C (top), C τ (middle), and their differences (bottom) in the x-z (y = 0) plane obtained from the CHO model, as functions of the rotational frequency ω (in units of ω0 ). behavior is also observed in the yrast band of 152 Dy: Figure 3.8 shows the comparison between the distributions of C and C τ for neutrons with σ̆y = −1 (y-simplex ry = −i) at rotational frequency ℏω = 0.9 MeV, and the difference C − C τ is only significant at the surface. This difference is less pronounced at lower rotational frequencies. One can see that the simplified NLF eliminates the large value of the NLF at the nuclear surface, and it is thus unnecessary to normalize the NLF as Cqs → Cqs ρqs / max ρqs [44].   Considering the interpretation of Dqsµ as the excess of the kinetic density due to the Pauli 2 2 principle, we are not surprised to see that ∇ρqsµ and j qsµ are non-negligible only at the surface where only a limited number of s.p. orbits are prominent and “localized.” Therefore, we can replace the old NLF with the simplified NLF in most cases, except perhaps some extremely dynamic processes where the current and density gradient can become appreciable 51 C Cτ C − Cτ 10 0.6 z (fm) 0 0.3 -10 0 -5 0 5 -5 0 5 -5 0 5 x (fm) Figure 3.8: C (left), C τ (middle), and their difference (right) in the x-z (y = 0) plane for neutrons with σ̆y = −1 (y-simplex ry = −1) in the SD yrast band of 152 Dy at rotational frequency ℏω = 0.9 MeV. inside the nucleus. In addition, this interpretation also helps the study of the Pauli energy in heavy-ion fusion reactions [51]. 3.3.3 Dependence on the choice of spin-quantization axis τ As mentioned in Sec. 3.1, the NLF Cqsµ (3.7) and its simplified version Cqsµ (3.14) depend on the choice of the spin quantization direction µ. This dependence is visualized in Figs. 3.9 (y = 0 plane) and 3.10 (x = 0 plane) for the rotating 152 Dy at frequency ℏω = 0.5 MeV. It is shown that the NLF depends slightly on the choice of the quantization axis µ as well as τ ≈C the plane selection. One can also notice the relation Cqsµ qsµ inside the nuclear volume, no matter which quantization axis µ and which cross section are selected. 3.3.4 Angular-momentum alignment in the CHO model Figure 3.7 has shown the NLFs of the CHO model in the y = 0 plane at different rotational frequencies; patterns in the x = 0 plane are similar. A pattern similar to that shown in Ref. [43] for the deformed HO model can be clearly observed at ω = 0, but as ω increases the pattern gradually becomes blurred. At ω = 0.2ω0 where the band crossing occurs, the 52 Cx C τx Cy C τy Cz C τz 1.0 10 n ↑ 0.6 0 0.3 0.8 -10 0 10 n ↓ 0.6 0 0.6 0.3 -10 z (fm) 10 p ↑ 0 0.6 0.4 0 0.3 -10 0 10 p ↓ 0.2 0.6 0 0.3 -10 0.0 0 0.0 -5 0 5 -5 0 5 0.2 -50.40 5 -5 00.65 -5 0 0.8 5 -5 0 5 1.0 x (fm) Figure 3.9: Cqsµ (3.7) and Cqsτ µ (3.14) in the x-z (y = 0) plane for three spin quantization directions µ = x, y, z, obtained from CHF calculations for the SD yrast band of 152 Dy at rotational frequency ℏω = 0.5 MeV. The symbols ↑ and ↓ represent sµ = +1 and −1, respectively. NLF changes dramatically: The number of maxima along the z axis increases due to the occupation of the [0, 0, 7] state; the number of maxima in the x direction, on the other hand, decreases since the [3, 0, 0] orbit becomes empty. Thus, the NLF is a good indicator for the shell structure and can easily visualize the effect of band crossing. To better display the evolution of the simplified NLF C τ induced by rotation, in Fig. 3.11 we show the difference ∆C τ (r; ω) ≡ C τ (r; ω) − C τ (r; ω = 0), (3.15) together with density variations ∆τ ≡ τ (r; ω) − τ (r; ω = 0) and τ TF ≡ τ TF (r; ω) − τ TF (r; ω = 0). Figure 3.11 shows a clear correspondence between the peaks of ∆C τ and 53 Cx C τx Cy C τy Cz C τz 1.0 10 n ↑ 0.6 0 0.3 0.8 -10 0 10 n ↓ 0.6 0 0.6 0.3 -10 z (fm) 10 p ↑ 0 0.6 0.4 0 0.3 -10 0 10 p ↓ 0.2 0.6 0 0.3 -10 0.0 0 0.0 -5 0 5 -5 0 5 0.2 -50.40 5 -5 00.65 -5 0 0.8 5 -5 0 5 1.0 y (fm) Figure 3.10: Similar to Fig. 3.9 but shown in the y-z (x = 0) plane. valleys (peaks) of ∆τ (∆τ TF ), which is consistent with Eq. (3.14). This phenomenon indi- cates that ∆τ and ∆τ TF oscillates in antiphase, causing a “constructive interference” when we compute their ratio. One can notice that this out-of-phase oscillation already exists in the patterns of τ and τ TF at ω = 0. Figure 3.12(a) shows τ , τ TF and C τ of the non-rotating HO model along the z axis (x = y = 0), as well as the density profile of the [0, 0, 6] orbit. We see that the valleys (peaks) of τ (τ TF or C τ ) roughly coincide with the maximum points of the s.p. density |ψ006 |2 , while lower s.p. levels provide a smooth background for the total τ (τ TF ). This observation is more evident in the one-dimensional (1D) HO model, as presented in Fig. 3.12(b); in this 1D HO model s.p. orbits with quantum number N ≤ 6 are occupied. The antiphase correspondence is expected because the kinetic density τ is based on the gradient of s.p. wave 54 1.0 ωτ = 0.0 ω = 0.05 ω = 0.1 ω = 0.15 ω = 0.2 10 C ∆C τ 0.2 0.6 0 0.0 0.8 0.3 -10 -0.2 0 -3 ×10 10 τ 0.6 0.06 ∆τ 5.0 z (fm) 0 0.03 0 0.4 -5.0 -10 0 -3 ×10 10 τ TF 0.06 ∆τ TF 8.0 0.2 0 0.03 0 -8.0 -10 0.0 0 0.0 -5 0 5 0.2 -5 00.45 -5 0 0.6 5 -5 0 50.8 -5 0 5 1.0 x (fm) Figure 3.11: C τ (top), τ (in fm−5 , middle) and τ TF (in fm−5 , bottom) in the x-z (y = 0) plane, obtained from the CHO model. The leftmost column shows the reference plots at ω = 0 while the other columns show the rotational dependence relative to the ω = 0 reference as a function of the rotational frequency ω (in units of ω0 ). functions while the τ TF depends on the particle density ρ and involves no derivatives. The oscillating patterns of τ and τ TF basically reflect the characteristic nodal structure of high- N s.p. orbits, allowing us to have a sense of the shell structure by counting the number of maxima; the NLF then successfully magnifies this nodal pattern thanks to the “constructive interference” between τ and τ TF . In addition, the nodal structure is closely associated with clustering, which has been discussed in many works, e.g., Refs. [166, 168, 169, 170]; and hence the NLF can also visualize clusters inside the nucleus. From the s.p. perspective, the density differences shown in Fig. 3.11 are related to the particle-hole (p-h) excitations across the Fermi level induced by the cranking operator ω L̂y , especially the low-energy transitions with ∆n1 = ±1, ∆n2 = 0, ∆n3 = ∓1. Crossings and mixtures of s.p. levels far below the Fermi energy, on the other hand, can barely impact densities. Figure 3.13 displays the kinetic-density variations generated by the p-h excitations 55 1.0 Cωτ = 0.0 ω = 0.2 ω = 0.4 ω = 0.6 ω = 0.8 10 n ↑ 0.6 0 0.3 0.8 -10 0 10 C τn ↓ 0.6 0 0.3 0.6 -10 z (fm) 10 C τp ↑ 0 0.6 0.4 0 0.3 -10 0 10 C τp ↓ 0.2 0.6 0 0.3 -10 0.0 0 0.0-5 0 5 0.2 -5 0 5 0.4 -5 0 5 0.6 -5 0 5 0.8 -5 0 5 1.0 x (fm) τ Figure 3.14: The NLF Cqσ̆ in the x-z (y = 0) plane as a function of ω (in MeV/ℏ), obtained y from CHF calculation for the SD yrast band of 152 Dy. The symbols ↑ and ↓ represent σ̆y = +1 and −1 (y-simplex ry = +i and −i), respectively. time-even and time-odd components. For example, the difference ∆τqσ̆y can be written as 1 1 ∆τqσ̆y (r; ω) = ∆τq (r; ω) + σ̆y Tq′ (r; ω), (3.16) 2 2 where the time-even term ∆τq (r; ω) = τq (r; ω) − τq (r; ω = 0) produces the same background in ∆τqσ̆y =±1 , while the time-odd term Tq′ (r; ω) = τqσ̆y =+1 (r; ω) − τqσ̆y =−1 (r; ω) results in the difference between the NLFs of different simplexes. For the connection between the NLFs and s.p. orbits in the rotating 152 Dy, we follow Sec. 3.2.2 and concentrate on the study of ∆τ . It is difficult to figure out individual p-h excitations 58 ω = 0.2 1.0 ∆Cτ ω = 0.4 ω = 0.6 ω = 0.8 10 n ↑ 0.1 0 0 0.8 -10 -0.1 10 ∆Cτn ↓ 0.1 0 0 0.6 -10 -0.1 z (fm) 10 ∆Cτp ↑ 0.1 0.4 0 0 -10 -0.1 10 ∆Cτp ↓ 0.2 0.1 0 0 -10 -0.1 0.0 0.0 -5 0 0.2 5 -5 00.45 -50.60 5 0.8 -5 0 5 1.0 x (fm) τ . The reference value of C τ at ω = 0 is shown Figure 3.15: Similar to Fig. 3.14 but for ∆Cqσ̆ y in the leftmost column of Fig. 3.14. in the realistic case, so we track individual s.p. levels near the Fermi energy from ℏω = 0 up to 0.1 MeV (see Fig. 3.1) and investigate the evolution of their kinetic densities. When the rotational frequency becomes larger than 0.1 MeV, it is difficult to identify contributions from different s.p. orbits due to strong level mixing. Figure 3.18 presents the kinetic-density variations ∆τ for different parity-signature combinations when ω grows from 0 to 0.1 MeV. These patterns can be approximately reproduced with contributions from a few occupied s.p. levels close to the Fermi energy, as shown in Fig. 3.19. Among neutron orbits with negative parity, rotation-aligned levels 71,2 with high N provide important contributions to 59 1.0 τ ω = 0.0 ω = 0.2 ω = 0.4 ω = 0.6 ω = 0.8 ×10 -3 10 n ↑ ∆τ n ↑ 0.04 1.5 0 0 0.02 0.8 -1.5 -10 0 -3 ×10 10 τn ↓ ∆τ n ↓ 0.04 1.5 0 0 0.6 0.02 -1.5 -10 z (fm) 10 τp ↑ 0 ∆τ p ↑ ×10 -3 0.04 1.5 0.4 0 0 0.02 -1.5 -10 0 -3 ×10 10 τp ↓ 0.2 ∆τ p ↓ 0.04 1.5 0 0 0.02 -1.5 -10 0.0 0 0.0 -5 0 5 0.2 -5 00.45 -5 0 0.6 5 -5 0 50.8 -5 0 5 1.0 x (fm) Figure 3.16: Similar to Fig. 3.15, but for ∆τqσ̆y (in fm−5 ). The reference value of τqσ̆y at ω = 0 is shown in the leftmost column. the vertical oscillation of ∆τn . Among neutron orbits with positive parity, the highest four occupied levels, [651]1/2, [642]5/2, [413]5/2, and [411]1/2, are closely lying and deformation- aligned; they are responsible for the patterns of ∆τn , especially the horizontal structures. For protons, rotation-aligned orbits 61,2,3,4 and [541]1/2 are most crucial and the quantum numbers of their dominant HO components explains their contributions to the difference ∆τp . One can now reach the same conclusion as Sec. 3.3.4: The nodal structures of ∆τ (∆τ TF and ∆C τ ) along the z axis are attributed to rotation-aligned s.p. orbits below the Fermi energy with large N and n3 , while their horizontal features are associated with deformation-aligned levels. 60 ω = 0.0 1.0 TF ω = 0.2 ω = 0.4 ω = 0.6 ω = 0.8 -3 ×10 10 τn ↑ ∆τnTF↑ 3.0 0.04 0 0 0.02 0.8 -3.0 -10 0 -3 ×10 10 τnTF↓ ∆τnTF↓ 3.0 0.04 0 0 0.6 0.02 -3.0 -10 z (fm) 10 τpTF↑ 0 ∆τpTF↑ ×10 -3 0.04 3.0 0.4 0 0 0.02 -3.0 -10 0 -3 ×10 10 τpTF↓ 0.2 ∆τpTF↓ 3.0 0.04 0 0 0.02 -3.0 -10 0.0 0 0.0 -5 0 5 0.2 -5 00.45 -5 0 0.6 5 -5 0 50.8 -5 0 5 1.0 x (fm) TF (in fm−5 ). Figure 3.17: Similar to Fig. 3.16, but for ∆τqσ̆ y 3.4 Summary In this chapter, the nucleon localization function has been studied in anisotropic, spin- unsaturated and spin-polarized rotating systems. The concept of the NLF is first gener- alized to the case of point-group symmetries with the help of the interpretation that the NLF measures the kinetic-density excess owing to the Pauli principle. Then we propose a simplified NLF that is easier to compute and interpret. It is shown that the “constructive interference” between the kinetic density τ and Thomas-Fermi kinetic density τ TF makes the NLF a powerful tool for the visualization of the nodal structure of high-lying s.p. orbits, which explains why the NLF pattern is closely connected with the shell structure and clus- tering. In the rotating 152 Dy, time-odd effects are clearly displayed in the NLFs of opposite 61 Chapter 4 Origin of reflection-asymmetric shapes This chapter focuses on the microscopic origin of ground-state reflection-asymmetric defor- mations of even-even nuclei and investigates this problem from the perspectives of the mul- tipole expansion and s.p. spectrum. Results discussed in the following has been published in Ref. [122], and main conclusions are presented for Ra and Yb isotopes. The multipole expansion of the total energy is discussed in Sec. 4.1, followed by the behaviors of various multipole components shown in Sec. 4.2. In Sec. 4.3 the origin of pear-like deformations are studied through the spectra of canonical s.p. states. 4.1 Multipole expansion of the EDF The multipole expansion of the total energy was first adopted in Ref. [172] and further studied in Ref. [173] to reveal the origin of quadrupole deformations. We generalize the concept of the expansion for the EDF, and employ it for reflection-asymmetric shapes. Since the octupole moment is the lowest nonzero moment related to the reflection-symmetry violation, the pear-like shape is often referred to as “octupole” deformation. In the Skyrme EDF presented in Sec. 2.1, each term inside the integrand can be written in the form of Γ(r)ϱ(r), where ϱ is some density. When the time-reversal symmetry is conserved, ϱ can be ρ, ρ̃, τ , ∇2 ρ, ∇ · J or Jij (i, j = x, y, z). For the term Ctτ ρtt3 (r)τtt3 (r), one can choose Γ = Ctτ ρtt3 and ϱ = τ , and it does not matter whether the coupling C τ is density-dependent. With the axial symmetry imposed, the field Γ and density ϱ can be 64 expanded with the help of spherical harmonics YLM as [13] X ξ(r) = ξ[L] (r)YL,M =0 (Ω), ξ = Γ or ϱ, (4.1) L where ξ[L] (r) = ∗ =0 (Ω). The corresponding energy becomes R dΩ ξ(r)YL,M Z XZ X (Γϱ) E (Γϱ) = dr Γ(r)ϱ(r) = dr Γ[L] ϱ[L] = E[L] , (4.2) L L where dr Γ[L] ϱ[L′ ] (L ̸= L′ ) vanishes thanks to the orthogonality of spherical harmonics. R In the isospin representation, the density that determines Γ is in the same isospin channel of ϱ, and we can calculate the multipole components of isoscalar (t = 0) and isovector (t = 1) energies to compare their contributions to the appearance of reflection-asymmetric shapes. The kinetic ED enters the isoscalar energy as the Skyrme parameterization we use in this chapter adopts the same proton and neutron masses. Like-particle pairing EDF (2.21) with Vp ̸= Vn , and the Coulomb interaction (2.22), however, are neither isoscalar nor isovector because they breaks the isospin symmetry. Without proton-neutron mixing, one can also write the EDF in terms of neutron and proton densities with the help of Eq. (2.9), and the total energy is then decomposed into neutron-neutron (n-n), proton-proton (p-p) and proton-neutron (p-n) components. The Coulomb energy contributes to the p-p energy, while neutron and proton pairing energies are included in n-n and p-p energies, respectively. One can also investigate the multipole components of n-n, p-p and p-n energies to examine their roles in the onset of reflection-asymmetric deformations. 65 shows the octupole deformabilities along the isotopic chain of Yb. One can see that there are no pear-shaped Yb ground states, and the corresponding multipole components given in Fig. 4.3(b) vary slowly as the neutron number N increases. To illustrate the cancellation between various multipolarities, Fig. 4.2(c) displays the cumulative sum X L ∆E[0−L] (β3 = 0.05) = ∆E[L′ ] (β3 = 0.05), (4.4) L′ =0 along the Ra isotopic chain. It is seen that a summation up to L = 7 is necessary to con- verge due to the strong cancellation between monopole and octupole components, especially for isotopes with negative total deformabilities, and that dotriacontapole (L = 5) compo- nents significantly contribute to the appearance of pear-like shapes. This observation is in agreement with previous works that show the important roles of high-order deformations in pear-shaped nuclei [174, 175, 176, 177, 178, 179, 180, 181]. 4.3 Relation to single-particle spectra From the s.p. perspective, the octupole deformation results from the coupling between close- lying opposite-parity orbits (parity doublets) with angular-momentum difference ∆ℓ = ∆j = 3. One can find such pairs of orbits across the Fermi energy when neutron or proton numbers are around Noct = 34, 56, 88 or 134, slightly larger than the magic numbers, where the unique-parity intruder shell (ℓ, j) becomes close to the normal-parity shell (ℓ − 3, j − 3). Refs. [178, 180] noted that the dotriacontapole (and even higher-order) coupling is also related to the ∆ℓ = ∆j = 3 excitations, consistent with the essential role of high-order multipolarities shown in the previous section. One should note that the s.p. argument regarding the appearance of reflection-asymmetric shapes is qualitative and cannot be used 68 to exactly determine which isotope is octupole-deformed. The quadrupole deformation makes the simple criterion on ∆ℓ and ∆j difficult to apply as ℓ and j are no longer good quantum numbers, and pairing correlations disfavor the emergence of deformations. Figure 4.4 displays the canonical s.p. energies of 224 Ra as functions of the quadrupole deformation β2 . HFB calculations for this figure are carried out with the reflection symmetry imposed (β3 = 0). The Fermi energies and equilibrium quadrupole deformations of a series of Ra isotopes are also marked. Around the Fermi energy, we find two pairs of shells, π1i13/2 ↔ π2f7/2 and ν1j15/2 ↔ ν2g9/2 , whose couplings drive the octupole shape. As the neutron number grows from 130 to 136, the quadrupole deformation emerges and breaks the spherical symmetry; then the octupole deformation appears and causes more level repulsion due to the interactions between opposite-parity pairs of states with the same angular-momentum projection Ω. After we pass the neutron number 136, the pear-like shape gradually becomes less favored as more s.p. levels originating from octupole-driving pairs of shells become occupied. For comparison, the s.p. diagram of 176 Yb is given in Fig. 4.5. Compared with 224 Ra, 176 Yb has a larger equilibrium quadrupole deformation but no close-lying opposite-parity pairs of levels with the same Ω value around the Fermi energy, and its ground state is thus reflection symmetric. 4.4 Summary The microscopic origin of the reflection-asymmetric ground-state deformations is investigated in this chapter, with the help of the multipole expansion of the EDF and the spectrum of canonical s.p. states. The monopole and octupole energy components significantly contribute to the appearance of pear-like deformations, but higher-multipolarity interactions are also 69 crucial due to the strong cancellation between ∆E[0] and ∆E[3] . We also see that the isoscalar t=0 or the proton-neutron part ∆E pn part ∆E[3] [3] plays a dominant role in ∆E[3] . From the s.p. perspective, the emergence of reflection-asymmetric shapes is mainly attributed to the coupling between parity doublets with ∆ℓ = ∆j = 3 around the Fermi energy; both octupole and dotriacontapole couplings between such pair of states are strong. 70 2g7/2 Chapter 5 Model calibration for beta-decay studies This chapter discusses the EDF calibration for beta-decay calculations. The physics model we employ in the calibration has been presented in Sec. 2.3, and the χ2 optimization is carried out to determine optimal model parameters. Sections 5.1 and 5.2 discuss model parameters and experimental observables considered in the fit. Then the basic framework of the χ2 optimization is presented in Sec. 5.3, with numerical details discussed in Sec. 5.4. Results are shown in Sec. 5.5 and a summary is given in Sec. 5.6. 5.1 Model parameters For beta-decay studies within the Skyrme-HFB-PNFAM framework, it is necessary to fit parameters that play important roles in the PNFAM but are not constrained in static HFB calculations. In the Skyrme interaction ED (2.16), only isovector (t = 1) terms with t3 = ±1 contribute to the proton-neutron induced fields in the PNFAM. Because of time-reversal symmetry breaking, H(odd) (2.18) that is bilinear in time-odd densities must be taken into account, but corresponding coupling constants are not constrained by the properties of even- even nuclei and thus need to be calibrated. The Skyrme parameterization we utilize for time-even Skyrme couplings and like-particle pairing strengths is UNEDF1-HFB [100], whose parameters have been well calibrated within the χ2 -optimization framework. Following Ref. [66], we preset some time-odd couplings and exclude them from the fit for simplicity. First, we assume that there is no density dependence j in the time-odd couplings we fit. Second, the local gauge invariance [147] can relate Ct and 73 j ∇j Ct∇J to time-even couplings as Ct = −Ctτ , Ct = Ct∇J . Third, we set C1∆s = 0 to avoid finite-size instabilities [182]. Fourth, as for time-odd tensor terms, we only take C1T into account while C1F and C1∇s are set to zero so that the number of free parameters is limited. In the end, the time-odd couplings that we need to calibrate are C1s and C1T . As for the particle-particle channel, we utilize the pairing ED in the isospin representation (2.20), where neutron-neutron, proton-proton and proton-neutron pairing correlations are all included. Both isoscalar and isovector terms in Eq. (2.20) contribute in the PNFAM, but their strengths are determined in different ways. We use the average of neutron-neutron and proton-proton pairing strengths as the isovector pairing strength V1 , i.e., V1 = Vp + Vn /2,  where Vp and Vn are defined in Eq. (2.21) and their values have been fitted together with time-even Skyrme couplings to reproduce selected properties of even-even nuclei. Since Vp and Vn will be the same as V1 when the isospin symmetry is strictly preserved, our choice of V1 partially fulfills the consistency requirement between static HFB and PNFAM calculations. The isoscalar pairing strength V0 , on the other hand, should be calibrated in the same manner as C1s and C1T . Although the isoscalar pairing is related to the Wigner effect in nuclei with N = Z [183], it is not yet reliably constrained in static HFB calculations for these nuclei. Besides EDF parameters, the effective axial-vector coupling constant gA is also included as a model parameter in the fit. It impacts the strengths of Gamow-Teller (GT) and forbidden transitions. For a free neutron its value is gA free ≈ 1.27 [184], but in nuclei it is quenched due to nuclear-medium effects and deficiencies in nuclear models; see Refs. [185, 186] for reviews on the quenching of the axial-vector coupling. It should be noted that the quenching in principle depends on the nucleus and transition type, and heavy systems usually require strong quenching. In this work, however, we adopt a universal gA parameter for simplicity. 74 In previous global beta-decay studies within the Skyrme-HFB-PNFAM framework [66, 67], an empirical value gA = 1.0 is used and the results are acceptable when compared with experimental data. Hence, this empirical value provides a starting point for our fits. It is a good practice to exploit dimensionless parameters with similar variations in the χ2 -optimization routine. The effective axial-vector coupling gA is already dimensionless, with its value constrained within the range of [0, 2]. A natural choice for the dimensionless isoscalar pairing strength is v0 = V0 / |V1 | ≤ 0. As for the time-odd Skyrme couplings, we transform them into dimensionless Landau-Migdal parameters that are defined by the Landau interaction, a residual-interaction form suitable for the studies of nuclear-matter properties and low-lying excitations [187, 188]. The transformation is [94, 189, 190]:   2 F 2 g0′ = N0 2C1s + 2C1T kF 2 + C1 kF , (5.1a) 3 g1′ = −2N0 C1T kF 2 − 2 N C F k2 , (5.1b) 3 0 1 F 1 h′0 = N0 kF 2 CF , 1 (5.1c) 3 2/3 where kF = 3π 2 ρc /2 is the Fermi momentum of nuclear matter at the saturation density ρc . The normalization factor is N0 = 2m∗ kF / π 2 ℏ2 , where  !−1  −1 2 ∂E 1 2 τ m∗ = = + C ρc (5.2) ℏ2 ∂τ00 ρ =ρc m ℏ2 0 00 is the isoscalar effective mass of symmetric nuclear matter [191, 192]. One can notice that g1′ and h′0 are purely determined by the couplings of tensor terms in the Skyrme EDF. In the fit we have h′0 = 0 as we set C1F = 0, and the stability condition of nuclear matter in the spin-isospin channel yields g1′ > −3 when h′0 = 0 [193], providing a constraint on g1′ . On the 75 Table 5.1: Four optimization schemes with different free and fixed parameter sets. Scheme Free parameters Fixed parameters A g0′ , v0 gA = 1, g1′ = 0 B g0′ , v0 , gA g1′ = 0 C g0′ , g1′ , v0 gA = 1 D g0′ , g1′ , v0 , gA None Table 5.2: GTR energies and their experimental errors (in MeV) taken from Refs. [197, 199, 200, 201] for the four nuclei selected in this work. No. Nucleus EGTR Exp error No. Nucleus EGTR Exp error 1 208 Pb 15.6 0.2 3 90 Zr 8.7 – 2 132 Sn 16.3 0.6 4 112 Sn 8.94 0.25 other hand, no constraint is placed on g0′ , but experiments on nuclear spin-isospin responses give a value of g0′ ∼ 1.6 [194, 195, 196, 197]1 , which is a good starting point for our fits. In summary, model parameters considered in the calibration are g0′ , g1′ , v0 and gA . Table 5.1 summarizes the four optimization schemes we use, where g0′ and v0 are always free parameters but g1′ and gA can be either fixed or free. The comparison between these schemes can answer if g1′ and gA can be well constrained and if they should be included in the fits. 5.2 Fit observables As shown in Tables 5.2 and 5.3, two types of observables are employed in the calibration: Gamow-Teller-resonance (GTR) energies EGTR and β − -decay half lives T1/2 . The half lives are given in ascending order, and their logarithms lg T1/2 (base 10) are adopted in the fit as the half lives can vary by several orders of magnitudes. Our data selection is similar to that of Ref. [66], but some data points have been excluded for various reasons. As for the resonance data, only GTRs of doubly-magic (208 Pb and 132 Sn) 1 It should be noted that experimental papers usually adopt the π + ρ + g ′ model with a different normal- ization factor 1/N0′ = 392 MeV · fm [198] , so a transformation of normalization is necessary here. 76 Table 5.3: Beta-minus-decay half lives and their experimental errors (in second) taken from Ref. [202] for the 25 nuclei selected in this work. Half lives are listed in ascending order. No. Nucleus T1/2 Exp error No. Nucleus T1/2 Exp error 5 98 Kr 0.043 0.004 18 166 Gd 4.8 1.0 6 58 Ti 0.058 0.009 19 156 Nd 5.26 0.2 7 102 Sr 0.069 0.006 20 204 Pt 10.3 1.4 8 82 Zn 0.166 0.011 21 74 Zn 95.6 1.2 9 48 Ar 0.475 0.04 22 52 Ti 102 6 10 60 Cr 0.49 0.01 23 180 Yb 144 30 11 126 Cd 0.515 0.017 24 114 Pd 145.2 3.6 12 114 Ru 0.54 0.03 25 242 U 1008 30 13 134 Sn 1.05 0.011 26 134 Te 2508 48 14 152 Ce 1.4 0.2 27 92 Sr 9399.6 61.2 15 78 Zn 1.47 0.15 28 156 Sm 33840 720 16 72 Ni 1.57 0.05 29 200 Pt 45360 1080 17 92 Kr 1.84 0.008 and semi-magic (90 Zr and 112 Sn) systems are considered, while transitional soft systems 76 Ge, 130 Te and 150 Nd are not well described by the mean-field model and thus excluded. The GTR of 48 Ca and spin-dipole resonances (SDRs) of 90 Zr and 208 Pb are also excluded because their experimental spectra do not exhibit clear resonance peaks (see Refs. [203, 204, 205] for the GTR of 48 Ca and Refs. [196, 206] for the SDRs of 90 Zr and 208 Pb). As for beta-decay data, systems with possible octupole ground states (148 Ba, 226 Rn) [31] are ruled out as the reflection symmetry is imposed in our calculations. 5.3 Least-squares fit Reference [68] provides a comprehensive guide and a compilation of examples for the appli- cation of the least-squares fit in nuclear physics. Here we briefly summarize conclusions that are useful for this work. The optimal model parameters are determined through minimizing 77 the weighted sum of squared residuals (errors): nd nd  sk (x) − dk 2  1 1 χ2 (x) 2 X X = εk (x) = , (5.3) nd − nx nd − nx wk k=1 k=1 where x ∈ Rnx is the (column) vector constituted by model parameters, nx = |x| is the number of model parameters, nd is the total number of observables, subscript k denotes the observable index, sk (x) is the model prediction, dk is the experimental value, and wk is known as the weight or adopted error. When there is more than one observable type, the weight is necessary to make the residual εk = [sk (x) − dk ] /wk dimensionless. For an inexact model, there is arbitrariness on the values of weights and one can adjust them to vary the importance of different observables [207]. However, the statistical assumption we adopt for uncertainty estimation can provide a guide for the weight determination. It is assumed that all weighted residuals εk are independent and follow the same normal distribution with expectation 0 and variance σ 2 , and the χ2 value at the optimal point x̂ is approximately the variance σ 2 . To satisfy this assumption, we should choose the weight wk close to the error of model sk , including theoretical, numerical and experimental errors, and then the χ2 (x̂) should be approximately 1. Although each point can have a distinct weight, the points of the same type usually share the same weight as their errors are expected to be close. One can note that what actually matters is the relative weights between observable types, because we can introduce a global scale factor s (Birge factor [208]) such that √ χ2 (x̂) → χ̃2 (x̂) = χ2 (x̂)/s = 1, wk → w̃k = wk s. (5.4) Then, for consistency between weights and residual distributions, the scaled weight for a 78 given observable type w̃typ should be close to nd s [sk (x̂) − dk ]2 , X rtyp = (5.5) ntyp (nd − nx ) k∈typ where ntyp is the number of points of the given type. With the linear expansions of weighted residuals εk (x) around the optimal point x̂, we can transform the nonlinear optimization problem to a linear one and employ the linear- regression framework. Rigorous mathematical discussions can be found in Ref. [209]; in the following we only list important conclusions without proof. Let x∗ be the true parameter vector. Then the difference (x̂−x∗ ) approximately follows a multivariate normal distribution: x̂ − x∗ ∼ N (0, Cov(x̂)). The covariance matrix is h i−1 Cov(x̂) ≈ χ2 (x̂) T J (x̂)J(x̂) , (5.6) ∂ε where the nd × nx Jacobian matrix J(x) is defined by Jkl = ∂xk . There exist several l approximate evaluations for the covariance matrix, but Eq. (5.6) is simpler, numerically cheaper and stabler than other choices [210]. The correlation matrix is then Cov(x̂)kl Rkl = , (5.7) σk σl where σk = Cov(x̂)kk is the standard deviation of x̂k . The (1 − α) confidence region for p x∗ is n o ∗ n x ∗ T −1 ∗ x ∈ R : (x − x̂) Cov (x̂)(x − x̂) ≤ nx Fnx , nd −nx , 1−α , (5.8) where Fnx , nd −nx , 1−α is the (1 − α) quantile of the F distribution with nx and (nd − nx ) 79 degrees of freedom. Suppose x = (xT1 , xT2 )T , where x1 ∈ Rn1 and x2 ∈ Rnx −n1 ; the (1 − α) confidence region for the parameter subset x∗1 is n o x∗1 ∈ Rn1 : (x∗1 − x̂1 )T Cov−1 (x̂1 )(x∗1 − x̂1 ) ≤ n1 Fn1 , nd −nx , 1−α , (5.9) where Cov(x̂1 ) is the upper left n1 × n1 submatrix of Cov(x̂). When n1 = 1 we obtain the n o (1−α) confidence interval for the k-th parameter: x∗k ∈ R : x∗k − x̂k ≤ σk tn −nx , 1−α/2 , d where tn −nx , 1−α/2 is the (1 − α/2) quantile of the t distribution with (nd − nx ) degrees of d freedom. As discussed in Ref. [192], the sensitivity matrix is h i−1 T S = J (x̂)J(x̂) J T (x̂), (5.10) ∂ x̂ whose matrix elements are Skl ≈ ∂εk , representing the variations of optimal parameters l when experimental data are slightly changed. To eliminate the impact of different param- eter scales, we introduce normalized model parameters yk = xk /σk ; correspondingly, the ∂ε normalized Jacobian matrix is Jˇkl (y) = ∂yk = σl Jkl (x), and the normalized sensitivity l matrix is h i−1  ∂ ŷk Škl = Skl /σk = JˇT (ŷ)J(ŷ) ˇ JˇT (ŷ) ≈ . (5.11) kl ∂εl The principal component analysis (PCA) [211, 212, 213] is a useful tool to explore the possibility of reducing the dimension of the parameter space. It has been employed to find the number of effective parameters in nuclear mass models [69]. The PCA provides an orthogonal transformation among model parameters such that the first few new parameters can explain most variation in data. In the PCA, we calculate the singular value decomposition (SVD) of 80 the normalized Jacobian matrix Jˇ at the optimal point ŷ to obtain principal components: ˇ kl V̌lm = sm Ǔkm , s1 ≥ s2 ≥ · · · > sn ≥ 0, X J(ŷ) x (5.12) l where V̌ ∈ Rnx ×nx and Ǔ ∈ Rnd ×nd are orthogonal matrices. The transformation V̌ defines a new set of parameters as zm = l yl V̌lm , and the singular value sm indicates the relevance P of zm in the fit. One can prove that the approximate covariance matrix of the new parameter h i set is Cov(ẑ)kl = V̌ T Cov(ŷ)V̌ ∝ s−2 k δkl . Therefore, ẑk and ẑl (k ̸= l) are uncorrelated, kl and the larger the singular value sk is, the better the parameter zk is constrained (the less soft the χ2 surface is in the corresponding direction). To determine the number of effective model parameters, we define a cumulative fraction: Pm s2k Sm = Pk=1 nx 2 , (5.13) k=1 sk which is the percentage of the variation in data that first m components (columns) of matrices Ǔ and V̌ account for. The number of effective parameters is the minimum m that satisfies Sm > Sth . A typical threshold value is Sth = 0.99, and the corresponding parameter set {z1 , z2 , · · · , zm } explains 99% variation. 5.4 Numerical details The program HFBTHO is utilized to calculate the HFB ground states of parent nuclei whose GTRs or beta-decay half lives are fit observables. To be consistent with the fit of UNEDF1- HFB [100], the number of HO shells is NHO = 20. The ground states of 208 Pb, 90 Zr and 132 Sn are spherical. Other HFB calculations start from a series of kickoff deformations 81 β2 = −0.4 + 0.1k (k = 0, 1, · · · , 8) with the HO-basis deformation taking the same β2 value; then the minimum-energy point corresponds to the ground state. Based on HFB results, the PyNFAM program calls the PNFAM routine to compute β − - transition strength distributions, and then evaluates the GTR energy or beta-decay half life; see Sec. 2.3 for details. In GTR calculations, the resonance is searched on a horizontal line in the complex ω plane with a constant imaginary part (smearing width) Γ = 0.5 MeV, while the grid spacing of the real part is 0.4 MeV. In half-life calculations, 60 sample points are used for the Gauss-Legendre integration on the full circular contour (2.74); PNFAM runs on half of the points are actually needed thanks to the relation S(ω ∗ ) = S ∗ (ω). In addition, a 15-point Chebyshev grid is adopted for the rational interpolation of phase-space factors. The least-squares fit routine we choose is POUNDERS (practical optimization using no derivatives for sums of squares) [214] in PETSc/TAO [215, 216, 217]. It is first employed for the fit of Skyrme parameterization UNEDF0 [192], and shown to be robust and efficient among derivative-free optimizers [218, 219]. In POUNDERS, weighted residuals εi are mod- eled by a quadratic function within a trust region to guide the χ2 minimization. To find the statistical quantities discussed in Sec. 5.3, the Jacobian matrix J is approximately evaluated via the central difference formula: ∂εk ε (x̂ + del ) − εk (x̂ − del ) Jkl (x̂) = ≈ k , (5.14) ∂xl x̂ 2d where el ∈ Rnx denotes a standard unit vector with a 1 in the l-th component and zeros elsewhere, and we choose d = 10−3 in our calculations. A Python interface has been de- veloped to connect the PyNFAM program and POUNDERS routine. This interface is also capable of performing calculations for a given group of parameter vectors, and the results 82 Table 5.4: Optimal parameter values obtained from fits with GTR energies only. Standard deviations are given in brackets. The first letter of the fit indicates the scheme (see Table 5.1). Schemes B and D are absent as gA does not affect the GTR energy. Root-mean-square errors (RMSEs) are also presented for comparison. Fit g0′ g1′ v0 RMSE (MeV) A-GTR 1.60553 (0.025) 0 −1.73430 (0.561) 0.118 C-GTR 1.47434 (0.779) 0.38129 (2.353) −1.70663 (0.699) 0.116 Table 5.5: Similar to Table 5.4 but with beta-decay half-life data only. Optimization schemes C and D with g1′ as a free parameter are absent because of severe numerical instability. Fit g0′ v0 gA RMSE A-β 2.42992 (0.503) −1.06817 (0.169) 1 0.924 B-β −0.07262 (0.212) −1.06339 (0.274) 0.12243 (0.059) 0.718 can then be used as training or test data to build a physics-model emulator for the Bayesian calibration. 5.5 Results and discussions 5.5.1 Optimizations involving one observable type In order to examine the role of different data types and to estimate their weights for a fit with the full data set, we first carry out fits with only one type of data and results are shown in Tables 5.4 (GTR data) and 5.5 (beta-decay half-life data). Optimization schemes B and D are not considered in Table 5.4 since gA does not affect the GTR location, while schemes C and D are not listed in Table 5.5 because in these schemes model parameters are pushed to unphysical regions where severe numerical instability is observed. In Fit A-GTR, we see that g0′ is well constrained by the GTR data but v0 is not, since only the GTR locations of semi-magic 90 Zr and 112 Sn weakly depend on v0 . On the other hand, in Fit A-β the half-life data yield a large uncertainty of g0′ and drive it far from the experimentally accepted value, while the standard deviation of v0 is small. Therefore, when 83 Table 5.6: Optimal parameter values obtained from fits with weights wGTR = 0.3 MeV and wβ = 1 (relative weight γ = 0.3 MeV). Standard deviations are given in brackets, and the χ2 values are also presented for comparison. Fit g0′ g1′ v0 gA χ2 A 1.59560 (0.039) 0 −0.99993 (0.178) 1 25.057 B 1.59184 (0.034) 0 −1.19745 (0.179) 0.50345 (0.143) 19.385 C 1.73245 (0.820) −0.37034 (2.143) −0.99920 (0.183) 1 25.019 D 2.72206 (0.422) −2.54125 (0.781) −1.23511 (0.179) 0.41168 (0.132) 17.788 Table 5.7: Correlation matrix (5.7) obtained from Fit A in Table 5.6. A g0′ v0 g0′ 1.000 v0 −0.135 1.000 Table 5.8: Similar to Table 5.7 but for Fit B given in Table 5.6. Correlations larger than 0.6 in absolute value are marked in italics. B g0′ v0 gA g0′ 1.000 v0 −0.097 1.000 gA −0.004 0.615 1.000 results in Table 5.6, we conclude that g1′ cannot be constrained by our data and should not be included as a free parameter in the fit. Figure 5.3 displays individual weighted residuals εk given by the four fits in Table 5.6, which are consistent with the χ2 values shown in the table. One can notice a general trend that the errors of decay half lives are significantly larger for long-lived nuclei (those with large data-point numbers). These nuclei usually have small beta-decay Q values, and Refs. [66, 67] show that their phase-space factors incur greater theoretical errors. In addition, a small Q value also means that the radius of the integration contour (2.74) is small and a very limited number of residues are enclosed, so one peak that is poorly determined in the calculation can lead to a large error in the half life. Hence, one may consider assigning different weights for short- and long-lived nuclei in future fits. 5.5.4 Number of effective model parameters The correlation matrix (5.7) can give us hints on the number of effective model parameters as two well correlated variables should be combined into one effective parameter. Tables 5.7, 5.8, 5.9 and 5.10 present the correlation matrices of the four fits given in Table 5.6. One can 87 Table 5.9: Correlation matrix (5.7) obtained from Fit C in Table 5.6. C g0′ g1′ v0 g0′ 1.000 g1′ −0.999 1.000 v0 0.057 −0.064 1.000 Table 5.10: Similar to Table 5.9 but for Fit D given in Table 5.6. D g0′ g1′ v0 gA g0′ 1.000 g1′ −0.996 1.000 v0 −0.259 0.256 1.000 gA −0.500 0.507 0.622 1.000 note that g0′ and g1′ are highly correlated in Fits C and D, similar to the case of Fit C-GTR in Table 5.4. Also, v0 and gA are moderately correlated in Fits B and D as both of them depend heavily on the half-life data. As discussed in Sec. 5.3, more concrete determination of effective parameters is provided by the PCA. Figure 5.4 displays squared singular values of normalized Jacobian matrices Jˇ obtained from Eq. (5.12) for the four fits listed in Table 5.6, while Fig. 5.5 shows the corresponding cumulative fractions Sm (5.13). In Fits A and B, the number of effective parameters is the same as that of free parameters. In Fits C and D, the 99% threshold is reached when the first 1 and 2 singular values are included in the cumulative quantity Sm (5.13), respectively. Parameters corresponding to the smallest 2 singular values should thus be removed from the set of effective parameters in Fits C and D. Figure 5.6 shows the squared principal components V̌kl2 corresponding to the first (l = 1) and second (l = 2) largest singular values in the four fits; they represent how the two most effective parameters determined by the PCA are constructed from the normalized parameter set {yk }. In Fit A, g0′ and v0 are equally important and democratically contribute to V̌k1 2 and V̌k22 . In Fit B, the first principal component is constituted by v and g that are 0 A 88 2 2 V̌k1 V̌k2 results of various schemes shows that g1′ , the parameter related to the tensor term, cannot be well constrained by data. The correlation matrices and principal component analysis help us reduce the dimension of the effective parameter space. Based on current results, there are a couple of possible extensions for future works. First, the density dependence of C1s should be taken into account in the fit to investigate its role in GTR and beta-decay calculations. Second, as noted in Ref. [220] the two-body weak current plays a crucial role in the quenching of gA , and Ref. [221] has included it in the PNFAM framework. In this new model, low-energy constants involved in the chiral effective field theory become model parameters to calibrate while gA should take the value of gAfree . Moreover, we consider adding new data or including other types of data (e.g., the GTR strength and beta-decay rate of a given order) to better constrain model parameters. One should notice that in the total decay rate, contributions from various transitions can compensate each other; thus, the total half life may not well constrain some model parameters and individual transitions can provide more information. In addition to the χ2 optimization, we are also working on the Bayesian model calibration with the help of the Kennedy-O’Hagen (KOH) framework [222]. The Bayesian framework avoids the assumptions employed in the χ2 optimization and directly produces the posterior distributions of model parameters through Monte Carlo sampling. An emulator is necessary for such sampling as the physics model we are using is not fast enough to finish sampling in a reasonable time. Our choice is the Gaussian-process (GP) emulator [223], a widely used statistical model for emulation. The Python toolkit discussed in this chapter can also be conveniently utilized to generate data for emulator building. 92 Chapter 6 New HFB solver HFBFFT in three-dimensional coordinate space This chapter describes the numerical implementations and benchmarks of HFBFFT, a new HFB solver for even-even nuclei in the 3D coordinate space. The theoretical framework of the new solver has been discussed in Sec. 2.2.2, and this chapter starts from its numerical implementation (Sec. 6.1). Then the benchmarks of the new solver against well-established HFB solvers are discussed in Sec. 6.2, with a brief summary presented in Sec. 6.3. The material discussed in this chapter has recently been published in Ref. [224]. 6.1 Numerical framework 6.1.1 Numerical realization on a grid HFBFFT is adapted from the HF+BCS solver Sky3D, whose numerical framework has been presented in Refs. [81, 123, 124]. This section briefly summarizes the numerical representation of densities, fields and wave functions in the two solvers, and discusses how we correct the problem of Hermiticity violation. All densities, fields and wave functions in HFBFFT are defined on a 3D Cartesian grid. Grid points along one (x, y or z) direction are equidistant, but one can in principle choose different numbers of points and different grid spacings in different directions. In our calcu- lations we use identical geometries along the three directions. For simplicity, the numerical framework will be explained for the 1D case in the following; one can straightforwardly generalize it to the 3D case. 93   The grid points in the x direction are xν = − Nx2+1 + ν δx, ν = 1, . . . , Nx , where Nx is the (even) number of grid points and δx is the grid spacing. On a grid, the action of a local operator U (x) on a wave function ψ(x) is a simple multiplication, i.e., U (xν )ψ(xν ). The action of the momentum operator requires numerical derivatives defined in the Fourier space. The Fourier technique has been shown to be precise and advantageous for large grids [80]. In the Fourier space, the discrete grid points kn are defined as:  Nx (n − 1)δk, n = 1, . . . ,   2 2π kn = , δk = . (6.1)  Nx Nx δx (n − Nx − 1)δk,  n= + 1, . . . , Nx 2 A coordinate-space wave function ψ(xν ) is connected to its Fourier-space counterpart ψ(ke n) via the discrete Fourier transform and its inverse Nx X ψen = exp (−ikn xν )ψν , (6.2a) ν=1 Nx 1 X ψν = exp (ikn xν )ψen , (6.2b) Nx n=1 where ψν ≡ ψ(xν ), ψen ≡ ψ(k e n ). The discrete Fourier transform and its inverse can be efficiently evaluated through the FFT routine provided by the FFTW3 library [225]. One should note that the use of the Fourier space indicates that the wave function ψ is periodic, i.e., ψ(x + Nx δx) = ψ(x). The long-range Coulomb field is not compatible with periodicity, so we deal with it in the same way as Sky3D, i.e., solve Poisson’s equation for an isolated charged distribution [81, 226, 227]. The appropriate integration method that works with the 94 FFT technique and the periodic boundary condition is the trapezoidal rule Z Nx δx Nx 2 X dx f (x) ≈ f (xν )δx, (6.3) − N2x δx ν=1 where all the points have equal weights. m Algorithm 6.1 presents the numerical evaluation of ddxmψ in the FFT scheme, while algo- h i rithm 6.2 shows the method for the calculation of position-varying differentiation dxd B(x) dψ . dx The position-varying derivative in the HF mean field ĥ comes from the functional term Ctτ ρtt3 τtt3 . Mathematical details about the two algorithms can be found in Ref. [228]. Here we only discuss issues related to the point kNx /2+1 and the product rule. One should note that kNx /2+1 in principle can be ± N2x δk. This arbitrariness does not impact the transforms (6.2) but leads to different ψ (m) when m is odd. We can equally split ψ(k e Nx /2+1 ) between positive and negative momenta, and they will cancel each other in the final result of an odd-order derivative, which is equivalent to setting ψ(ke Nx /2+1 ) = 0 in Algorithm 6.1. On the one hand, this ensures that the derivative of a real-valued function is still real; on the other hand, it means that the second derivative is not equivalent to two consecutive first derivatives in this framework. As noted in Refs. [124, 228], the FFT-based differentiation is not compatible with the product rule. However, in Sky3D the position-varying derivative is still computed via the m Algorithm 6.1: Compute the 1D differentiation ψ (m) ≡ ddxmψ , m = 1, 2, 3, · · · . 1: Compute Fourier transform ψ en = FFT[ψν ] via Eq. (6.2a). (m) 2: If m is odd, set ψeNx /2+1 = 0. Then compute ψen = (ikn )m ψen . (m) (m) 3: Compute inverse Fourier transform ψν = FFT−1 [ψen ] via Eq. (6.2b). 95 h i Algorithm 6.2: Compute the 1D position-varying differentiation χ ≡ d dx B(x) dψ dx . 1: Compute Fourier transform ψen = FFT[ψν ] via Eq. (6.2a). 2: Save ψ e Nx /2+1 → Ψ where Ψ is a temporary variable. e e 3: Set ψ e(1) Nx /2+1 = 0 and compute ψn = ikn ψn (see Algorithm 6.1). e e (1) en(1) ] via Eq. (6.2b). 4: Compute inverse Fourier transform ψν = FFT−1 [ψ (1) 5: Calculate ϕν = Bν ψν with Bν = B(xν ), then compute ϕen = FFT[ϕν ]. 2 e(1) (1)   6: Compute ϕ = ik ϕ and set ϕ = − 1 PNx B Nx δk Ψ. n n n e e N /2+1 x N x ν=1 ν 2 e −1 (1) 7: Compute inverse Fourier transform χν = FFT [ϕen ]. product rule d2 ψ   d dψ dB dψ B(x) = + B 2, (6.4) dx dx dx dx ∂x h i which leads to Hermiticity violation. To restore Hermiticity, one should compute d dx B dψ dx with two consecutive first derivatives, but special treatment must be performed for the point kNx /2+1 to ensure that we return to Algorithm 6.1 with m = 2 when B = 1. As suggested in Ref. [228], we keep the term ψeNx /2+1 in the two first derivatives and then average the results of kNx /2+1 = ± N2x δk for the symmetry in the Fourier space. One can prove that this treatment gives Algorithm 6.2 (Algorithm 3 in [228]), the algorithm implemented in HFBFFT. 6.1.2 Iteration scheme in coordinate space The self-consistent HFB equations (2.40, 2.41, 2.44) in the canonical basis are solved via the damped gradient iteration [76, 78, 79, 80, 81]. The main iteration scheme is presented in Algorithm 6.3. The iteration starts from a number of HF+BCS steps (default is 30 HF+BCS steps), and the HF+BCS calculation is initialized by a 3D HO wave function. The initial HO potential can be spherical, axially deformed or triaxially deformed, but triaxiality is not 96 considered in examples shown this work. Both Sky3D and HFBFFT are parallelized with OpenMP and MPI. The parallelization design for Sky3D can be found in Ref. [123], and it is adapted in HFBFFT to perform scalable HFB calculations. In steps 1 and 4 ∼ 6, each thread / process is responsible for a subset of s.p. wave functions. Matrix operations in steps 7 ∼ 9 and 11 are performed with BLAS [132] + LAPACK [229] (sequential or OpenMP) or ScaLAPACK [230] (MPI or OpenMP/MPI hybrid) libraries. For good scaling in MPI parallelization, the 2D block cyclic distribution is employed to divide matrices among processes in steps 7 ∼ 9 and 11. Steps 2, 3 and 10 are not parallelized, i.e., all the threads / processes do the same work at the same time. 6.1.3 Sub-iteration scheme in the configuration space To accelerate convergence, we develop the sub-iteration scheme (step 11 in Algorithm 6.3) in which damped gradient iterations are performed in the configuration space. In the con- figuration space, the canonical s.p. wave function is expanded in some basis {φn } as XΩ ψα = φn cnα , (6.12) n=1 (0) and we choose an expansion basis such that cnα = δnα when entering the sub-iteration scheme. The corresponding damped gradient iteration is    (new)  δ X X  cnα = Ô cnα −  Hα,nm cmα − cnβ λβα   hnn − h11 +E0 m  β   (6.13)  δ  cnβ λ− X = Ô cnα − βα  ,  hnn − h11 + E0 β 97 Algorithm 6.3: Damped gradient iteration scheme of HFBFFT. 1: Given canonical s.p. states and their occupation fractions {ψα , vα , α = 1, ..., Ω}, compute local densities ρ, τ , J and ρ̃ on the grid. 2: Linear mix new densities with old ones for better convergence: ϱ(new) = (1 − γ)ϱ(old) + γϱψ , ϱ = ρ, τ or ρ̃, (6.5) where subscript ψ denotes the density directly computed from wave functions and γ is the mixture ratio whose default value is 0.2. ˆ 3: Calculate HF mean field ĥ and pairing field h̃ based on densities. ˆ 4: Apply ĥ and h̃ on all the wave functions ψα : ĥψα → Ψα , h̃ψˆ (6.6) α→Ψ e α, 5: Compute canonical s.p. energies and pairing gaps: hαα = ⟨ψα |Ψα ⟩, h̃αα = ⟨ψα |Ψ e α ⟩. (6.7) 6: Evaluate the action of the generalized Hamiltonian Ĥ on ψα and overwrite Ψα : Ĥα ψα = vα2 Ψα + uα vα Ψ e α → Ψα . (6.8) 7: Compute X Ψα − ψβ λβα → Ψα , (6.9) β where the matrix of Lagrange multipliers is taken into account. 8: Apply the damping operator D̂ and orthonormalization Ô (new) n o x0 ψα = Ô ψα − D̂Ψα , D̂ = , (6.10) vα2 (T̂ + E0 ) + 21 uα vα h̃0 h i where x0 and E0 are adjustable numerical parameters, and h̃0 = max h̃n (r), h̃p (r) . Empirical values x0 = 0.45, E0 = 100 MeV are employed. 9: Reevaluate the action of the generalized Hamiltonian and compute the matrix of La- grange multipliers ⟨ψβ |Ĥα |ψα ⟩ + ⟨ψα |Ĥβ |ψβ ⟩∗ λβα = . (6.11) 2 10: With new hαα and h̃αα , compute new occupations amplitudes vα via Eq. (2.45). 11: Perform iterations in the configuration space (see Sec. 6.1.3). 12: If convergence is achieved (see Sec. 6.1.4), exit the iteration; otherwise, return to step 1. 98 where Hα,nm = ⟨φn |Ĥα |φm ⟩, λβα = 12 ∗ Hα,nm + Hβ,nm cmα and P  mn cnβ 1D E 1X D E λ−βα = ψβ Ĥα − Ĥβ ψα = c∗nβ φn Ĥα − Ĥβ φm cmα . (6.14) 2 2 mn Default values of E0 and δ are 10 MeV and 2, respectively. The gap equation (2.44) is solved between two damped gradient steps (6.13), but densities, the HF Hamiltonian ĥ and pairing ˆ are not updated in the sub-iteration scheme. In one word, the program tries potential h̃ to find a unitary transformation among canonical s.p. states and a new set of occupation amplitudes to further minimize the HFB Routhian in the sub-iteration scheme. Let N be the total number of grid points. The time complexity of the FFT algorithm is O(N log N ), and thus the complexity of one coordinate-space iteration is O(ΩN log N ). On the other hand, the complexity of one sub-iteration step is O(Ω2 ). Since Ω ≪ N (several hundreds vs. at least 104 ), one iteration in the configuration space is much less computationally expensive than that on a full 3D grid, so a number of configuration-space iterations can be performed two coordinate-space steps to achieve fast convergence. The best combination of iterations in coordinate and configuration spaces should be determined by numerical experiments; the default choice is to perform 100 iterations in the configuration space between two coordinate-space steps. 6.1.4 Convergence criterion In a gradient-descent framework, a natural choice for the convergence check is the norm of the gradient vector. The α-th component of the gradient vector in the coordinate space is ψβ λβα , whose projection in the configuration space is given by ⟨ψβ |Ĥα ψα − P Ĥα ψα − β ψβ λβα ⟩ = λ− − βα . The matrix elements of λ can be combined into one convergence P β 99 measure: v 1 X u u 1 X − 2 ∆S ≡ t 2 λβα . (6.15) 2 Ωq q∈{n,p} α,β∈q The iteration stops once ∆S is smaller than a specific threshold. This convergence measure, however, does not always work for the pure HF case with no pairing present. Starting from the Hermiticity relation ⟨ψα |ĥ|ψβ ⟩ = ⟨ψβ |ĥ|ψα ⟩∗ , one can show that λ− αβ vanishes when both ψα and ψβ are occupied (vα = vβ = 1) or unoccupied (vα = vβ = 0). Hence, in the pure HF calculation the quantity ∆S measures overlaps of occupied and unoccupied orbits, which reaches zero at the HF solution. Without pairing included, it is allowed that the size of the active s.p. space equals the number of particles (Ωn = N , Ωp = Z) and that all the s.p. levels are occupied. Under this circumstance, ∆S always stays zero and is thus not an appropriate measure of convergence; it instead becomes a measure of Hermiticity violation in our numerical implementations (see Secs. 6.1.1 and 6.2.1). r Therefore, we should adopt the convergence measure α ⟨ψα |ĥ |ψα ⟩ − ⟨ψα |ĥ|ψα ⟩ for 2 2 P the pure HF calculation. 6.1.5 Strategies for pairing Since the HF solution with no pairing is a valid but probably unstable HFB solution, the iteration can be easily locked in a no-pairing solution for a long time; this phenomenon is called pairing breakdown. There exist many strategies to avoid this problem, and here we use a method similar to simulated annealing [231]. The iteration begins with an enhanced pairing strength which is gradually reduced to the physical value as the iteration goes on:   (effective) (physical) max(Nenh − n, 0) Vpair = Vpair ηenh +1 , (6.16) Nenh 100 where n is the iteration number. Default values of ηenh and Nenh are 2 and 400, respectively. It is worth mentioning that the pairing annealing strategy also mitigates the point collapse issue [77]. Here we only discuss the origin of the point collapse in a qualitative way. Assume that there exists a s.p. state ψα with hαα ≫ ϵF and thus vα ≈ 0, uα ≈ 1. Then we have uα vα ≫ vα2 and Ĥα ≈ uα vα h̃, ˆ which means that the generalized Hamiltonian becomes a local operator with no derivative when h̃ ˆ is derived from a density-dependent δ interaction. The corresponding solution is ψα (r) ∝ ∼ δ(r − r min ), where r min is the minimum point of the pairing potential. Then the s.p. state α acquires a huge kinetic energy and thus a huge canonical s.p. energy hαα , which drives the occupation fraction even closer to zero (vα → 0) and forms a positive feedback. This phenomenon is called the point collapse since the corresponding wave function collapses into a spatial point around r min . The enhanced pairing strength at the early stage of iterations ensures that the occupation amplitude vα cannot be easily deadlocked at a small value and thus mitigates the point collapse. Otherwise, it will be difficult to recover correct pairing if some s.p. states are kicked out of the pairing window at an early stage far prior to convergence. As for the pairing cutoff, the soft cutoff scheme discussed in Sec. 2.2.1 is adopted. In HFBFFT the cutoff weight w given by Eq. (2.31) is a function of e = hαα − ϵF and we choose ∆Ecut = 0.1Ecut by default. 6.2 Benchmarks Following the procedure presented in Ref. [232], we have benchmarked our new solver HF- BFFT against several well established solvers, including HFBTHO, Sky1D (spherical systems only) and Sky2D. The Skyrme parameterization SLy4 [98] in the particle-hole channel and the mixed density-dependent δ interaction (ρref = 0.32 fm−3 in Eq. (2.21)) in the particle- 101 Table 6.1: Parameters adopted in different solvers for benchmark calculations: HFBTHO – the number of HO shells NHO and the axial deformation of the HO basis β2 ; Sky1D, Sky2D and HFBFFT – the number of points Ni and grid spacing δi in the direction of i; Common parameter – the pairing cutoff energy Ecut . Solver Parameters HFBTHO NHO = 25, β2 = 0 (spherical systems) or 0.2, Ecut = 60 MeV Sky1D Nr = 141, δr = 0.15 fm, Ecut = 15 MeV Sky2D Nr = Nz = 31, δr = δz = 0.7 fm, Ecut = 15 MeV HFBFFT Nx = Ny = Nz = 48, δx = δy = δz = 0.8 fm, Ecut = 15 MeV particle channel are employed for benchmarks. The nuclei we use for benchmarks are: (i) spherical closed-shell nuclei – 208 Pb, 132 Sn, (ii) spherical superfluid nucleus – 120 Sn, (iii) axially deformed superfluid nuclei – 102 Zr, 110 Zr, and (iv) superfluid nucleus with a superde- formed fission isomer – 240 Pu. All the nuclear systems used for benchmarks have axial and reflection symmetries, but they are also rather different: 110 Zr is a weakly bound system while other nuclei in (i), (ii) and (iii) are well bound; 120 Sn and 102 Zr have no proton pairing while static pairing correlations exist for both protons and neutrons in 110 Zr and 240 Pu. Table 6.1 summarizes some parameters utilized in HFBTHO, Sky1D, Sky2D and HF- BFFT. The reflection symmetry is imposed in both HFBTHO and Sky2D. The numbers of active neutron and proton s.p. states in three coordinate-space solvers (Sky1D, Sky2D and HFBFFT) for nuclei in items (i), (ii) and (iii) are Ωn = 176 and Ωp = 126, respectively; they are large enough as the total energy does not vary significantly (< 10 keV) when we increase (Ωn , Ωp ) to (200, 150). For 240 Pu, however, the active s.p. space should be expanded to accommodate its large number of nucleons: We choose (Ωn , Ωp ) = (300, 200) for the ground state and (Ωn , Ωp ) = (400, 300) for the fission isomer. The box size adopted in HFBFFT is also large enough so that proton and neutron densities are small enough (< 10−7 fm−3 ) at the boundary. 102 Table 6.2: Total energies Etot (in MeV) and ∆S (in MeV) obtained from HFBFFT without and with Hermiticity restoration. Digits that do not coincide are marked in bold. Hermiticity broken Hermiticity restored Etot ∆S Etot ∆S 132 Sn −1103.5429 3.44 × 10−7 −1103.5423 1.60 × 10−15 208 Pb −1635.6817 3.16 × 10−7 −1635.6807 1.20 × 10−15 120 Sn −1018.3310 3.44 × 10−7 −1018.3305 4.01 × 10−7 110 Zr −893.8578 4.59 × 10−7 −893.8574 5.54 × 10−7 102 Zr −859.4696 4.94 × 10−7 −859.4692 3.93 × 10−7 6.2.1 Effect of Hermiticity restoration As mentioned in Sec. 6.1, Algorithm 6.2 for the position-varying derivative avoids the Her- miticity breaking caused by the incompatibility between the product rule and FFT-based differentiation. The effect of this Hermiticity restoration is shown in Table 6.2, where one can compare the results of some closed- and open-shell nuclei before and after the restoration. For closed-shell 132 Sn and 208 Pb the HFB problem is reduced to the pure HF case, and thus we choose (Ωn , Ωp ) = (50, 82) for 132 Sn and (Ωn , Ωp ) = (82, 126) for 208 Pb. The Hermiticity violation is then demonstrated by their non-vanishing ∆S in the third column of Table 6.2, and the magnitudes of their ∆S indicate the order of errors brought by the Hermiticity breaking. After the restoration, the values of ∆S obtained from HF calculations reach almost zero as Hermiticity is well preserved in the new implementation. On the other hand, the impact of the Hermiticity breaking on the total energy is usually a few keV and thus insignificant. Moreover, the convergence measure ∆S of a superfluid system does not vary significantly before and after the restoration, and is also larger than the error resulting from the Hermiticity violation. Hence, the Hermiticity breaking is not a critical issue for static HFB calculations. But it should not be neglected in a time-dependent framework where such errors can accumulate through time steps and become appreciable. 103 Table 6.3: Total energies Etot of 132 Sn and 208 Pb obtained from HFBTHO, HFBFFT, Sky1D and Sky2D. Various energy components are also listed: Subscript of densities denotes contributions from corresponding ED terms in Eq. (2.17), “Coul” the Coulomb energy, and “kin” the kinetic energy. Digits that do not coincide with HFBFFT are marked in bold. HFBTHO HFBFFT Sky1D Sky2D 132 Sn Etot −1103.49 −1103.54 −1103.57 −1103.56 Ekin,n 1637.71 1637.97 1638.01 1638.02 Ekin,p 808.44 808.57 808.59 808.56 Eρρ −4876.26 −4877.02 −4877.04 −4877.07 Eρτ 821.49 821.70 821.73 821.72 Eρ∇2 ρ 248.11 248.23 248.25 248.23 Eρ∇·J −84.40 −84.43 −84.44 −84.43 ECoul 341.42 341.44 341.44 341.43 208 Pb Etot −1635.46 −1635.68 −1635.70 −1635.70 Ekin,n 2528.42 2529.13 2529.16 2529.12 Ekin,p 1336.71 1337.06 1337.07 1337.08 Eρρ −7845.66 −7847.54 −7847.63 −7847.57 Eρτ 1329.79 1330.20 1330.22 1330.20 Eρ∇2 ρ 315.12 315.29 315.29 315.30 Eρ∇·J −96.42 −96.45 −96.45 −96.45 ECoul 796.56 796.63 796.63 796.62 6.2.2 Benchmarks without pairing We start our benchmarks from systems without pairing correlations, i.e., closed-shell nu- clei. Table 6.3 shows the energies of doubly magic 132 Sn and 208 Pb obtained from various solvers; contributions from different functional terms are also given for comparison. These two nuclei are spherical so we can use 1D, 2D and 3D solvers for them. One can see that the results of three coordinate-space solvers agree with each other quite well, and the small differences are mainly attributed to their different grid geometries. The agreement between the results of HFBFFT and HFBTHO is acceptable but not perfect; the main reason is the slow convergence of the total energy (especially the kinetic energy) with respect to the number of HO shells [233, 234]. This is clearly shown by the large discrepancies of Ekin and 104 Eρρ between HFBTHO and HFBFFT. According to Refs. [233, 235], the relation between the total energy Etot obtained from a HO-basis-based solver and the number of HO shells N can be approximated by Etot (L) = E∞ + a0 e−2k∞ L , where L ≡ 2(N + 3/2 + 2)b, b is the p oscillator length of the HO basis, and a0 , k∞ and E∞ should be obtained by fitting. Then E∞ is the energy corresponding to an infinitely large model space. The fit for 208 Pb yields E∞ = −1635.786 MeV, which is closer to the energy given by HFBFFT. 6.2.3 Benchmarks with pairing Due to different numerical representations adopted in different solvers, their discretized quasi- particle continua are also different. In addition, inconsistent cutoff schemes are adopted in these solvers, as we prefer a low cutoff energy and a large smearing factor in grid-based solvers due to the dense continuum in the coordinate-space representation. Thus, pairing strengths Vq (q ∈ {p, n}) are not portable and must be renormalized for benchmarks. Our strategy is to adjust pairing strengths in grid-based solvers to match the spectral pairing gap ∆q [70, 72, 236] given by HFBTHO. The pairing gap is defined in the canonical basis as P 2 h̃αα α∈q wα vα ∆q = P 2 , (6.17) α∈q wα vα where wα is the pairing cutoff weight (2.31). In HFBTHO, the neutron pairing strength Vn is tuned to reproduce the average experimental neutron pairing gap of 120 Sn, which is ∆n = 1.25 MeV, while the proton pairing strength Vp takes the same value as Vn . Another quantity based on which we can tune pairing strengths is the sum of kinetic and pairing energies, Ẽkin = Ekin + Epair , as it is less sensitive to the pairing cutoff energy than its com- ponents [102, 237]. But the pairing gap is preferred because it is connected with experimental 105 Table 6.4: Results of spherical 120 Sn obtained from HFBTHO, HFBFFT, Sky1D and Sky2D, including total energies Etot , some energy components, average pairing gaps ∆ and total root-mean-square (rms) radii rrms (in fm). All energies are in MeV. Pairing strengths Vn and Vp adopted in these calculations are equal and listed in the last row. Digits that do not coincide with HFBFFT are marked in bold and numbers used for the pairing renormalization are in italics. HFBTHO HFBFFT Sky1D Sky2D 120 Sn Etot −1018.77 −1018.78 −1018.92 −1018.74 Ekin,n 1340.51 1339.17 1339.14 1338.72 Ekin,p 830.75 831.25 831.31 831.28 Epair,n −12.48 −9.29 −9.14 −9.02 Epair,p 0.00 0.00 0.00 0.00 Ẽkin,n 1328.03 1329.88 1330.01 1329.70 ∆n 1.25 1.25 1.25 1.25 ∆p 0.00 0.00 0.00 0.00 rrms 4.67 4.67 4.67 4.67 Vn , Vp −284.57 −361.80 −367.30 −372.35 observables and usually yields better total-energy agreement among different solvers. Tables 6.4, 6.5 and 6.6 show the results obtained with different solvers for spherical, deformed and superdeformed systems. As shown in these tables, the values of various ob- servables (total energy, total root-mean-square radius and quadrupole moments) in different columns agree fairly well after the renormalization, although energy components differ sig- nificantly. It is worth noting that the calculation for a deformed nucleus should start from the wave functions of a HO (in grid-based solvers) or Woods-Saxon (in HFBTHO) potential with an appropriate deformation. Besides, the fission isomer of 240 Pu in principle should be determined by locating the local minimum on the potential-energy curve obtained from quadrupole-moment constrained calculations (see Fig. 7 in Ref. [232] for an example). This can be easily carried out in HFBTHO and Sky2D; in HFBFFT, however, the constrained calculation has not been implemented yet and the fission isomer is found by initializing iterations with various HO deformations. 106 Table 6.5: Results of axially deformed 102 Zr and 110 Zr obtained from HFBTHO, HFBFFT and Sky2D. Besides quantities shown in Table 6.4, quadrupole moments Q20 (in fm2 ) are also listed here. Pairing strengths Vn and Vp adopted in these calculations are listed in the last two rows. Digits that do not coincide with HFBFFT are marked in bold and numbers used for the pairing renormalization are in italics. 102 Zr 110 Zr HFBTHO HFBFFT Sky2D HFBTHO HFBFFT Sky2D Etot −859.65 −859.69 −859.67 −893.97 −894.01 −894.01 Ekin,n 1202.02 1200.96 1201.97 1368.08 1367.86 1367.13 Ekin,p 651.25 651.22 651.27 632.03 632.16 632.05 Epair,n −3.39 −2.50 −2.39 −3.18 −2.30 −2.19 Epair,p −1.97 −1.42 −1.38 0.00 0.00 0.00 Ẽkin,n 1198.63 1199.53 1199.58 1364.90 1365.56 1364.94 Ẽkin,p 649.28 649.79 649.89 632.03 632.16 632.05 ∆n 0.69 0.69 0.69 0.64 0.64 0.64 ∆p 0.56 0.56 0.56 0.00 0.00 0.00 rrms 4.58 4.58 4.58 4.73 4.73 4.74 Q20,n 639 639 640 789 791 796 Q20,p 411 411 411 444 445 447 Vn −284.57 −367.00 −378.40 −284.57 −371.00 −384.80 Vp −284.57 −372.00 −384.70 −284.57 −371.00 −384.80 Table 6.6: Similar to Table 6.5, but for the ground state and fission isomer of 240 Pu. 240 Pu ground state fission isomer HFBTHO HFBFFT HFBTHO HFBFFT Sky2D Etot −1802.11 −1802.43 −1797.00 −1797.35 −1797.35 Ekin,n 2938.92 2939.94 2922.56 2923.45 2923.43 Ekin,p 1520.95 1521.46 1525.25 1525.52 1525.33 Epair,n −3.11 −2.30 −3.52 −2.60 −2.48 Epair,p −1.54 −1.22 −2.85 −2.19 −2.07 Ẽkin,n 2935.81 2937.64 2919.03 2920.85 2920.55 Ẽkin,p 1519.40 1520.25 1522.39 1523.33 1523.25 ∆n 0.44 0.44 0.47 0.47 0.47 ∆p 0.33 0.33 0.46 0.46 0.46 rrms 5.93 5.93 6.40 6.40 6.40 Q20,n 1784 1782 5063 5072 5071 Q20,p 1166 1165 3336 3344 3343 Vn −284.57 −360.00 −284.57 −369.00 −384.60 Vp −284.57 −355.00 −284.57 −360.00 −375.80 107 6.3 Summary We have developed a new HFB solver HFBFFT in the 3D coordinate space, using the canonical-basis formalism and damped gradient method. The development of the new pro- gram is based on Sky3D, a well-optimized, highly-parallelized HF+BCS program; the paral- lelization framework of HFBFFT is similar to that of Sky3D. A number of implementations has been done to ensure correct results and quick convergence, including the sub-iteration method, soft pairing cutoff, pairing annealing, and the new algorithm to restore Hermiticity in FFT-based numerical derivatives. We analyze a variety of nuclei with different deforma- tions and s.p. structures to benchmark the new solver against HFBTHO, Sky2D, and (for spherical systems) Sky1D. In the benchmarks, pairing strengths are adjusted according to the pairing gap so that we can compare results obtained with different continuum structures. We expect this benchmark procedure to be useful for the development of other DFT solvers. As a 3D coordinate-space solver, HFBFFT performs better than HO-basis-based solvers for the study of deformed and weakly bound systems. As a solver with no spatial symmetry imposed, we will use HFBFFT to study triaxially deformed or reflection-asymmetric ground states in the future. Besides, there are a number of features we plan to add to HFBFFT to make it more versatile, such as the deformation constraint, the blocking procedure for odd-A and odd-odd systems, and the pairing regularization that removes the dependence of pairing strengths on the cutoff. We are also going to further optimize the performance of HFBFFT for modern supercomputers using GPU architectures. 108 Chapter 7 Conclusions As shown in this dissertation, the nuclear-DFT framework has broad applications to various nuclear physics problems. Research works discussed in this dissertation develop the frame- work, and employ it to study nuclear ground-state properties, collective rotation and beta decays. We also demonstrate that the combination of the nuclear DFT and tools developed in other fields (e.g., the localization function developed in the electronic DFT, statistics, and high-performance computing) can be helpful. Chapter 3 shows the NLF patterns in two rotating systems, the CHO model and SD 152 Dy. These two examples demonstrate the “constructive-interference” mechanism that ex- plains the usefulness of the NLF in the visualization of the internal structure; they also reveal the close connection between NLF patterns and high-lying occupied s.p. orbits. Chapter 4 discusses the origin of reflection-asymmetric ground-state shapes. The roles that different multipole components of the energy play in the onset of pear-like deformations are exam- ined for Ra and Yb isotopic chains, and high-order components are small but important due to the strong cancellation between monopole and octupole parts. Besides, the coupling between close-lying orbits with ∆ℓ = ∆j = 3 across the Fermi energy is shown to be re- sponsible for the reflection asymmetry. Chapter 5 presents the procedure and results of the model calibration for beta-decay calculations. The time-odd Skyrme couplings, isoscalar pairing strength, and effective axial-vector coupling are fitted within the χ2 -optimization framework; their uncertainties and correlations are also obtained in the fit, with the number of effective parameters determined by the PCA. The χ2 optimization paves the way for the 109 Bayesian model calibration, and can be utilized to provide high-quality beta-decay inputs with quantified uncertainties for r-process simulations. Finally, chapter 6 finally discusses the numerical implementation and benchmark of the new HFB solver HFBFFT. This new solver uses the canonical-HFB formalism in the 3D coordinate-space representation, so it is an efficient and reliable tool for the studies of weakly bound and large-deformed nuclei. HFBFFT is highly performant and well parallelized, and its correctness is ensured by careful benchmarks against other HFB solvers. We expect to add more features to HFBFFT and further optimize it for future applications. Generally speaking, there are three main ingredients necessary in the nuclear-DFT re- search, namely, (i) universal and well calibrated EDF parameterizations, (ii) versatile and highly performant solvers, and (iii) useful approaches to extract physics information. Re- search projects presented in this dissertation have significantly contributed to all these three aspects. 110 APPENDIX 111 Appendix List of my contributions 1. Tong Li, Mengzhi Chen, Chunli Zhang, Witold Nazarewicz, and Markus Kortelainen. Nucleon localization function in rotating nuclei. Phys. Rev. C 102, 044305 (2020). • Rederived the localization function with no term omitted. • Performed cranked HF calculations for rotating 152 Dy. • Carried out the analysis and prepared figures regarding the results of rotating 152 Dy and the cranked harmonic-oscillator model. • Wrote the first paper draft. 2. Mengzhi Chen, Tong Li, Jacek Dobaczewski, and Witold Nazarewicz. Microscopic origin of reflection-asymmetric nuclear shapes. Phys. Rev. C 103, 034303 (2021). • Prepared Figs. 11 and 12 to facilitate discussions on s.p. levels. • Checked equations in Sec. III and proofread the draft. 3. Mengzhi Chen, Tong Li, Bastian Schuetrumpf, Paul-Gerhard Reinhard, and Witold Nazarewicz. Three-dimensional Skyrme Hartree-Fock-Bogoliubov solver in coordinate- space representation. Comput. Phys. Commun. 276, 108344 (2022). • Performed a large portion of work regarding the code development. • Wrote the last two paragraphs in Sec. 2.2 and the whole Sec. 3.6. • Checked all the equations and proofread the draft. 112 4. Evan M. Ney, Jonathan Engel, Tong Li, and Nicolas Schunck. Global description of β − decay with the axially deformed Skyrme finite-amplitude method: Extension to odd-mass and odd-odd nuclei. Phys. Rev. C 102, 034326 (2020). • Improved the performance of PNFAM by replacing loops with BLAS routines. 5. Model calibration for beta-decay calculations (unpublished). • Added new codes to PyNFAM to connect the physics model with POUNDERS and to provide data for emulator building. • Produced physics-model outputs for the GP emulator building. • Checked the codes for Bayesian model calibration and corresponding results. 6. Presentations • “Model calibration for beta-decay calculations”, TRIUMF theory seminar, online seminar hosted by TRIUMF, Dec 2021. • “Skyrme EDF parameter calibration for beta-decay calculations”, 2021 NUCLEI collaboration meeting, online, Jun 2021. • “Nucleon localization function in rotating nuclei”, 2021 APS April meeting, online, Apr 2021. • “Nucleon localization function in rotating systems”, 2020 NUCLEI collaboration meeting, online, Jun 2020. • “Hartree-Fock-Bogoliubov solver in natural orbit representation using Fast Fourier Transformation”, joint presentation with Mengzhi Chen, 2018 NUCLEI collabo- ration meeting, University of Tennessee, Knoxville, TN, May 2018. 113 BIBLIOGRAPHY 114 BIBLIOGRAPHY [1] M. Bender, P.-H. Heenen, and P.-G. Reinhard. Self-consistent mean-field models for nuclear structure. Rev. Mod. Phys., 75(1):121–180, January 2003. [2] N. Schunck, editor. Energy Density Functional Methods for Atomic Nuclei. IOP Pub- lishing, Bristol, January 2019. [3] P. Ring and P. Schuck. The Nuclear Many-Body Problem. Springer Science & Business Media, March 2004. [4] L. Coraggio, A. Covello, A. Gargano, N. Itaco, and T. T. S. Kuo. Shell-model cal- culations and realistic effective interactions. Prog. Part. Nucl. Phys., 62(1):135–182, January 2009. [5] E. Caurier, G. Martínez-Pinedo, F. Nowacki, A. Poves, and A. P. Zuker. The shell model as a unified view of nuclear structure. Rev. Mod. Phys., 77(2):427–488, June 2005. [6] P. Möller, J. R. Nix, W. D. Myers, and W. J. Swiatecki. Nuclear Ground-State Masses and Deformations. At. Data Nucl. Data Tables, 59(2):185–381, March 1995. [7] P. Möller, J. R. Nix, and K. L. Kratz. Nuclear properties for astrophysical and radioactive-ion-beam applications. At. Data Nucl. Data Tables, 66(2):131–343, July 1997. [8] K. Pomorski and J. Dudek. Nuclear liquid-drop model and surface-curvature effects. Phys. Rev. C, 67(4):044316, April 2003. [9] A. Bhagwat, X. Viñas, M. Centelles, P. Schuck, and R. Wyss. Microscopic-macroscopic approach for binding energies with the Wigner-Kirkwood method. Phys. Rev. C, 81(4):044321, April 2010. [10] A. Bhagwat, X. Viñas, M. Centelles, P. Schuck, and R. Wyss. Microscopic-macroscopic approach for binding energies with the Wigner-Kirkwood method. II. Deformed nuclei. Phys. Rev. C, 86(4):044316, October 2012. [11] H. Hergert. A Guided Tour of ab initio Nuclear Many-Body Theory. Front. Phys., 8, 2020. [12] D. Vautherin and D. M. Brink. Hartree-Fock calculations with Skyrme’s interaction. I. Spherical nuclei. Phys. Rev. C, 5(3):626–647, March 1972. 115 [13] D. Vautherin. Hartree-Fock calculations with Skyrme’s interaction. II. Axially de- formed nuclei. Phys. Rev. C, 7(1):296–316, January 1973. [14] Y. M. Engel, D. M. Brink, K. Goeke, S. J. Krieger, and D. Vautherin. Time- dependent Hartree-Fock theory with Skyrme’s interaction. Nucl. Phys. A, 249(2):215– 238, September 1975. [15] E. Perlińska, S. G. Rohoziński, J. Dobaczewski, and W. Nazarewicz. Local density approximation for proton-neutron pairing correlations: Formalism. Phys. Rev. C, 69(1):014316, January 2004. [16] D. Gogny. In J. de Boer and H. J. Mang, editors, Proceedings of the International Conference on Nuclear Physics, page 48. North-Holland, 1973. [17] D. Gogny. In G. Ripka and M. Porneuf, editors, Proceedings of the International Conference on Nuclear Self-Consistent Fields, page 333. North-Holland, 1975. [18] J. Dechargé and D. Gogny. Hartree-Fock-Bogolyubov calculations with the d1 effective interaction on spherical nuclei. Phys. Rev. C, 21(4):1568–1593, April 1980. [19] J. D. Walecka. A theory of highly condensed matter. Ann. Phys., 83(2):491–529, April 1974. [20] J. Boguta and A. R. Bodmer. Relativistic calculation of nuclear matter and the nuclear surface. Nucl. Phys. A, 292(3):413–428, December 1977. [21] H. A. Jahn, E. Teller, and F. G. Donnan. Stability of polyatomic molecules in degen- erate electronic states I-Orbital degeneracy. Proc. R. Soc. A: Math. Phys. Eng. Sci., 161(905):220–235, July 1937. [22] H. A. Jahn and W. H. Bragg. Stability of polyatomic molecules in degenerate electronic states II-Spin degeneracy. Proc. R. Soc. A: Math. Phys. Eng. Sci., 164(916):117–131, January 1938. [23] P.-G. Reinhard and E. W. Otten. Transition to deformed shapes as a nuclear Jahn- Teller effect. Nucl. Phys. A, 420(2):173–192, May 1984. [24] W. Nazarewicz. Microscopic origin of nuclear deformations. Nucl. Phys. A, 574(1):27– 49, July 1994. [25] P. A. Butler and W. Nazarewicz. Intrinsic reflection asymmetry in atomic nuclei. Rev. Mod. Phys., 68(2):349–421, April 1996. [26] P. A. Butler. Pear-shaped atomic nuclei. Proc. R. Soc. A: Math. Phys. Eng. Sci., 476(2239):20200202, July 2020. 116 [27] V. Spevak and N. Auerbach. Parity mixing and time reversal violation in nuclei with octupole deformations. Phys. Lett. B, 359(3):254–260, October 1995. [28] V. Spevak, N. Auerbach, and V. V. Flambaum. Enhanced T -odd, P -odd electro- magnetic moments in reflection asymmetric nuclei. Phys. Rev. C, 56(3):1357–1369, September 1997. [29] J. Engel, J. L. Friar, and A. C. Hayes. Nuclear octupole correlations and the en- hancement of atomic time-reversal violation. Phys. Rev. C, 61(3):035502, February 2000. [30] J. Dobaczewski, J. Engel, M. Kortelainen, and P. Becker. Correlating Schiff Moments in the Light Actinides with Octupole Moments. Phys. Rev. Lett., 121(23):232501, December 2018. [31] Y. Cao, S. E. Agbemava, A. V. Afanasjev, W. Nazarewicz, and E. Olsen. Landscape of pear-shaped even-even nuclei. Phys. Rev. C, 102(2):024311, August 2020. [32] A. Bohr and B. R. Mottelson. Nuclear Structure, volume II. W. A. Benjamin, 1975. [33] Z. Szymański. Fast Nuclear Rotation. Clarendon Press, January 1983. [34] W. Nazarewicz. The nuclear collective motion. In J. M. Arias and M. Lozano, editors, An Advanced Course in Modern Nuclear Physics, Lecture Notes in Physics, pages 102–140. Springer, Berlin, Heidelberg, 2001. [35] S. Frauendorf. Spontaneous symmetry breaking in rotating nuclei. Rev. Mod. Phys., 73(2):463–514, June 2001. [36] A. D. Becke and K. E. Edgecombe. A simple measure of electron localization in atomic and molecular systems. J. Chem. Phys., 92(9):5397–5403, May 1990. [37] A. Savin, O. Jepsen, J. Flad, O. K. Andersen, H. Preuss, and H. G. von Schnering. Electron localization in solid-state structures of the elements: The diamond structure. Angew. Chem. Int. Ed. Engl., 31(2):187–188, 1992. [38] A. Savin, R. Nesper, S. Wengert, and T. F. Fässler. ELF: The electron localization function. Angew. Chem. Int. Ed. Engl., 36(17):1808–1832, 1997. [39] T. Burnus, M. A. L. Marques, and E. K. U. Gross. Time-dependent electron localization function. Phys. Rev. A, 71(1):010501, January 2005. [40] P.-G. Reinhard, J. A. Maruhn, A. S. Umar, and V. E. Oberacker. Localization in light nuclei. Phys. Rev. C, 83(3):034312, March 2011. 117 [41] B. Schuetrumpf, C. Zhang, and W. Nazarewicz. Clustering and pasta phases in nuclear density functional theory. In Nuclear Particle Correlations and Cluster Physics, pages 135–153. World Scientific, February 2017. [42] J.-P. Ebran, E. Khan, T. Nikšić, and D. Vretenar. Localization and clustering in atomic nuclei. J. Phys. G: Nucl. Part. Phys., 44(10):103001, August 2017. [43] B. Schuetrumpf and W. Nazarewicz. Cluster formation in precompound nuclei in the time-dependent framework. Phys. Rev. C, 96(6):064608, December 2017. [44] C. L. Zhang, B. Schuetrumpf, and W. Nazarewicz. Nucleon localization and fragment formation in nuclear fission. Phys. Rev. C, 94(6):064323, December 2016. [45] G. Scamps and C. Simenel. Impact of pear-shaped fission fragments on mass- asymmetric fission in actinides. Nature, 564(7736):382–385, December 2018. [46] G. Scamps and C. Simenel. Effect of shell structure on the fission of sub-lead nuclei. Phys. Rev. C, 100(4):041602, October 2019. [47] Z. Matheson, S. A. Giuliani, W. Nazarewicz, J. Sadhukhan, and N. Schunck. Cluster radioactivity of 294 118 Og176 . Phys. Rev. C, 99(4):041304, April 2019. [48] J. Sadhukhan, S. A. Giuliani, Z. Matheson, and W. Nazarewicz. Efficient method for estimation of fission fragment yields of r-process nuclei. Phys. Rev. C, 101(6):065803, June 2020. [49] J. Sadhukhan, S. A. Giuliani, and W. Nazarewicz. Theoretical description of fission yields: Toward a fast and efficient global model. Phys. Rev. C, 105(1):014619, January 2022. [50] P. Jerabek, B. Schuetrumpf, P. Schwerdtfeger, and W. Nazarewicz. Electron and nucleon localization functions of Oganesson: Approaching the Thomas-Fermi limit. Phys. Rev. Lett., 120(5):053001, January 2018. [51] A. S. Umar, C. Simenel, and K. Godbey. Pauli energy contribution to the nucleus- nucleus interaction. Phys. Rev. C, 104(3):034619, September 2021. [52] N. R. Council. Connecting quarks with the cosmos: Eleven science questions for the new century. The National Academies Press, Washington, DC, 2003. [53] E. M. Burbidge, G. R. Burbidge, W. A. Fowler, and F. Hoyle. Synthesis of the Elements in Stars. Rev. Mod. Phys., 29(4):547–650, October 1957. [54] A. G. W. Cameron. Stellar evolution, nuclear astrophysics, and nucleogenesis. Tech- nical report, Canada, June 1957. 118 [55] F. Käppeler, R. Gallino, S. Bisterzo, and W. Aoki. The s process: Nuclear physics, stellar models, and observations. Rev. Mod. Phys., 83(1):157–193, April 2011. [56] M. Arnould, S. Goriely, and K. Takahashi. The r-process of stellar nucleosynthesis: Astrophysics and nuclear physics achievements and mysteries. Phys. Rep., 450(4):97– 213, September 2007. [57] C. J. Horowitz, A. Arcones, B. Côté, I. Dillmann, W. Nazarewicz, I. U. Roed- erer, H. Schatz, A. Aprahamian, D. Atanasov, A. Bauswein, T. C. Beers, J. Bliss, M. Brodeur, J. A. Clark, A. Frebel, F. Foucart, C. J. Hansen, O. Just, A. Kankainen, G. C. McLaughlin, J. M. Kelly, S. N. Liddick, D. M. Lee, J. Lippuner, D. Martin, J. Mendoza-Temis, B. D. Metzger, M. R. Mumpower, G. Perdikakis, J. Pereira, B. W. O’Shea, R. Reifarth, A. M. Rogers, D. M. Siegel, A. Spyrou, R. Surman, X. Tang, T. Uesaka, and M. Wang. r-process nucleosynthesis: connecting rare-isotope beam facilities with the cosmos. J. Phys. G: Nucl. Part. Phys., 46(8):083001, July 2019. [58] J. Engel, M. Bender, J. Dobaczewski, W. Nazarewicz, and R. Surman. β decay rates of r-process waiting-point nuclei in a self-consistent approach. Phys. Rev. C, 60(1):014302, June 1999. [59] M. Mumpower, J. Cass, G. Passucci, R. Surman, and A. Aprahamian. Sensitivity studies for the main r process: β-decay rates. AIP Adv., 4(4):041009, April 2014. [60] M. Mumpower, R. Surman, and A. Aprahamian. Variances in r-process predictions from uncertain nuclear rates. J. Phys.: Conf. Ser., 599:012031, April 2015. [61] M. R. Mumpower, R. Surman, G. C. McLaughlin, and A. Aprahamian. The impact of individual nuclear properties on r-process nucleosynthesis. Prog. Part. Nucl. Phys., 86:86–126, January 2016. [62] T. Shafer, J. Engel, C. Fröhlich, G. C. McLaughlin, M. Mumpower, and R. Surman. β decay of deformed r-process nuclei near A = 80 and A = 160, including odd-A and odd-odd nuclei, with the skyrme finite-amplitude method. Phys. Rev. C, 94(5):055802, November 2016. [63] T. Nakatsukasa, T. Inakura, and K. Yabana. Finite amplitude method for the solution of the random-phase approximation. Phys. Rev. C, 76(2):024318, August 2007. [64] P. Avogadro and T. Nakatsukasa. Finite amplitude method for the quasiparticle random-phase approximation. Phys. Rev. C, 84(1):014314, July 2011. [65] M. T. Mustonen, T. Shafer, Z. Zenginerler, and J. Engel. Finite-amplitude method for charge-changing transitions in axially deformed nuclei. Phys. Rev. C, 90(2):024308, August 2014. 119 [66] M. T. Mustonen and J. Engel. Global description of β − decay in even-even nuclei with the axially-deformed skyrme finite-amplitude method. Phys. Rev. C, 93(1):014304, January 2016. [67] E. M. Ney, J. Engel, T. Li, and N. Schunck. Global description of β − decay with the axially deformed skyrme finite-amplitude method: Extension to odd-mass and odd-odd nuclei. Phys. Rev. C, 102(3):034326, September 2020. [68] J. Dobaczewski, W. Nazarewicz, and P.-G. Reinhard. Error estimates of theoretical models: a guide. J. Phys. G: Nucl. Part. Phys., 41(7):074001, May 2014. [69] V. Kejzlar, L. Neufcourt, W. Nazarewicz, and P.-G. Reinhard. Statistical aspects of nuclear mass models. J. Phys. G: Nucl. Part. Phys., 47(9):094001, July 2020. [70] J. Dobaczewski, H. Flocard, and J. Treiner. Hartree-Fock-Bogolyubov description of nuclei near the neutron-drip line. Nucl. Phys. A, 422(1):103–139, June 1984. [71] S. T. Belyaev, A. V. Smirnov, S. V. Tolokonnikov, and S. A. Fayans. Pairing in nuclei in the coordinate representation. Sov. J. Nucl. Phys., 45(5):783–792, 1987. [72] J. Dobaczewski, W. Nazarewicz, T. R. Werner, J. F. Berger, C. R. Chinn, and J. Dechargé. Mean-field description of ground-state properties of drip-line nuclei: Pair- ing and continuum effects. Phys. Rev. C, 53(6):2809–2840, June 1996. [73] J. Dobaczewski and W. Nazarewicz. Hartree-Fock-Bogoliubov Solution of the Pairing Hamiltonian in Finite Nuclei. In Fifty Years of Nuclear BCS, pages 40–60. World Scientific, May 2012. [74] A. Bulgac and Y. Yu. Renormalization of the Hartree-Fock-Bogoliubov Equations in the Case of a Zero Range Pairing Interaction. Phys. Rev. Lett., 88(4):042504, January 2002. [75] M. Bender, R. Bernard, G. Bertsch, S. Chiba, J. Dobaczewski, N. Dubray, S. A. Giuliani, K. Hagino, D. Lacroix, Z. Li, P. Magierski, J. Maruhn, W. Nazarewicz, J. Pei, S. Péru, N. Pillet, J. Randrup, D. Regnier, P.-G. Reinhard, L. M. Robledo, W. Ryssens, J. Sadhukhan, G. Scamps, N. Schunck, C. Simenel, J. Skalski, I. Stetcu, P. Stevenson, S. Umar, M. Verriere, D. Vretenar, M. Warda, and S. Åberg. Future of nuclear fission theory. J. Phys. G: Nucl. Part. Phys., 47(11):113002, October 2020. [76] P.-G. Reinhard, M. Bender, K. Rutz, and J. A. Maruhn. An HFB scheme in natural orbitals. Z. Phys. A, 358(3):277–278, September 1997. [77] N. Tajima. Canonical-basis solution of the Hartree-Fock-Bogoliubov equation on a three-dimensional Cartesian mesh. Phys. Rev. C, 69(3):034305, March 2004. 120 [78] P.-G. Reinhard and R. Y. Cusson. A comparative study of Hartree-Fock iteration techniques. Nucl. Phys. A, 378(3):418–442, April 1982. [79] C. Bottcher, M. R. Strayer, A. S. Umar, and P.-G. Reinhard. Damped relaxation tech- niques to calculate relativistic bound states. Phys. Rev. A, 40(8):4182–4189, October 1989. [80] V. Blum, G. Lauritsch, J. A. Maruhn, and P.-G. Reinhard. Comparison of coordinate- space techniques in nuclear mean-field calculations. J. Comput. Phys., 100(2):364–376, June 1992. [81] J. A. Maruhn, P.-G. Reinhard, P. D. Stevenson, and A. S. Umar. The TDHF code Sky3D. Comput. Phys. Commun., 185(7):2195–2216, July 2014. [82] P. Hohenberg and W. Kohn. Inhomogeneous Electron Gas. Phys. Rev., 136(3B):B864– B871, November 1964. [83] W. Kohn and L. J. Sham. Self-Consistent Equations Including Exchange and Corre- lation Effects. Phys. Rev., 140(4A):A1133–A1138, November 1965. [84] J. E. Drut, R. J. Furnstahl, and L. Platter. Toward ab initio density functional theory for nuclei. Prog. Part. Nucl. Phys., 64(1):120–168, January 2010. [85] J. Erler, N. Birge, M. Kortelainen, W. Nazarewicz, E. Olsen, A. M. Perhac, and M. Stoitsov. The limits of the nuclear landscape. Nature, 486(7404):509–512, June 2012. [86] J. Erler, P. Klüpfel, and P.-G. Reinhard. Self-consistent nuclear mean-field models: example Skyrme–Hartree–Fock. J. Phys. G: Nucl. Part. Phys., 38(3):033101, January 2011. [87] T. H. R. Skyrme. CVII. The nuclear surface. Philos. Mag. J. Theor. Exp. Appl. Phys., 1(11):1043–1054, November 1956. [88] T. H. R. Skyrme. The effective nuclear potential. Nucl. Phys., 9(4):615–634, January 1958. [89] J. W. Negele and D. Vautherin. Density-Matrix Expansion for an Effective Nuclear Hamiltonian. Phys. Rev. C, 5(5):1472–1493, May 1972. [90] J. W. Negele and D. Vautherin. Density-matrix expansion for an effective nuclear Hamiltonian. II. Phys. Rev. C, 11(3):1031–1041, March 1975. [91] F. Stancu, D. M. Brink, and H. Flocard. The tensor part of Skyrme’s interaction. Phys. Lett. B, 68(2):108–112, May 1977. 121 [92] T. Lesinski, M. Bender, K. Bennaceur, T. Duguet, and J. Meyer. Tensor part of the Skyrme energy density functional: Spherical nuclei. Phys. Rev. C, 76(1):014312, July 2007. [93] M. Bender, K. Bennaceur, T. Duguet, P.-H. Heenen, T. Lesinski, and J. Meyer. Tensor part of the Skyrme energy density functional. II. Deformation properties of magic and semi-magic nuclei. Phys. Rev. C, 80(6):064302, December 2009. [94] V. Hellemans, P.-H. Heenen, and M. Bender. Tensor part of the Skyrme energy density functional. III. Time-odd terms at high spin. Phys. Rev. C, 85(1):014326, January 2012. [95] J. Dobaczewski, W. Nazarewicz, and P.-G. Reinhard. Pairing interaction and self- consistent densities in neutron-rich nuclei. Nucl. Phys. A, 693(1):361–373, October 2001. [96] P. Klüpfel, P.-G. Reinhard, T. J. Bürvenich, and J. A. Maruhn. Variations on a theme by Skyrme: A systematic study of adjustments of model parameters. Phys. Rev. C, 79(3):034310, March 2009. [97] J. Bartel, P. Quentin, M. Brack, C. Guet, and H. B. Håkansson. Towards a better parametrisation of Skyrme-like effective forces: A critical study of the SkM force. Nucl. Phys. A, 386(1):79–100, September 1982. [98] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer. A Skyrme parametrization from subnuclear to neutron star densities Part II. Nuclei far from stabilities. Nucl. Phys. A, 635(1):231–256, May 1998. [99] K. W. Schmid and P.-G. Reinhard. Center-of-mass projection of Skyrme-Hartree-Fock densities. Nucl. Phys. A, 530(2):283–302, August 1991. [100] N. Schunck, J. D. McDonnell, J. Sarich, S. M. Wild, and D. Higdon. Error analysis in nuclear density functional theory. J. Phys. G: Nucl. Part. Phys., 42(3):034024, February 2015. [101] M. Kortelainen, J. McDonnell, W. Nazarewicz, P.-G. Reinhard, J. Sarich, N. Schunck, M. V. Stoitsov, and S. M. Wild. Nuclear energy density optimization: Large deforma- tions. Phys. Rev. C, 85(2):024304, February 2012. [102] P. J. Borycki, J. Dobaczewski, W. Nazarewicz, and M. V. Stoitsov. Pairing renor- malization and regularization within the local density approximation. Phys. Rev. C, 73(4):044319, April 2006. [103] P. Bonche, H. Flocard, P. Heenen, S. J. Krieger, and M. S. Weiss. Self-consistent study of triaxial deformations: Application to the isotopes of Kr, Sr, Zr and Mo. Nucl. Phys. A, 443(1):39–63, September 1985. 122 [104] S. J. Krieger, P. Bonche, H. Flocard, P. Quentin, and M. S. Weiss. An improved pairing interaction for mean field calculations using Skyrme potentials. Nucl. Phys. A, 517(2):275–284, October 1990. [105] C. Bloch and A. Messiah. The canonical form of an antisymmetric tensor and its application to the theory of superconductivity. Nucl. Phys., 39:95–106, December 1962. [106] W. Younes and D. Gogny. Microscopic calculation of 240 Pu scission with a finite-range effective force. Phys. Rev. C, 80(5):054313, November 2009. [107] A. Staszczak, M. Stoitsov, A. Baran, and W. Nazarewicz. Augmented Lagrangian method for constrained nuclear density functional theory. Eur. Phys. J. A, 46(1):85– 90, October 2010. [108] T. Nakatsukasa, K. Matsuyanagi, M. Matsuo, and K. Yabana. Time-dependent density-functional description of nuclear dynamics. Rev. Mod. Phys., 88(4):045004, November 2016. [109] M. V. Stoitsov, J. Dobaczewski, W. Nazarewicz, and P. Ring. Axially deformed solution of the Skyrme–Hartree–Fock–Bogolyubov equations using the transformed harmonic oscillator basis. The program HFBTHO (v1.66p). Comput. Phys. Commun., 167(1):43– 63, April 2005. [110] M. V. Stoitsov, N. Schunck, M. Kortelainen, N. Michel, H. Nam, E. Olsen, J. Sarich, and S. Wild. Axially deformed solution of the Skyrme-Hartree-Fock-Bogoliubov equa- tions using the transformed harmonic oscillator basis (II) hfbtho v2.00d: A new version of the program. Comput. Phys. Commun., 184(6):1592–1604, June 2013. [111] R. N. Perez, N. Schunck, R. D. Lasseri, C. Zhang, and J. Sarich. 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. Comput. Phys. Commun., 220:363–375, November 2017. [112] P. Marević, N. Schunck, E. M. Ney, R. Navarro Pérez, M. Verriere, and J. O’Neal. 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 pro- gram. Computer Physics Communications, 276:108367, July 2022. [113] J. Dobaczewski and J. Dudek. Solution of the Skyrme-Hartree-Fock equations in the Cartesian deformed harmonic oscillator basis I. The method. Comput. Phys. Commun., 102(1):166–182, May 1997. 123 [114] J. Dobaczewski and J. Dudek. Solution of the Skyrme—Hartree—Fock equations in the Cartesian deformed harmonic oscillator basis II. The program HFODD. Comput. Phys. Commun., 102(1):183–209, May 1997. [115] J. Dobaczewski and J. Dudek. Solution of the Skyrme–Hartree–Fock equations in the Cartesian deformed harmonic-oscillator basis. (III) HFODD (v1.75r): a new version of the program. Comput. Phys. Commun., 131(1):164–186, September 2000. [116] J. Dobaczewski and P. Olbratowski. Solution of the Skyrme–Hartree–Fock–Bogolyubov equations in the Cartesian deformed harmonic-oscillator basis. (IV) HFODD (v2.08i): a new version of the program. Comput. Phys. Commun., 158(3):158–191, April 2004. [117] J. Dobaczewski and P. Olbratowski. Solution of the Skyrme–Hartree–Fock–Bogolyubov equations in the Cartesian deformed harmonic-oscillator basis. (V) hfodd(v2.08k). Comput. Phys. Commun., 167(3):214–216, May 2005. [118] J. Dobaczewski, W. Satuła, B. G. Carlsson, J. Engel, P. Olbratowski, P. Powałowski, M. Sadziak, J. Sarich, N. Schunck, A. Staszczak, M. Stoitsov, M. Zalewski, and H. Zduńczuk. Solution of the Skyrme–Hartree–Fock–Bogolyubov equations in the Cartesian deformed harmonic-oscillator basis.: (VI) hfodd (v2.40h): A new version of the program. Comput. Phys. Commun., 180(11):2361–2391, November 2009. [119] N. Schunck, J. Dobaczewski, J. McDonnell, W. Satuła, J. A. Sheikh, A. Staszczak, M. Stoitsov, and P. Toivanen. Solution of the Skyrme–Hartree–Fock–Bogolyubov equa- tions in the Cartesian deformed harmonic-oscillator basis.: (VII) hfodd (v2.49t): A new version of the program. Comput. Phys. Commun., 183(1):166–192, January 2012. [120] N. Schunck, J. Dobaczewski, W. Satuła, P. Bączyk, J. Dudek, Y. Gao, M. Konieczka, K. Sato, Y. Shi, X. B. Wang, and T. R. Werner. Solution of the Skyrme- Hartree–Fock–Bogolyubov equations in the Cartesian deformed harmonic-oscillator ba- sis. (VIII) hfodd (v2.73y): A new version of the program. Comput. Phys. Commun., 216:145–174, July 2017. [121] J. Dobaczewski, P. Bączyk, P. Becker, M. Bender, K. Bennaceur, J. Bonnard, Y. Gao, A. Idini, M. Konieczka, M. Kortelainen, L. Próchniak, A. M. Romero, W. Satu\la, Y. Shi, T. R. Werner, and L. F. Yu. Solution of universal nonrelativistic nuclear DFT equations in the Cartesian deformed harmonic-oscillator basis. (IX) HFODD (v3.06h): a new version of the program. J. Phys. G: Nucl. Part. Phys., 48(10):102001, September 2021. [122] M. Chen, T. Li, J. Dobaczewski, and W. Nazarewicz. Microscopic origin of reflection- asymmetric nuclear shapes. Phys. Rev. C, 103(3):034303, March 2021. [123] M. Afibuzzaman, B. Schuetrumpf, and H. M. Aktulga. Scalable nuclear density func- tional theory with Sky3D. Comput. Phys. Commun., 223:34–44, February 2018. 124 [124] B. Schuetrumpf, P.-G. Reinhard, P. D. Stevenson, A. S. Umar, and J. A. Maruhn. The TDHF code Sky3D version 1.1. Comput. Phys. Commun., 229:211–213, August 2018. [125] P.-G. Reinhard, HFB solvers Sky1D and Sky2D, unpublished, 2021. [126] P.-G. Reinhard, B. Schuetrumpf, and J. A. Maruhn. The Axial Hartree–Fock + BCS Code SkyAx. Comput. Phys. Commun., 258:107603, January 2021. [127] N. Hinohara, M. Kortelainen, and W. Nazarewicz. Low-energy collective modes of deformed superfluid nuclei within the finite-amplitude method. Phys. Rev. C, 87(6):064309, June 2013. [128] H. Behrens and W. Bühring. Electron Radial Wave Functions and Nuclear Betadecay. Clarendon Press, 1982. [129] L. N. Trefethen. Approximation Theory and Approximation Practice. SIAM, January 2013. [130] M. Abramowitz and I. A. Stegun. Handbook of Mathematical Functions: With Formu- las, Graphs, and Mathematical Tables. Courier Corporation, January 1965. [131] E. Yüksel, N. Paar, G. Colò, E. Khan, and Y. F. Niu. Gamow-Teller excitations at finite temperature: Competition between pairing and temperature effects. Phys. Rev. C, 101(4):044305, April 2020. [132] L. S. Blackford, J. Demmel, J. Dongarra, I. Duff, S. Hammarling, G. Henry, M. Heroux, L. Kaufman, A. Lumsdaine, A. Petitet, R. Pozo, K. Remington, and R. C. Whaley. An updated set of basic linear algebra subprograms (BLAS). ACM Trans. Math. Softw., 28(2):135–151, June 2002. [133] T. Li, M. Z. Chen, C. L. Zhang, W. Nazarewicz, and M. Kortelainen. Nucleon local- ization function in rotating nuclei. Phys. Rev. C, 102(4):044305, October 2020. [134] A. D. Becke. Hartree–Fock exchange energy of an inhomogeneous electron gas. Int. J. Quantum Chem., 23(6):1915–1922, 1983. [135] J. Dobaczewski, B. G. Carlsson, and M. Kortelainen. The Negele–Vautherin density- matrix expansion applied to the Gogny force. J. Phys. G: Nucl. Part. Phys., 37(7):075106, May 2010. [136] C. Zhang. Interplay between single-particle and collective motion within nuclear density functional theory. PhD thesis, Michigan State University, United States – Michigan, 2016. 125 [137] S. G. Rohoziński, J. Dobaczewski, and W. Nazarewicz. Self-consistent symmetries in the proton-neutron Hartree-Fock-Bogoliubov approach. Phys. Rev. C, 81(1):014313, January 2010. [138] P. Fuentealba, E. Chamorro, and J. C. Santos. Chapter 5 Understanding and using the electron localization function. In A. Toro-Labbé, editor, Theoretical and Com- putational Chemistry, volume 19 of Theoretical Aspects of Chemical Reactivity, pages 57–85. Elsevier, January 2007. [139] C. F. v. Weizsäcker. Zur Theorie der Kernmassen. Z. Physik, 96(7):431–458, July 1935. [140] A. L. Goodman. Self-consistent symmetries of the Hartree-Fock-Bogoliubov equations in a rotating frame. Nucl. Phys. A, 230(3):466–476, September 1974. [141] J. Dobaczewski, J. Dudek, S. G. Rohoziński, and T. R. Werner. Point symmetries in the Hartree-Fock approach. I. Densities, shapes, and currents. Phys. Rev. C, 62(1):014310, June 2000. [142] W. Satuła, J. Dobaczewski, J. Dudek, and W. Nazarewicz. Additivity of quadrupole moments in superdeformed bands: Single-particle motion at extreme conditions. Phys. Rev. Lett., 77(26):5182–5185, December 1996. [143] S. T. Clark, G. Hackman, R. V. F. Janssens, R. M. Clark, P. Fallon, S. N. Floor, G. J. Lane, A. O. Macchiavelli, J. Norris, S. J. Sanders, and C. E. Svensson. Empirical Investigation of Extreme Single-Particle Behavior of Nuclear Quadrupole Moments in Highly Collective A ∼ 150 Superdeformed Bands. Phys. Rev. Lett., 87(17):172503, October 2001. [144] Y. R. Shimizu, E. Vigezzi, and R. A. Broglia. Role of static and dynamic pairing cor- relations in the superdeformed band of 152 Dy. Phys. Lett. B, 198(1):33–36, November 1987. [145] C. Baktash, B. Haas, and W. Nazarewicz. Identical Bands in Deformed and Superde- formed Nuclei. Annu. Rev. Nucl. Part. Sci., 45:485–541, 1995. [146] J. Dobaczewski, W. Nazarewicz, and T. R. Werner. Closed shells at drip-line nuclei. Phys. Scr., T56:15–22, January 1995. [147] J. Dobaczewski and J. Dudek. Time-odd components in the mean field of rotating superdeformed nuclei. Phys. Rev. C, 52(4):1827–1839, October 1995. [148] P. Bonche, H. Flocard, and P.-H. Heenen. Microscopic study of superdeformation in the A = 150 mass region. Nucl. Phys. A, 598(2):169–186, February 1996. 126 [149] A. V. Afanasjev and H. Abusara. Time-odd mean fields in covariant density functional theory: Rotating systems. Phys. Rev. C, 82(3):034329, September 2010. [150] M. Bender, J. Dobaczewski, J. Engel, and W. Nazarewicz. Gamow-Teller strength and the spin-isospin coupling constants of the Skyrme energy functional. Phys. Rev. C, 65(5):054322, May 2002. [151] T. Bengtsson, I. Ragnarsson, and S. Åberg. The role of high-N orbits in superdeformed states. Phys. Lett. B, 208(1):39–44, July 1988. [152] W. Nazarewicz, R. Wyss, and A. Johnson. Structure of superdeformed bands in the A ≈ 150 mass region. Nucl. Phys. A, 503(2):285–330, October 1989. [153] F. S. Stephens. Coriolis effects and rotation alignment in nuclei. Rev. Mod. Phys., 47(1):43–65, January 1975. [154] S. G. Nilsson and I. Ragnarsson. Shapes and Shells in Nuclear Structure. Cambridge University Press, Cambridge, 1995. [155] D. Glas, U. Mosel, and P. G. Zint. The cranked harmonic oscillator in coordinate space. Z Physik A, 285(1):83–87, March 1978. [156] T. Troudet and R. Arvieu. Shapes of nuclear configurations in a cranked harmonic oscillator model. Ann. Phys., 134(1):1–44, June 1981. [157] W. Nazarewicz and J. Dobaczewski. Dynamical symmetries, multiclustering, and oc- tupole susceptibility in superdeformed and hyperdeformed nuclei. Phys. Rev. Lett., 68(2):154–157, January 1992. [158] Z. Szymański and W. Nazarewicz. Rotating pseudo-oscillator scheme: pseudo-spin symmetry and identical bands. Phys. Lett. B, 433(3):229–235, August 1998. [159] M. Radomski. Nuclear rotational current and velocity fields, the cranking model, and transverse electroexcitation of the first excited state of 12 C. Phys. Rev. C, 14(5):1704– 1708, November 1976. [160] P. Gulshani and D. J. Rowe. Quantum mechanics in rotating frames. II. The lattice structure of current circulations for a rotating single-particle fluid. Can. J. Phys., 56(4):480–484, April 1978. [161] J. Kunz and U. Mosel. Flow patterns for collective rotations in heavy nuclei. Nucl. Phys. A, 323(2):271–284, July 1979. [162] J. Fleckner, J. Kunz, U. Mosel, and E. Wuest. Current distributions for rotating nuclei. Nucl. Phys. A, 339(2):227–237, April 1980. 127 [163] M. Durand, P. Schuck, and J. Kunz. Semiclassical description of currents in normal and superfluid rotating nuclei. Nucl. Phys. A, 439(2):263–288, June 1985. [164] H. Laftchiev, D. Samsœn, P. Quentin, and I. N. Mikhailov. Equivalence of pairing cor- relations and intrinsic vortical currents in rotating nuclei. Phys. Rev. C, 67(1):014307, January 2003. [165] J. Bartel, K. Bencheikh, and P. Quentin. Currents, spin densities and mean-field form factors in rotating nuclei: a semi-classical approach. Int. J. Mod. Phys. E, 13(01):225– 233, February 2004. [166] A. V. Afanasjev and H. Abusara. From cluster structures to nuclear molecules: The role of nodal structure of the single-particle wave functions. Phys. Rev. C, 97(2):024329, February 2018. [167] K. Petrík and M. Kortelainen. Thouless-Valatin rotational moment of inertia from linear response theory. Phys. Rev. C, 97(3):034321, March 2018. [168] M. Freer, R. R. Betts, and A. H. Wuosmaa. Relationship between the deformed harmonic oscillator and clustering in light nuclei. Nucl. Phys. A, 587(1):36–54, May 1995. [169] J. M. Yao, N. Itagaki, and J. Meng. Searching for a 4α linear-chain structure in excited states of 16 O with covariant density functional theory. Phys. Rev. C, 90(5):054307, November 2014. [170] J.-P. Ebran, E. Khan, T. Nikšić, and D. Vretenar. Density functional theory studies of cluster states in nuclei. Phys. Rev. C, 90(5):054329, November 2014. [171] M. Gennari and P. Navrátil. Nuclear kinetic density from ab initio theory. Phys. Rev. C, 99(2):024305, February 2019. [172] J. Dobaczewski, W. Nazarewicz, J. Skalski, and T. Werner. Nuclear Deformation: A Proton-Neutron Effect? Phys. Rev. Lett., 60(22):2254–2257, May 1988. [173] J. Dobaczewski and J. Skalski. Quadrupole collective models from the Hartree-Fock standpoint. Phys. Rev. C, 40(2):1025–1039, August 1989. [174] P. Rozmej, S. Cwiok, and A. Sobiczewski. Is octupole deformation sufficient to describe the properties of “octupolly” unstable nuclei? Phys. Lett. B, 203(3):197–199, March 1988. [175] A. Sobiczewski, Z. Patyk, S. Cwiok, and P. Rozmej. Study of the potential energy of “octupole”-deformed nuclei in a multidimensional deformation space. Nucl. Phys. A, 485(1):16–30, August 1988. 128 [176] J. L. Egido and L. M. Robledo. Microscopic study of the octupole degree of freedom in the radium and thorium isotopes with Gogny forces. Nucl. Phys. A, 494(1):85–101, March 1989. [177] J. L. Egido and L. M. Robledo. A self-consistent approach to the ground state and lowest-lying negative-parity state in the barium isotopes. Nucl. Phys. A, 518(3):475– 495, November 1990. [178] S. Ćwiok and W. Nazarewicz. Ground-state shapes and spectroscopic properties of Z ∼ 58, N ∼ 88 nuclei. Nucl. Phys. A, 496(2):367–384, May 1989. [179] S. Ćwiok and W. Nazarewicz. Reflection-asymmetric shapes in transitional odd-A th isotopes. Phys. Lett. B, 224(1):5–10, June 1989. [180] S. Ćwiok and W. Nazarewicz. Reflection-asymmetric shapes in odd-A actinide nuclei. Nucl. Phys. A, 529(1):95–114, July 1991. [181] W. Nazarewicz. Static multipole deformations in nuclei. Prog. Part. Nucl. Phys., 28:307–330, January 1992. [182] N. Schunck, J. Dobaczewski, J. McDonnell, J. Moré, W. Nazarewicz, J. Sarich, and M. V. Stoitsov. One-quasiparticle states in the nuclear energy density functional theory. Phys. Rev. C, 81(2):024316, February 2010. [183] W. Satuła, D. J. Dean, J. Gary, S. Mizutori, and W. Nazarewicz. On the origin of the Wigner energy. Phys. Lett. B, 407(2):103–109, August 1997. [184] Particle Data Group. Review of Particle Physics. Prog. Theor. Exp. Phys., 2020(8):083C01, August 2020. [185] J. T. Suhonen. Value of the Axial-Vector Coupling Strength in β and ββ Decays: A Review. Front. Phys., 0, 2017. [186] J. Engel and J. Menéndez. Status and future of nuclear matrix elements for neutrinoless double-beta decay: a review. Rep. Prog. Phys., 80(4):046301, March 2017. [187] A. B. Migdal. Theory of finite Fermi systems and application to atomic nuclei. Wiley- Interscience Publ., 1967. [188] I. S. Towner. Quenching of spin matrix elements in nuclei. Phys. Rep., 155(5):263–377, November 1987. [189] K.-F. Liu, H. Luo, Z. Ma, Q. Shen, and S. A. Moszkowski. Skyrme-Landau param- eterization of effective interactions (I). Hartree-Fock ground states. Nucl. Phys. A, 534(1):1–24, November 1991. 129 [190] T. Lesinski, K. Bennaceur, T. Duguet, and J. Meyer. Isovector splitting of nucleon effective masses, ab initio benchmarks and extended stability criteria for Skyrme energy functionals. Phys. Rev. C, 74(4):044315, October 2006. [191] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer. A Skyrme parametrization from subnuclear to neutron star densities. Nucl. Phys. A, 627(4):710– 746, December 1997. [192] M. Kortelainen, T. Lesinski, J. Moré, W. Nazarewicz, J. Sarich, N. Schunck, M. V. Stoitsov, and S. Wild. Nuclear energy density optimization. Phys. Rev. C, 82(2):024313, August 2010. [193] S. O. Bäckman, O. Sjöberg, and A. D. Jackson. The role of tensor forces in Fermi liquid theory. Nucl. Phys. A, 321(1):10–24, May 1979. [194] F. Osterfeld. Nuclear spin and isospin excitations. Rev. Mod. Phys., 64(2):491–557, April 1992. [195] T. Wakasa, M. Ichimura, and H. Sakai. Unified analysis of spin isospin responses of nuclei. Phys. Rev. C, 72(6):067303, December 2005. [196] T. Wakasa, M. Okamoto, M. Dozono, K. Hatanaka, M. Ichimura, S. Kuroita, Y. Maeda, H. Miyasako, T. Noro, T. Saito, Y. Sakemi, T. Yabe, and K. Yako. Complete sets of polarization transfer observables for the 208 Pb(⃗p, ⃗n) reaction at 296 MeV and Gamow- Teller and spin-dipole strengths for 208 Pb. Phys. Rev. C, 85(6):064606, June 2012. [197] J. Yasuda, M. Sasano, R. G. T. Zegers, H. Baba, D. Bazin, W. Chao, M. Dozono, N. Fukuda, N. Inabe, T. Isobe, G. Jhang, D. Kameda, M. Kaneko, K. Kisamori, M. Kobayashi, N. Kobayashi, T. Kobayashi, S. Koyama, Y. Kondo, A. J. Kraszna- horkay, T. Kubo, Y. Kubota, M. Kurata-Nishimura, C. S. Lee, J. W. Lee, Y. Mat- suda, E. Milman, S. Michimasa, T. Motobayashi, D. Muecher, T. Murakami, T. Naka- mura, N. Nakatsuka, S. Ota, H. Otsu, V. Panin, W. Powell, S. Reichert, S. Sak- aguchi, H. Sakai, M. Sako, H. Sato, Y. Shimizu, M. Shikata, S. Shimoura, L. Stuhl, T. Sumikama, H. Suzuki, S. Tangwancharoen, M. Takaki, H. Takeda, T. Tako, Y. Togano, H. Tokieda, J. Tsubota, T. Uesaka, T. Wakasa, K. Yako, K. Yoneda, and J. Zenihiro. Extraction of the Landau-Migdal parameter from the Gamow-Teller Giant Resonance in 132 Sn. Phys. Rev. Lett., 121(13):132501, September 2018. [198] M. Ichimura, H. Sakai, and T. Wakasa. Spin–isospin responses via (p,n) and (n,p) reactions. Prog. Part. Nucl. Phys., 56(2):446–531, April 2006. [199] H. Akimune, I. Daito, Y. Fujita, M. Fujiwara, M. B. Greenfield, M. N. Harakeh, T. Inomata, J. Jänecke, K. Katori, S. Nakayama, H. Sakai, Y. Sakemi, M. Tanaka, and M. Yosoi. Direct proton decay from the Gamow-Teller resonance in 208 Bi. Phys. Rev. C, 52(2):604–615, August 1995. 130 [200] C. Gaarde, J. Rapaport, T. N. Taddeucci, C. D. Goodman, C. C. Foster, D. E. Bainum, C. A. Goulding, M. B. Greenfield, D. J. Hören, and E. Sugarbaker. Excitation of giant spin-isospin multipole vibrations. Nucl. Phys. A, 369(2):258–280, October 1981. [201] K. Pham, J. Jänecke, D. A. Roberts, M. N. Harakeh, G. P. A. Berg, S. Chang, J. Liu, E. J. Stephenson, B. F. Davis, H. Akimune, and M. Fujiwara. Fragmentation and split- ting of Gamow-Teller resonances in Sn(3 He,t)Sb charge-exchange reactions, A=112– 124. Phys. Rev. C, 51(2):526–540, February 1995. [202] Evaluated nuclear structure data file (ENSDF), URL: https://www.nndc.bnl.gov/ ensdf/. [203] B. D. Anderson, J. N. Knudson, P. C. Tandy, J. W. Watson, R. Madey, and C. C. Foster. Observation of a T> Gamow-Teller state in 48 Ca(p, n) 48 Sc at 160 MeV. Phys. Rev. Lett., 45(9):699–702, September 1980. [204] B. D. Anderson, T. Chittrakarn, A. R. Baldwin, C. Lebo, R. Madey, P. C. Tandy, J. W. Watson, B. A. Brown, and C. C. Foster. Particle-hole strength excited in the 48 Ca(p,n)48 Sc reaction at 134 and 160 MeV: Gamow-Teller strength. Phys. Rev. C, 31(4):1161–1172, April 1985. [205] K. Yako, M. Sasano, K. Miki, H. Sakai, M. Dozono, D. Frekers, M. B. Greenfield, K. Hatanaka, E. Ihara, M. Kato, T. Kawabata, H. Kuboki, Y. Maeda, H. Matsubara, K. Muto, S. Noji, H. Okamura, T. H. Okabe, S. Sakaguchi, Y. Sakemi, Y. Sasamoto, K. Sekiguchi, Y. Shimizu, K. Suda, Y. Tameshige, A. Tamii, T. Uesaka, T. Wakasa, and H. Zheng. Gamow-Teller Strength Distributions in 48 Sc by the 48 Ca(p, n) and 48 Ti(n, p) Reactions and Two-Neutrino Double-β Decay Nuclear Matrix Elements. Phys. Rev. Lett., 103(1):012503, July 2009. [206] K. Yako, H. Sagawa, and H. Sakai. Neutron skin thickness of 90 Zr determined by charge exchange reactions. Phys. Rev. C, 74(5):051303, November 2006. [207] J. Toivanen, J. Dobaczewski, M. Kortelainen, and K. Mizuyama. Error analysis of nuclear mass fits. Phys. Rev. C, 78(3):034306, September 2008. [208] R. T. Birge. The Calculation of Errors by the Method of Least Squares. Phys. Rev., 40(2):207–227, April 1932. [209] G. A. F. Seber and C. J. Wild. Nonlinear Regression. John Wiley & Sons, Inc., February 1989. [210] J. R. Donaldson and R. B. Schnabel. Computational Experience With Confidence Re- gions and Confidence Intervals for Nonlinear Least Squares. Technometrics, 29(1):67– 82, February 1987. 131 [211] K. Pearson. LIII. On lines and planes of closest fit to systems of points in space. Lond. Edinb. Dublin Philos. Mag. J. Sci., 2(11):559–572, November 1901. [212] T. Hastie, R. Tibshirani, and J. Friedman. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer Series in Statistics. Springer, New York, 2009. [213] I. T. Jolliffe. Choosing a Subset of Principal Components or Variables. In Principal Component Analysis, pages 111–149. Springer, New York, NY, 2002. [214] S. M. Wild. Chapter 40: POUNDERS in TAO: Solving Derivative-Free Nonlinear Least-Squares Problems with POUNDERS. In Advances and Trends in Optimization with Engineering Applications, MOS-SIAM Series on Optimization, pages 529–539. Society for Industrial and Applied Mathematics, April 2017. [215] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. M. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang. PETSc Web page, URL: https://petsc.org/, 2021. [216] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang. PETSc/TAO users manual. Technical Report ANL-21/39 - Revision 3.16, Argonne National Labo- ratory, 2021. [217] S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith. Efficient management of parallelism in object oriented numerical software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, Modern software tools in scientific computing, pages 163–202. Birkhäuser Press, 1997. [218] S. M. Wild, J. Sarich, and N. Schunck. Derivative-free optimization for parameter esti- mation in computational nuclear physics. J. Phys. G: Nucl. Part. Phys., 42(3):034031, February 2015. [219] R. Bollapragada, M. Menickelly, W. Nazarewicz, J. O’Neal, P.-G. Reinhard, and S. M. Wild. Optimization and supervised machine learning methods for fitting numerical physics models without derivatives. J. Phys. G: Nucl. Part. Phys., 48(2):024001, Jan- uary 2021. 132 [220] P. Gysbers, G. Hagen, J. D. Holt, G. R. Jansen, T. D. Morris, P. Navrátil, T. Pa- penbrock, S. Quaglioni, A. Schwenk, S. R. Stroberg, and K. A. Wendt. Discrepancy between experimental and theoretical β-decay rates resolved from first principles. Nat. Phys., 15(5):428–431, May 2019. [221] E. M. Ney, J. Engel, and N. Schunck. Two-body weak currents in heavy nuclei. Phys. Rev. C, 105(3):034349, March 2022. [222] M. C. Kennedy and A. O’Hagan. Bayesian calibration of computer models. J. R. Stat. Soc. Ser. B Stat. Methodol., 63(3):425–464, 2001. [223] C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. Adaptive Computation and Machine Learning series. MIT Press, Cambridge, MA, USA, November 2005. [224] M. Chen, T. Li, B. Schuetrumpf, P.-G. Reinhard, and W. Nazarewicz. Three- dimensional Skyrme Hartree-Fock-Bogoliubov solver in coordinate-space representa- tion. Comput. Phys. Commun., 276:108344, July 2022. [225] M. Frigo and S. G. Johnson. The Design and Implementation of FFTW3. Proc. IEEE, 93(2):216–231, February 2005. [226] R. W. Hockney. The potential calculation and some applications. Methods Comput. Phys., 9:135–211, 1970. [227] J. W. Eastwood and D. R. K. Brownrigg. Remarks on the solution of poisson’s equation for isolated systems. J. Comput. Phys., 32(1):24–38, July 1979. [228] S. G. Johnson. Notes on FFT-based differentiation, URL: https://math.mit.edu/ ~stevenj/fft-deriv.pdf, May 2011. [229] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK users’ guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, third edition, 1999. [230] L. S. Blackford, J. Choi, A. Cleary, E. D’Azevedo, J. Demmel, I. Dhillon, J. Dongarra, S. Hammarling, G. Henry, A. Petitet, K. Stanley, D. Walker, and R. C. Whaley. ScaLA- PACK users’ guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, 1997. [231] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical recipes 3rd edition: The art of scientific computing. Cambridge University Press, 2007. 133 [232] J. C. Pei, M. V. Stoitsov, G. I. Fann, W. Nazarewicz, N. Schunck, and F. R. Xu. Deformed coordinate-space Hartree-Fock-Bogoliubov approach to weakly bound nuclei and large deformations. Phys. Rev. C, 78(6):064306, December 2008. [233] R. J. Furnstahl, G. Hagen, and T. Papenbrock. Corrections to nuclear energies and radii in finite oscillator spaces. Phys. Rev. C, 86(3):031301, September 2012. [234] S. Binder, A. Ekström, G. Hagen, T. Papenbrock, and K. A. Wendt. Effective field theory in the harmonic oscillator basis. Phys. Rev. C, 93(4):044332, April 2016. [235] S. N. More, A. Ekström, R. J. Furnstahl, G. Hagen, and T. Papenbrock. Universal properties of infrared oscillator basis extrapolations. Phys. Rev. C, 87(4):044326, April 2013. [236] M. Bender, K. Rutz, P.-G. Reinhard, and J. A. Maruhn. Pairing gaps from nuclear mean-field models. Eur. Phys. J. A, 8(1):59–75, July 2000. [237] T. Papenbrock and G. F. Bertsch. Pairing in low-density Fermi gases. Phys. Rev. C, 59(4):2052–2055, April 1999. 134