APPLICATION OF NUCLEAR DENSITY FUNCTIONAL THEORY TO EXOTIC NUCLEI By Mengzhi Chen 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 NUCLEAR DENSITY FUNCTIONAL THEORY TO EXOTIC NUCLEI By Mengzhi Chen Nuclear density functional theory (DFT) is the method of choice to study the nuclear properties of medium-mass and heavy nuclei. This dissertation employs the Skryme Hartree- Fock-Bogoliubov (HFB) approach to study nuclear reflection-asymmetric deformations and collective rotation. Nuclear ground states with stable reflection-asymmetric shapes, predicted by theory, have been confirmed experimentally. To explore the microscopic origin of reflection-asymmetric nuclear shapes, we applied the density expansion method to decompose the total HFB energy into different multipolarities. We demonstrated that the reflection-asymmetric deformation is driven by the isoscalar part of the interaction energy. We also confirmed the importance of high-multipolarity fields for stabilizing reflection-asymmetric deformations. The nucleon localization function (NLF) has been successfully applied to characterize nuclear shell structure and collective motion. In our work, we extended the application of NLF to study the nuclear response to fast rotation. By solving the cranked harmonic- oscillator and comparing it with cranked Hartree-Fock results, we defined the simplified localization measure and demonstrated its usefulness as an indicator of nuclear rotation. The above nuclear DFT calculations were performed using existing HFB solvers. How- ever, the current HFB solvers are deficient in the study of exotic nuclei whose properties are strongly affected by the quasiparticle continuum space. For this purpose, we developed a three-dimensional Skyrme-HFB solver HFBFFT in the coordinate-space representation using the canonical basis approach. We implemented the soft energy cutoff and pairing annealing to solve the problem of pairing collapse; a sub-iteration method to improve the convergence; and an algorithm to restore the Hermiticity of differential operators brought by Fourier-transform-based differentiation. The accuracy and performance of HFBFFT were tested by benchmarking it against other HFB codes for a set of nuclei. ACKNOWLEDGMENTS My deepest gratitude goes first to my advisor, Witold Nazarewicz, for his academic guidance and support. He is always open to new ideas and discussions. He showed me how to do research, how to write a scientific paper, and provided me opportunities to attend conferences and summer schools. I would like to thank my other committee members, Shanker Balasubramaniam, B Alex Brown, Dean Lee, and Brian O’Shea, for presenting my annual committee meetings. They always care about my progress and provide great questions and helpful advice. I am grateful to my collaborators, Tong Li, Paul-Gerhard Reinhard, Bastian Schuetrumpf, Jacek Bobaczewski, Chunli Zhang, and Markus Kortelainen, for the great projects we com- pleted together. I would also like to thank other members of our research group. They are Sylvester Agbe- mava, Yuchen Cao, Eric Flynn, Kevin Fossez, Pablo Giuliani, Samuel Giuliani, Kyle God- bey, Yannen Jaganathen, Rahul Jain, Daniel Lay, Xingze Mao, Zachary Matheson, Nicholas Michel, Futoshi Minato, Léo Neufcourt, Erik Olsen, Jimmy Rotureau, Mookyong Son, Simin Wang, and Joshua Wylie. I enjoy the time we share our ideas, achievements, and snacks during the group meeting. Special thanks to Pablo Giuliani and Kyle Godbey for reading and providing constructive comments on my manuscript. I am fortunate to have many excellent graduate students and visiting scholars at MSU who enriched my life: Long Cheng, Xi Dong, Hao Lin, Rongzheng He, He Huang, Shiyun Hu, Jian Liu, Felix Ndayisabye, Pierre Nzabahimana, Avik Sarkar, Jason Surbrook, Tommy Tsang, Kang Yu, Boyao Zhu, and Kuan Zhu. Thank you for your friendship. In my first and second years, I shared my apartment with my roommates Ankang Yang and Zhe Zhang iv whose company I will miss forever. I would like to express gratitude towards several professors: H. Metin Aktulga, Morten Hjorth-Jensen, Filomena Nunes, Scott Pratt, Jianliang Qian, Yuyin Xie, Vladimir Zelevinsky. From their classes and discussions, I learned a lot of valuable and helpful knowledge. I appreciate the help from other staff and faculty members at NSCL and MSU physics department, especially Kim Crosslan, Elizabeth Deliyski, Gillian Olsen, Artemis Spyrou, and Remco Zegers. I am also indebted to my girlfriend Mingxu. Thanks for your love and company during the hard COVID period. I am grateful to my parents for their love throughout my whole life. Only with their unconditional support can I focus on my interest and finish my graduate student career. v TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x LIST OF ALGORITHMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Overview of nuclear properties . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Nuclear shapes and deformations . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Nuclear rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4 Exotic Nuclei . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.5 Organization of this dissertation . . . . . . . . . . . . . . . . . . . . . . . . . 4 Chapter 2 Nuclear Density Functional Theory . . . . . . . . . . . . . . . . . 5 2.1 General formalism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 The Skyrme energy density functional . . . . . . . . . . . . . . . . . . . . . 6 2.3 Hartree-Fock-Bogoliubov theory . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3.1 Formalism in canonical basis . . . . . . . . . . . . . . . . . . . . . . . 12 2.4 Nuclear DFT solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Chapter 3 Microscopic origin of reflection-asymmetric nuclear shapes . . 17 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2 Multipole expansion of densities and HFB energy . . . . . . . . . . . . . . . 19 3.2.1 Multipole decomposition . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.2.2 Isospin and neutron-proton energy decomposition . . . . . . . . . . . 24 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3.1 Multipole expansion of the deformation energy . . . . . . . . . . . . . 28 3.3.2 Isospin and neutron-proton structure of the octupole deformation energy 30 3.3.3 Reflection-asymmetric deformability along isotopic chains . . . . . . . 31 3.3.4 Relation to shell structure . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4 Summary: Origin of reflection-asymmetric deformations . . . . . . . . . . . . 38 Chapter 4 Nucleon localization function in rotating nuclei . . . . . . . . . 43 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.2 Cranked Harmonic Oscillator calculations . . . . . . . . . . . . . . . . . . . . 44 4.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.3.1 General considerations . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.3.2 Simplified nucleon localization function . . . . . . . . . . . . . . . . . 46 4.3.3 Angular-momentum alignment: Cranked harmonic-oscillator analysis 48 4.4 Summary: Nucleon localization function in rotating nuclei . . . . . . . . . . 51 vi Chapter 5 Development of Three-dimensional Skyrme HFB solver . . . . 52 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.2 Numerical representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.2.1 Numerical realization on a 3D coordinate-space grid . . . . . . . . . . 53 5.2.2 Solution by accelerated gradient iteration . . . . . . . . . . . . . . . . 56 5.2.3 Sub-iterations in configuration space . . . . . . . . . . . . . . . . . . 59 5.2.4 Soft cutoff on pairing-active space . . . . . . . . . . . . . . . . . . . . 60 5.2.5 Strategies to avoid premature pairing breakdown . . . . . . . . . . . 61 5.2.6 Hermiticity restoration . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.2.7 Numerical realization in harmonic-oscillator basis . . . . . . . . . . . 65 5.3 Benchmarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.3.1 Parameter determination . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.3.2 Pairing renormalization . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.3.3 Energy shift by Hermiticity restoration . . . . . . . . . . . . . . . . . 68 5.3.4 Doubly magic nuclei: 132 Sn and 208 Pb . . . . . . . . . . . . . . . . . 69 5.3.5 Spherical superfluid nucleus: 120 Sn . . . . . . . . . . . . . . . . . . . 72 5.3.6 Axially deformed nuclei: 102,110 Zr . . . . . . . . . . . . . . . . . . . . 73 5.3.7 Superdeformed heavy nucleus: 240 Pu . . . . . . . . . . . . . . . . . . 74 5.4 Summary: Development of Three-dimensional Skyrme HFB solver . . . . . . 77 Chapter 6 Conclusions and Outlook . . . . . . . . . . . . . . . . . . . . . . . 79 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 vii LIST OF TABLES Table 5.1: Total energies Etot (in MeV) and ∆S (in MeV) for five nuclei calculated with HFBFFT without and with the Hermiticity restoration. The dig- its which do not coincide before and after the Hermiticity restoration are marked in bold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 Table 5.2: Energy contributions (in MeV) to the binding energy for 132 Sn computed with HFBTHO, HFBFFT, Sky1D, and Sky2D. The digits which do not coincide with HFBFFT are marked in bold. . . . . . . . . . . . . . . . . 69 Table 5.3: Energies for 208 Pb from HFBTHO (computed with different numbers of HO shells N ), HFBFFT, and Sky1D. All energies are in MeV. The digits which do not coincide with HFBFFT are marked in bold. . . . . . . . . . 70 Table 5.4: Results of HFB + SLy4 calculations for 120 Sn using HFBTHO, HFBFFT, Sky1D, and Sky2D. Two neutron pairing renormalization variants are con- sidered, by adjusting the neutron pairing strengths in HFBFFT, Sky1D, and Sky2D to reproduce the HFBTHO values of Ẽkin n and ∆n . All ener- gies are in MeV. The radius rrms is in fm. The digits which do not coincide with HFBFFT are marked in bold. . . . . . . . . . . . . . . . . . . . . . 72 Table 5.5: Results of HFB + SLy4 calculations for 110 Zr with HFBTHO, HFBFFT and Sky2D. Two neutron pairing renormalization variants are considered, by adjusting the neutron pairing strengths in HFBFFT and Sky2D to reproduce the HFBTHO values of Ẽkin n and ∆n . All energies are in MeV. p,n The radius rrms is in fm and quadrupole moments Q20 are in fm2 . The HFB proton pairing vanishes in this nucleus. The digits which do not coincide with HFBFFT are marked in bold. . . . . . . . . . . . . . . . . 74 Table 5.6: Results of HFB + SLy4 calculations for 102 Zr using HFBTHO, HFBFFT and Sky2D. The pairing renormalization is carried out by adjusting the proton and neutron pairing strengths in HFBFFT and Sky2D to repro- duce the HFBTHO values of ∆n and ∆p . All energies are in MeV. The p,n radius rrms is in fm and quadrupole moments Q20 are in fm2 . The digits which do not coincide with HFBFFT are marked in bold. . . . . . . . . . 75 viii Table 5.7: Results of HFB + SLy4 calculations for 240 Pu ground state and fission isomer using HFBTHO, HFBFFT and Sky2D. The pairing strengths in HFBFFT and Sky2D were adjusted to reproduce the spectral pairing gaps obtained in HFBTHO for the g.s. and f.i. separately. The s.p. space for HFBFFT is defined by means of (Ωn , Ωp ). All energies are in MeV, p,n rrms is in fm, and Q20 are in fm2 . The digits which do not coincide with HFBFFT are marked in bold. . . . . . . . . . . . . . . . . . . . . . . . . 76 ix LIST OF FIGURES Figure 3.1: The deformation energies, ∆E(β3 ) = E(β3 ) − E(β3 = 0), as functions of (0) β3 for 224 Ra (dashed lines) and 146 Ba (solid lines) calculated at β2 with the SLy4 (circles) and UNEDF2 (triangles) EDFs. . . . . . . . . . . . . . 27 Figure 3.2: Convergence of Ediff (λ) (3.26) for 224 Ra computed with SLy4 at β3 =0.05 (dashed line) and 0.15 (solid line). . . . . . . . . . . . . . . . . . . . . . 28 Figure 3.3: Multipole components, ∆E[λ] (β3 ) = E[λ] (β3 ) − E[λ] (β3 = 0), of the total deformation energy shown in Fig. 3.1, plotted for λ = 0 − 3 as functions (0) of the octupole deformation β3 at β2 . Upper (lower) panels show results for 224 Ra (146 Ba) obtained with the SLy4 (left) and UNEDF2 (right) EDFs. 29 Figure 3.4: Similar to Fig. 3.3 but for different isospin and neutron-proton components of the octupole energy ∆E[3] . . . . . . . . . . . . . . . . . . . . . . . . . 30 (0) Figure 3.5: Equilibrium quadrupole deformations β2 as functions of N for the iso- topic chains of (a) Ba, (b) Ra, (c) U, and (d) Yb computed with the SLy4 EDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Figure 3.6: Similar to Fig. 3.5 but for the deformation energy ∆E = E(β3 = 0.05) − E(β3 = 0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 Figure 3.7: Similar to Fig. 3.5 but for the deformation energies ∆E[λ] = E[λ] (β3 = 0.05) − E[λ] (β3 = 0) for λ = 0 − 3. . . . . . . . . . . . . . . . . . . . . . 34 Figure 3.8: Similar to Fig. 3.7 but for the deformation energies ∆E = E(β3 = 0.05) − E(β3 = 0) with multipole components summed up from λ = 0 to λmax . The values of λmax are listed in the legend. The regions of deformed isotopes exhibiting reflection-asymmetric instability in Fig. 3.6 are marked by shading. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 Figure 3.9: Similar to Fig. 3.8 but for the cumulative sum involving odd-λ multipoles only. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 Figure 3.10: Similar to Fig. 3.8 but for the cumulative sum involving even-λ multipoles only. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 x Figure 3.11: Single-particle (canonical) neutron (top) and proton (bottom) SLy4-HFB levels as functions of β2 (β3 = 0) for 176 Yb. Solid (dashed) lines indicate positive- (negative-) parity levels. Fermi levels εF at N = 106 and Z = 70 are marked by dash-dotted lines. The equilibrium deformation of 176 Yb is indicated by a vertical dotted line. . . . . . . . . . . . . . . . . . . . . 41 Figure 3.12: Similar to Fig. 3.11 but for 224 Ra. Fermi levels for even-even Ra isotopes with N = 130−144 are marked by circles. They have been shifted accord- ing to the position of the spherical 2g9/2 neutron and 1h9/2 proton shell. The equilibrium deformation of 224 Ra is indicated by a vertical dotted line. 42 Figure 4.1: Single-particle Routhians of the SD CHO model belonging to the super- shells Nshell = 6 and 7. The CHO quantum numbers [n1 , n2 , n3 ] are given in brackets. Positive-parity and negative-parity states are indicated by solid and dashed lines, respectively. q The rotational frequency ω is ex- pressed in units of ω0 = 3 ωz ω⊥ 2 while the Routhians E in units of ℏω . z Each level is 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 marked by the arrow. . . . . . . . . . . . . . 45 Figure 4.2: Current density j in the x-z (y = 0) plane, calculated in the CHO model with 60 particles in a SD HO well for four values of rotational frequency ω (in units of ω0 ). The magnitude |j| (in fm−4 ) is shown by color and line thickness. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 Figure 4.3: C (top), C τ (middle), and their difference (bottom) in the x-z (y = 0) plane, calculated in the CHO model with 60 particles in a SD HO well for five values of rotational frequency ω (in units of ω0 ). . . . . . . . . . . . 48 Figure 4.4: C τ (top), τ (in fm−5 , middle) and τ TF (in fm−5 , bottom) in the x-z (y = 0) plane, calculated in the CHO model with 60 particles in a SD HO well. The first column shows the reference plots at ω = 0 while the other columns show the rotational dependence relative to the ω = 0 reference as a function of ω (in units of ω0 ). . . . . . . . . . . . . . . . . . . . . . 49 Figure 4.5: Changes in the kinetic-energy density τ due to p-h excitations (at ω = 0) from the SD shell Nshell = 6 to the next supershell Nshell = 7 in Fig. 4.1. These excitations are induced in the CHO description of a 60-particle system by the cranking term. The rightmost panel shows the uniform average of individual p-h contributions. . . . . . . . . . . . . . . . . . . 50 Figure 5.1: Etot as a function of L for 208 Pb. The HFBTHO results are marked by red dots. The blue curve is fitted according to Eq. (5.26). . . . . . . . . . 71 xi LIST OF ALGORITHMS h i Algorithm 1 Compute the one-dimensional position-varing differetiation d dx B(x) dψ dx . 64 xii Chapter 1 Introduction 1.1 Overview of nuclear properties Nuclear physics aims at understanding atomic nuclei and their structure and interactions. The study of nuclear physics starts from the discovery of radioactivity by Henri Becquerel in 1896 [1]. The nucleus was first discovered by Ernest Rutherford and his assistants through the famous gold foil experiment [2]. In the beginning, the size and charge properties were studied by scattering experiments. After the discovery of the neutron by Chadwick in 1932 [3], scientists could calculate the nuclear binding energies by comparing the nuclear mass with the individual masses of the proton and neutron constituents. The first model being proposed to explain global properties of the nucleus was the liquid drop model (LDM) by Weizsäcker in 1935 [4]. The LDM successfully explains many nuclear features, including the general trend of nuclear binding energy and nuclear fission. However, the LDM is a classical approach that cannot describe quantum effects. For example, it is unable to explain the existence of enhanced binding energy for certain numbers of protons and neutrons which are known as magic numbers. Therefore, microscopic theoretical approaches considering the quantum mechanical effects have been developed. They can be roughly divided into three groups: the ab initio methods based on inter-nuclear forces [5–7], the nuclear shell model (configuration interaction method) [8–10] and the nuclear density 1 functional theory (DFT) based on the self-consistent mean field approach [11]. In the nuclear DFT approach, the nucleonic densities are the fundamental degrees of freedom instead of the nucleons. Compared with the other two methods, the advantage of nuclear DFT lies in the region of medium-mass and heavy nuclei. In my dissertation, I studied the origin of octupole deformations and nuclear rotation in heavy nuclei. 1.2 Nuclear shapes and deformations The atomic nucleus can exhibit various shapes beyond the simple spherical shape, such as prolate, oblate, and pear-shaped. Nuclear deformations result from the nuclear Jahn–Teller effect, which is a spontaneous geometric distortion mechanism for lowering the ground-state (g.s.) energy [12, 13]. The first experimental evidence of quadrupole deformation was observed from the optical spectra hyperfine structure, indicating the presence of nuclear quadrupole moments. [14, 15]. The microscopic origin of nuclear quadrupole deformations was studied in Refs. [16– 18]. Based on self-consistent Hartree-Fock (HF) calculations, he authors concluded that the quadrupole deformation can be traced back to the neutron-proton interaction. There is an ample experimental evidence confirming the existence of g.s. reflection- asymmetric nuclear shapes [19, 20]. According to theoretical calculations [21], such pear- shaped systems exists in different regions of the nuclear chart. Pear-shaped nuclei are signif- icant in searching for permanent atomic electric-dipole moments (EDM) because the atomic EDM can be enhanced by several orders of magnitude for an odd nucleus with static octupole deformation [22]. A non-zero EDM indicates a time-reversal (or equivalently charge-parity) violation and its magnitude can constrain the extensions to the standard model [20, 23]. 2 In our project [24], we extended the density expansion method [18] to study the micro- scopic origin of reflection-asymmetric nuclear shapes. The details of this work are presented in Ch. 3 1.3 Nuclear rotation Nuclear collective motion, such as rotations and vibrations, provides rich information about the nuclear structure and response to external fields. The observation of rotational bands in atomic nuclei has provided us with many insights into nuclear deformations, and the underlying shell structure [25–28]. In nuclear structure research, the nucleon localization function (NLF) is a useful tool to identify clusters in light nuclei [29–31] and nuclear reactions [32, 33], fission [34–38] and nuclear pasta phases in neutron stars [30]. We used the NLF to study the nuclear response to fast rotation in [39]. To this end, I employed the cranked harmonic-oscillator model compare the results with the cranked HF method. 1.4 Exotic Nuclei Exotic nuclei with extreme neutron-to-proton ratios are crucial for theoretical nuclear struc- ture research as their properties provide critical information on nuclear interactions, many- body techniques, and astrophysical scenarios. However, because of their weak binding, their quasiparticle excitations are often affected by the low-lying scattering space (a.k.a. particle continuum), which increases the necessary computational effort. For such nuclei, nucleonic pairing must be handled within the full Hartree-Fock-Bogoliubov (HFB) scheme instead of the simpler Bardeen-Cooper-Schrieffer (BCS) approximation [40–44]. In addition, the asso- 3 ciated self-consistent densities are usually very extended in space, which requires large basis sets or large coordinate-space boxes. Both requirements become particularly demanding for self-consistent methods, especially if one aims at symmetry-unrestricted calculations (i.e., without imposing space reflection or axial or spherical symmetries). In our paper [45], we proposed a reliable and efficient computational scheme to solve the HFB equations on a three-dimensional (3D) Cartesian coordinate-space grid. We developed the DFT solver HFBFFT in the coordinate-space representation using the canonical basis approach. Also, we benchmarked HFBFFT for several spherical and deformed nuclei against other solvers. The novel features of the solver and benchmark results are discussed in Ch. 5 1.5 Organization of this dissertation This dissertation is organized as follows. The general nuclear DFT formalism and solvers are introduced in Ch. 2. The microscopic study of reflection-asymmetric nuclear shapes is presented in Ch. 3. In Ch. 4, applications of the NLF to rotating nuclei is discussed. Ch. 5 presents the work of 3D Skyrme Hartree-Fock-Bogoliubov (HFB) solver HFBFFT. Finally, conclusions and outlook are given in Ch. 6. 4 Chapter 2 Nuclear Density Functional Theory The density functional theory (DFT) was developed as a computational method to investi- gate electronic many-body systems. Its extension to nuclear physics describes protons and neutrons in terms of nucleonic densities and currents. In nuclear DFT, the degrees of freedom of the nuclear many-body problem can be reduced significantly. Because of this, the nuclear DFT method has a robust scalability with the atomic mass number. As a result, all nuclear systems across the whole nuclear landscape can be described by this method. The advantages of the nuclear DFT are highlighted in the heavy-mass region where other methods such as the current ab initio approaches cannot be easily applied except for some special cases [46]. 2.1 General formalism The DFT originates from the framework of the two Hohenberg–Kohn theorems [47]. The first theorem states that the ground state (g.s.) properties of a many-body system can be uniquely determined by particle density ρ(r) which is only defined by three spatial coordinates. This theorem dramatically reduces the complexity of a N −body system to a manageable level. The second theorem defines the critical property that the density which minimizes the g.s. energy ρmin (r) is the actual density of the many-fermion system. There are analogies but also significant differences between electronic and nuclear sys- 5 tems. Nuclear systems are self-bound with no external potential. Moreover, the nuclear interaction is short-ranged, spin-dependent, and fairly complex formula when compared with the simple two-body Coulomb interaction. Therefore, the main ingredient of nuclear DFT is the effective interaction expressed in terms of an energy density functional (EDF) of local/non-local densities, currents, and their derivatives. Examples of widely used nuclear EDFs include the zero-range Skryme EDF [48], the finite-range Gogny EDF [49–51], and covariant EDFs [52–55]. The Skyrme functional is widely applied for its simple expression and high predictive power [56–58]. Therefore, the Skyrme EDFs are used in my work. 2.2 The Skyrme energy density functional The HFB theory describes a many-fermion system in terms of an orthonormal set of single particle (s.p.) wave functions ψα and occupation amplitudes vα , i.e., {ψα , vα , α = 1, ..., Ω} , (2.1) where Ω denotes the size of the active s.p. space. The occupation amplitude vα can take p values in the interval [0, 1]. The emptiness amplitude is uα = 1 − vα2 . Any self-consistent mean-field theory starts from expressing the energy of the system in terms of s.p. wave functions and occupation amplitudes (2.1). Among many EDFs, the Skyrme functional, originally based on the Skyrme interaction [48], is commonly used to study global nuclear properties, such as ground-state energies, deformations, and low-lying excitations [56–58]. The Skyrme EDF can be splitted into a sum of time-even and time-odd 6 terms (ignoring here isospin index for simplicity): HSkyrme (r) = E even + E odd , (2.2) E even = C ρ ρ2 + C τ ρτ + C ∆ρ ρ∆ρ + C J ρ∇ · J + C J J2 , (2.3) E odd = C s s2 + C ∆s s∆s + C T s · T + C j j 2 + C ∇j s · (∇ × j), (2.4) where in the time-even terms, ρ, τ are the particle and kinetic energy density respectively; J is the spin-orbit density and J is the spin-current tensor. In the time-odd terms, s, and T are the spin density and spin kinetic energy density, respectively, and j is the momentum density. Definitions and discussion of these densities can be found in , e.g., Ref. [59]. The C parameters in front of bi-linear densities are the coupling constants defining the Skyrme EDF. In many cases, one is interested in stationary states of even-even nuclei. In this specific instance, the EDF depends only on three energy densities ρq , τq , and J q : vα2 |ψα (r, s)|2 , XX ρq (r) = α∈q s vα2 |∇ψα (r, s)|2 , XX τq (r) = α∈q s vα2 ψα∗ (r, s)∇×σ ss′ ψα (r, s′ ), XX J q (r) = −i (2.5) α∈q ss′ where q ∈ {p, n} stands for protons or neutrons and s, s′ = ±1/2 label the two spinor components of the wave functions. Additionally, pairing EDFs require the pairing density X X ξq (r) = uα vα (−2s) ψα (r, −s)ψα (r, s), (2.6) α∈q s 7 where the first summation includes a cutoff in the pairing space. For a stationary state of an even-even nucleus, the conjugate s.p. state α can be assumed to be the time-reversed state of α, which leads to a simplified expression uα vα |ψα (r, s)|2 . XX ξq (r) = (2.7) α∈q s The Skyrme EDF is well described in several works [11,56,57]. Within the Skyrme HFB scheme, the total energy is a functional of the local densities: Etot = ESkyrme [ρ, τ, J ] + ECoul [ρp ] + Epair [ρ, ξ], (2.8) where Z  2  ℏ ρ τ J ESkyrme = d3 r 2 ∆ρ τ + C ρ + C ρτ + C ρ∆ρ + C ρ∇ · J , (2.9) 2m  1 e2 ρp (r)ρp r ′ 3e2 3 3 4/3 Z  Z ECoul = 3 d rd r 3 ′ − d r3 ρp , (2.10) 2 |r − r ′ | 4 π   1 ρ Z 3 2 X Epair = Vpair,q d r|ξq | 1 − . (2.11) 4 ρ0,pair q∈{p,n} ESkryme is a functional of densities ρ, τ , and J . Most of time-even coupling constants C in Eq. 2.9 are real numbers except for C ρ which is density-dependent: ρ ρ C ρ = C 0 + CT ρ γ . (2.12) Therefore, the Skyrme functionals of stationary even-even nuclei can be described by 11 8 parameters: ρ ρ {C0 , CT , C τ , C ∆ρ , C J }t=0,1 and γ. (2.13) The Coulomb energy ECoul is composed of direct and exchange terms as functionals of proton density ρp . The pairing energy Epair in form of 2.11 is a functional of densities ρ and ξ; this is motivated by a density-dependent δ interaction. The pairing strengths Vpair,q are adjustable parameters (see Sec. 5.3.2 for detailed discussion). The Eq. 2.11 covers two limiting cases. The first case is a pure contact interaction, also called the volume pairing, which is recovered when ρ0,pair → ∞. The second case corresponds to a value near matter equilibrium density ρ0,pair = 0.16 fm−3 , which localizes pairing around the nuclear surface. Adjustment of ρ0,pair as a free parameter delivers a form of the pairing functional which stays in between the extremes of volume and surface pairing [60, 61]. 2.3 Hartree-Fock-Bogoliubov theory The g.s. wave functions and energy of a nuclear system can be obtained by Hartree-Fock and Hartree-Fock-Bogoliubov theories based on variational principles [62]. The HFB theory generalizes the HF theory to the case of pairing correlations. The HFB g.s. wave functions can be written as a quasi-particle (q.p.) vacuum Y † |Ψ⟩ = β̂k |0⟩ (2.14) k † where |0⟩ is the s.p. vacuum, β̂k and β̂k are q.p. operators. They are related with particle 9 † operators ĉl and ĉl via the Bogoliubov linear transformation † † X β̂k = Ulk ĉl + Vlk ĉl , l (2.15) ∗ ĉ † + Vlk∗ ĉl . X β̂k = Ulk l l In the above equations, both indices l and k run over the whole one-particle configuration space. We can also rewrite the Eq. 2.15 in the matrix form        β̂k   U† V†  ĉk  ĉ †  k  =   = W   (2.16) † † † β̂k V T UT ĉk ĉk by defining the Bogoliubov unitary transformation matrix W † . The Bloch and Messiah theorem [63] states that W can be decomposed into three matrices with a special form     D 0  Ū V̄  C 0  W=   . (2.17) 0 D ∗ V̄ Ū 0 C ∗ It is worth noting that the first matrix defines a unitary transformation of the particle operators among themselves † † X αk = Dlk ĉl . (2.18) L † The new particle operators αk and αk define a new basis named the canonical basis in which the Bogoliubov transformation becomes a Bardeen-Cooper-Schrieffer (BCS) transformation. More properties of the HFB theory in the canonical basis will be discussed in the section 2.3.1. In the HFB theory, the density matrix ρ̂ and the pairing tensor κ̂ uniquely determine 10 the wave function |Ψ⟩. They are defined as † ρll′ = ⟨Ψ|ĉ ′ ĉl |Ψ⟩ = V ∗ V T , (2.19) l κll′ = ⟨Ψ|ĉl′ ĉl |Ψ⟩ = −U V † . (2.20) The pairing density ξl = (−2s)κll (2.6) is usually more convenient to use in actual calcula- tions for time-even q.p. states. A two-body Hamiltonian H of a fermion system can be written as † † † X X H= til ĉi ĉj + v̄ijlk ĉi ĉj ĉl ĉk . (2.21) il ijlk The first term is the kinetic energy, and the second represents the two-body interaction. The total energy E is an energy functional of ρ̂ and κ̂ 1 E[ρ, κ] = Tr[hρ̂] − Tr[∆κ̂∗ ], (2.22) 2 where h and ∆ are the HF and pairing mean-fields respectively 1X 1X hij = tij + v̄iljk ρkl , ∆ij = v̄ijlk κlk . (2.23) 2 2 kl kl Variation of the total energy with respect to ρ̂ and κ̂ results in the HFB equations:      h − λ ∆  U  U     = E   (2.24) −∆∗ −h∗ + λ V V where λ is the Lagrange multiplier to constrain the average particle number. The self- 11 consistency of the HFB equations is apparent because the fields h and ∆ are determined by the eigenvectors U and V . Therefore, the HFB equations are non-linear and need to be solved iteratively. 2.3.1 Formalism in canonical basis The canonical basis is defined as the basis of s.p. states ψα in which the one-body density P matrix ρ̂ is diagonal, i.e., ρ̂ = α |ψα ⟩nα ⟨ψα |, where nα , an eigenvalue of ρ̂, represents the canonical-state occupation. The numerical HFB scheme in the canonical basis was presented in [64] and improved in [65]. For the relation between the standard matrix formulation and the canonical formulation of HFB, see Refs. [42, 66]. In the canonical basis, the HFB mean- field state takes the BCS-like form: uα + vα â+ + Y |Φ⟩ = α âα |0⟩ (2.25) α>0 where |0⟩ is the vacuum state, â+ α is the creation operator of ψα , and α the conjugate partner to state α that corresponds to the same eigenvalue of ρ̂. In practice, one deals with two types of fermions: protons and neutrons. To keep the notation simple, in the following, we assume that the isospin quantum number is included in the quantum label α of the canonical state. The HFB equations are derived variationally by minimizing the HFB Routhian: vα2 − X X X  R = Etot − ϵF,q λαβ ⟨ψβ |ψα ⟩ − δαβ , (2.26) q∈{p,n} α∈q αβ with respect to ψα and vα . In Eq. (2.26) ϵF is the Fermi energy which is also the Lagrange pa- 12 rameter for the particle-number constraint, and the λ̂ is the matrix of Lagrangian multipliers that guarantee the orthonormality of canonical wave functions. Since ⟨ψβ |ψα ⟩ = ⟨ψα |ψβ ⟩∗ , it is required that the matrix λ̂ is Hermitian so that the number of its independent elements coincides with the total number of independent constraints. Variation of the Skryme and Coulomb energies with regard to the s.p. wave function yields the HF Hamiltonian ĥ:  δ ESkyrme + ECoul † = vα2 ĥψα . (2.27) δψα By the chain rule for derivatives, (2.27) can be reduced to the variation with respect to the densities, which delivers explicit expressions for ĥ [11,57] in terms of of local densities in the standard fashion of nuclear EDFs [11]. The variation of the pairing energy with respect to the s.p. wave function gives δEpair ˆ 2 ′ † = uα vα h̃ψα + vα ĥ ψα . (2.28) δψα The first term is related to the variation with respect to the pairing density, which yields the pairing potential [66]   1 ρ h̃q (r) = Vpair,q ξq 1 − , q ∈ {p, n}. (2.29) 2 ρ0,pair The second term is the pairing-rearrangement term, brought by the density dependence of the pairing functional. For simplicity, we treat the rearrangement term ĥ′ as part of the HF Hamiltonian ĥ. The pairing potential h̃q (r) is local. From Eq. 2.23, we can obtain the 13 state-dependent pairing gap ∆αα = ⟨ψα |h̃qα |ψα ⟩ , (2.30) where the subscript qα is the isospin of state α. Another aspect of the pairing is determined by the gap equation, which is obtained from the energy variation with respect to vα :  2  vα 0 = 4vα (hαα − ϵF,qα ) + 2 − uα ∆αα , (2.31) uα where hαα are the diagonal matrix elements of the HF Hamiltonian. The HF Hamiltonian together with the pairing potential constitutes the main ingredients of the HFB equations. With the orthonormality of canonical states taken into account, the constrained varia- † tion of the total energy with respect to ψα yields the mean-field equations: P Ĥα ψα = β ψβ λβα , (2.32) where ˆ Ĥα = vα2 ĥ + uα vα h̃, (2.33a) 1 λβα = ⟨ψβ |Ĥα + Ĥβ |ψα ⟩. (2.33b) 2 The mean-field equations (2.32,2.33) and the gap equations (2.31) constitute the self-consistent HFB equations in the canonical basis. In (2.33a) Ĥα is a state-dependent one-body Hamiltonian composed of the HF Hamil- tonian and the pairing potential. The full matrix λ̂ needs to be taken into account because Ĥα is state-dependent [64, 65]. In contrast, pure HF or HF+BCS calculations only require 14 diagonal matrix elements λαα , which are also known as s.p. energies. The Hermiticity of λ̂ is enforced by explicit symmetrization in Eq. (2.33b). It can be shown by multiplying both † sides of Eq. (2.32) by ψβ that the final solution should obey the symmetry conditions 1  0= λ− βα ≡ ⟨ψβ |Ĥα |ψα ⟩ − ⟨ψβ |Ĥβ |ψα ⟩ . (2.34a) 2 One can combine these into one condition: v 1 X u u 1 X − 2 0 = ∆S ≡ t 2 λβα . (2.34b) 2 Ωq q∈{p,n} α,β∈q The actual size of ∆S serves as a check for the convergence of the HFB solution. It is to be noted that λ̂− vanishes when both s.p. states ψα and ψβ are fully occupied (vα = vβ = 1) or unoccupied (vα = vβ = 0) since ⟨ψα |ĥ|ψβ ⟩ = ⟨ψβ |ĥ|ψα ⟩∗ . Thus, for a pure HF calculation, ∆S measures the overlap between occupied and unoccupied orbits, which should be zero at the self-consistent solution. 2.4 Nuclear DFT solvers Since the HFB equations are self-consistent, they need to be solved iteratively. Over the years, a number of HFB solvers have been developed; a recent summary can be found in Table 2 of Ref. [67]. These solvers can be divided into two families. The codes belonging to the first group are based on the expansion of s.p. wave functions in a finite set of basis functions such as the harmonic oscillator (HO) eigenfunctions. Solvers of this type include: HFBTHO [68,69], which solves the axial (2D) HFB equations in the axial HO or the transformed HO basis; HFODD [70–72], which solves the 3D HFB equations in the 15 Cartesian HO basis without assumption of self-consistent symmetries; and HFBPTG [73], which solves the axial HFB equations in the Pöschl-Teller-Ginocchio basis. The basis-expansion method is efficient and has been successfully employed in large-scale calculations [58]. However, when it is applied to weakly-bound nuclei, the performance of this method deteriorates as huge configuration spaces are required to describe the asymptotic behavior of HFB solutions. Here, the approach of choice is the HFB framework formulated in the coordinate-space representation [40, 42, 43]. The coordinate-space solvers constitute the second family of HFB codes. Examples of such solvers are: HFBRAD [74] that solves spherically symmetric HFB problem using finite differences; HFB-AX [75] is a 2D solver based on B-splines; SkyAx [76] is a highly optimized 2D HF + BCS code using the fast Fourier transform (FFT) method to compute derivatives; Sky3D [77, 78] is a 3D extension of SkyAx; the predecessor of SkyAx and Sky3D is a 1D spherical HF+BCS code using five-point finite differences which was published first in [79] and has meanwhile been developed into a full spherical HFB code Sky1D [80]; the HFB extension of SkyAx is Sky2D [80]; EV8 solves the Skyrme HF+BCS equations using the imaginary time method on a 3D mesh that is limited to one octant by imposing time- reversal and spatial symmetries [81, 82]; MOCCa [83, 84] is a Skyrme-HFB extension of EV8; MADNESS-HFB [85] is a 3D HFB solver based on multi-resolution analysis and multi-wavelet expansion; LISE is a 3D HFB solver [86] employing the discrete variable representation (or Lagrange-mesh method) and fast Fourier transforms; and there are also 3D HFB solvers based on the contour integral of the Green’s function using the shifted Krylov subspace method [87, 88]. 16 Chapter 3 Microscopic origin of reflection-asymmetric nuclear shapes 3.1 Introduction While the most majority of atomic nuclei have g.s. shapes that are either spherical or ellip- soidal (prolate or oblate), some isotopes exhibit pear-like shape deformations that intrinsi- cally break reflection symmetry. Characteristic features of nuclear spectra, nuclear moments, and electromagnetic matrix elements provide experimental support for such shapes [19, 20]. Low-energy negative-parity excitations in pear-shaped even-even nuclei are commonly at- tributed to octupole collective modes. As a result of this reason, pear-shaped nuclei are commonly referred to as "octupole-deformed". There are two regions of g.s. reflection-asymmetric shapes that have been experimentally established over the years: the neutron-deficient actinides around 224 Ra and the neutron-rich lanthanides around 146 Ba. Nuclear theory systematically predicts these nuclei to be pear- shaped (for a recent summary of theoretical results, see Ref. [21]). Other regions of pear- shaped nuclei predicted by theory, i.e., lanthanide nuclei around 200 Gd as well as actinide and superheavy nuclei with 184 < N < 206 are too neutron rich to be accessible by experiment [21,58,89–91]. Deformation energy associated with reflection-symmetry breaking shapes are, 17 on general, substantially lower than those associated with stable ellipsoidal shapes [92, 93]. As a result, beyond-mean-field approaches are required for a quantitative description of octupole-deformed nuclei; see, for example, Refs. [94–97]. According to the s.p. picture, the emergence of pear-shaped deformations can be at- tributed to the mixing of opposite-parity s.p. shells [98, 99]. In the macroscopic-microscopic (MM) approach, the macroscopic energy favors spherical shapes. Therefore, stable refection- asymmetric shape deformations obtained in the MM method [93, 100] can be traced back to the shape polarization originating from proton and neutron s.p. levels interacting via parity- breaking fields. The results are usually explained in terms of the deformation-driving proton or neutron shell effects because shell corrections are computed independently for protons and neutrons. The proton-neutron interactions are indirectly considered in the macroscopic energy with the assumption of identical proton and neutron shape deformation parameters, which follow those of the macroscopic term. In general, in the description based on the mean-field approach, nuclear shape defor- mations result from a coupling between collective surface vibrations of the nucleus and valence nucleons. The nuclear Jahn-Teller effect [13, 16] can be used to explain such a particle-vibration coupling mechanism [101]. The tendency towards deformation is partic- ularly strong if the Fermi level lies just between close-lying s.p. states. In this instance, the system could become unstable with respect to the mode that couples these states. Simple estimates of the particle-vibration coupling (Jahn-Teller vibronic coupling) for the quadrupole mode (multipolarity λ = 2) [102, 103] demonstrate that its contribution to the mass quadrupole moment at low energies doubles the quadrupole moment of valence nucle- ons. The HF analysis [17, 18] confirmes this estimate. In particular, it has been shown that the attractive isoscalar quadrupole-quadrupole term, which can be well approximated by the 18 neutron-proton quadrupole interaction, provides the largest contribution to the quadrupole deformation energy. When it comes to reflection-asymmetric deformations, the octupole mode (multipolarity λ = 3) has the strongest particle-vibration coupling. This coupling generates a vibronic Jahn-Teller interaction between close-lying opposite-parity s.p. orbits that may result in a static reflection-asymmetric shape. Such pairs of states can be discovered just above closed shells for g.s. configurations of atomic nuclei and involve a unique-parity intruder shell (ℓ, j) and a normal-parity shell (ℓ − 3, j − 3) around particle numbers Noct = 34, 56, 88, and 134 [19]. Indeed self-consistent calculations systematically predict pear shapes for nuclei having proton and neutron numbers close to Noct . To understand the origin of reflection-asymmetric g.s. deformations, in this study, we extend the quadrupole-energy analysis of Refs. [17, 18] to odd-multipolarity shapes. To this end we decompose the total HFB energy into isoscalar, isovector, neutron-neutron (nn), proton-proton (pp), and neutron-proton (np) contributions of different multipolarities. 3.2 Multipole expansion of densities and HFB energy In self-consistent mean-field approaches [11, 56, 62] with EDFs based on two-body functional generators, the total energy of a nucleus can expressed as: E = Tr(T ρ) + 12 Tr(Γρ) + 12 Tr(Γ̃ξ). (3.1) Here T is the kinetic energy operator, Γ and Γ̃ are mean fields in particle-hole (p-h) and particle-particle (p-p) channels, respectively, and ρ and ξ are one-body p-h and p-p density 19 matrices, respectively. (Instead of using the standard pairing tensor [62], here we use the “tilde” representation of the p-p density matrix [40].) The mean fields Γ and Γ̃ are defined as δE T +Γ= ′ , (3.2) δρ δE Γ̃ = ′ , (3.3) δξ where δ ′ denotes the variation of the total energy that neglects the dependence of the func- tional generators on density, that is, the mean fields (3.2) and (3.3) do not contain rear- rangement terms [11]. 3.2.1 Multipole decomposition As observed from Ref. [18], the density matrices and mean fields can be split into different multipole components as ρ = ρ[0] + ρ[1] + ρ[2] + ρ[3] + . . . (3.4a) ξ = ξ[0] + ξ[1] + ξ[2] + ξ[3] + . . . , (3.4b) Γ = Γ[0] + Γ[1] + Γ[2] + Γ[3] + . . . (3.4c) Γ̃ = Γ̃[0] + Γ̃[1] + Γ̃[2] + Γ̃[3] + . . . , (3.4d) where ρ[λ] , ξ[λ] , Γ[λ] , and Γ̃[λ] are rank-λ rotational components of ρ, ξ, Γ, and Γ̃, respectively. Traces appearing in Eq. (3.1) are invariant with respect to unitary transformations, and, in particular, with regard to spatial rotations. Therefore, the trace act like a multipolarity filter projecting the total energy onto a rotational scalar. In this way, when the multipole 20 expansions (3.4) are inserted in the expression for the total energy (3.1), only diagonal terms survive: E = E[0] + E[1] + E[2] + E[3] + . . . , (3.5) where E[λ] = 21 Tr(Γ[λ] ρ[λ] ) + 12 Tr(Γ̃[λ] ξ[λ] ). (3.6) In the above equation, we include the kinetic energy in the monopole energy E[0] since T is a scalar operator which implies Ekin = Tr(T ρ) ≡ Tr(T ρ[0] ). Therefore we can define E[0] = Ekin + 21 Tr(Γ[0] ρ[0] ) + 21 Tr(Γ̃[0] ξ[0] ). (3.7) When parity symmetry is conserved, only even-λ multipolarities appear in expansions (3.4) and (3.5). In Refs. [17,18], this allowed for the analysis in terms of the monopole (λ = 0), quadrupole (λ = 2), and higher even-λ components. In our work, we study broken-parity self-consistent states with an emphasis on the reflection-asymmetric (odd-λ) components of the expansion. It is worth noting that the integral of the isoscalar dipole density ρ[1] , i.e. the entire isoscalar dipole moment, vanishes by construction because our multipole expansion is specified with respect to the nucleus’s center of mass. Nevertheless, the dipole density ρ[1] and dipole energy E[1] can still be nonzero. In the spherical s.p. basis, the expansions (3.4) can be achieved by the angular-momentum coupling of basis wave functions. An explicit basis transformation is required since the HFB equations are usually solved in a deformed basis. Moreover, the direct angular-momentum coupling does not benefit from the fact that Skyrme EDFs only depend on (quasi)local densities, which is the property that significantly simplifies the HFB problem. Inspired by 21 the latter observation, in this work, we determine the multipole expansions of (quasi)local densities and (quasi)local mean fields directly in the coordinate space. With axial symmetry assumed, particle density ρ(r) can be decomposed as [104] X ρ(r) = ρ[λ] (r)YJ,M =0 (Ω), (3.8) J where Z ρ[λ] (r) = dΩρ(r)YJ,M ∗ =0 (Ω). (3.9) An identical decomposition can be carried out for all isoscalar (t = 0) and isovector (t = 1) (quasi)local p-h densities [105] ϱt ≡ {ρt , τt , ∆ρt , Jt , ∇ · Jt }, plus local neutron (q = n) and proton (q = p) pairing densities ξq . The p-h densities depend on neutron and proton densities in the usual way: ϱ0 = ϱn + ϱp , ϱ1 = ϱn − ϱp . (3.10) Our strategy is to use the energy-density expression for the time-even total energy (3.1), Z  2  ℏ even E= d3 r τ (r) + E (r) + Ẽ(r) , (3.11) 2m 0 where the standard Skyrme energy densities read [105, 106]: E even (r) = Eteven (r), X (3.12a) t=0,1 X Ẽ(r) = Ẽq (r), (3.12b) q=p,n 22 where E even is defined in the Eq. 2.3 and ρ(r) γ 2     Ẽq (r) = 1V 1 − V1 ξq (r). (3.13) 4 q ρ0 For simplicity, the Coulomb energy is not included in Eq. (3.11). It is convenient to rewrite the energy densities (2.3) in terms of local p-h and p-p potentials as Eteven (r) = Vt (r)ρt (r) + X Vtij (r)Jtij (r), (3.14a) ij Ẽq (r) = Ṽq (r)ξq (r), (3.14b) where ρ ∆ρ Vt (r) = Ct ρt (r) + Ct ∆ρt (r) + Ctτ τt (r) + Ct∇J ∇ · Jt (r), (3.15a) Vtij (r) = CtJ Jtij (r), (3.15b) ρ(r) γ     1 Ṽq (r) = 4 Vq 1 − V1 ξq (r), (3.15c) ρ0 with indices i, j denoting the components of the spin-current tensor density Jtij (r) in three dimensions. In analogy to Eqs. (3.8) and (3.9), we determine the multipole expansions of the local potentials (3.15). In this way, the total energy (3.5) can be decomposed into multipole 23 components:    Z X   E[λ] = d3 r  X Vt[λ] (r)ρt[λ] (r) + Vtij[λ] (r)Jtij[λ] (r)   t=0,1 ij # (3.16) X + Ṽq[λ] (r)ξq[λ] (r) . q=p,n Finally, the same strategy can also be applied to the Coulomb energy, which contributes to the multipole terms of Eq. (3.5) through the multipole expansions of both direct and exchange potentials: Z h i Coul E[λ] = d3 r 12 V[λ] dir (r) + 3 V exc (r) ρ 4 [λ] p[λ] (r), (3.17) where ρp (r′ ) Z V dir (r) = e2 d3 r ′ , (3.18) |r − r ′ | h i1 V (r) = −e π ρp (r) 3 . exc 2 3 (3.19) 3.2.2 Isospin and neutron-proton energy decomposition The total energy can be expressed in the isospin and the neutron-proton schemes. In the isospin scheme, the total energy can be written as E = E t=0 + E t=1 + E Coul + E pair , (3.20) 24 where Z Et = Ekin δt0 + d3 rEteven (r), (3.21a) X Z E pair = d3 r Ẽq (r). (3.21b) q=p,n Note that the kinetic energy Ekin is included in the isoscalar energy E t=0 . The Coulomb energy E Coul is separated out because the Coulomb interaction breaks the isospin symmetry. The pairing functional is not isospin invariant either, as the neutron and proton pairing strengths are different. By decomposing the isoscalar and isovector p-h densities ϱt into the neutron and proton components (3.10), the total energy can be expressed in the neutron-proton scheme [18]: E = Ekin + E nn + E pp + E np . (3.22) ′ In Eq. (3.22), the individual E qq components (q, q ′ = n or p): Z ′ h i E qq = d3 r Eqq even (r) + δ Ẽ (r) , ′ qq ′ q (3.23) are defined through the energy densities E even ′ and Ẽq , which are bilinear in the densities ϱq qq or ξq . In this scheme, the Coulomb energy E Coul is included in the proton-proton energy E pp . As previously discussed, all of the energy terms entering the isospin and neutron-proton decompositions can be expanded into multipoles. 25 3.3 Results The systems we studied are even-even barium (Ba), radium (Ra) and uranium (U) isotopes. They are predicted to have stable pear-like shapes at certain neutron numbers [21]. For comparison, we also calculate ytterbium (Yb) isotopes which have stable quadrupole but no reflection-asymmetric deformations. We performed axial HFB calculations using the code HFBTHO (v3.00) [69] for two Skyrme EDFs given by SLy4 [107] and UNEDF2 [108] parametrizations. We used the mixed-pairing strengths of Vn = −325.25 MeV and Vp = −340.06 MeV (SLy4) and Vn = −231.30 MeV and Vp = −255.04 MeV (UNEDF2). For UNEDF2, we did not apply the Lipkin-Nogami treatment of pairing; instead, we took the neutron pairing strength Vn to reproduce the average experimental neutron pairing gap for 120 Sn, ∆n = 1.245 MeV. The proton pairing strength Vp was adjusted proportionally based on the default values of Vn and Vp . In the first step, we performed parity-conserving calculations by constraining the oc- tupole deformation to zero and determined the corresponding equilibrium quadrupole defor- (0) (0) mation β2 . At the fixed value of β2 , we varied β3 from 0.0 to 0.25. In the HFBTHO code, multipole constraints are actually applied to quadrupole (Q20 ) and octupole (Q30 ) moments related to β2 and β3 through r ! 16π 3 β2 = Q20 / AR02 , 5 4π r ! (3.24) 16π 3 β3 = Q30 / AR03 , 7 4π 26 where A is the mass number, R0 = 1.2 fm×A1/3 , and D E 2 2 Q20 = 2z − x − y ,2 D  E (3.25) Q30 = z 2z 2 − 3x2 − 3y 2 . 146 Ba UNEDF2 224 Ra SLy4 ΔE (MeV) β3 Figure 3.1: The deformation energies, ∆E(β3 ) = E(β3 ) − E(β3 = 0), as functions of β3 for 224 Ra (dashed lines) and 146 Ba (solid lines) calculated at β (0) with the SLy4 (circles) and 2 UNEDF2 (triangles) EDFs. Figure 3.1 shows reflection-asymmetric deformation energies ∆E(β3 ) = E(β3 )−E(β3 = 0) determined for 224 Ra and 146 Ba obtained in this way. We see that UNEDF2 yields a higher octupole deformability than SLy4 for both nuclei. This is consistent with the results in Ref. [21]. 27 2 10 1 10 β3=0.15 Ediff (MeV) 0 10 β3=0.05 −1 10 −2 10 224Ra, SLy4 −3 10 0 1 2 3 4 5 6 7 8 9 λ Figure 3.2: Convergence of Ediff (λ) (3.26) for 224 Ra computed with SLy4 at β3 =0.05 (dashed line) and 0.15 (solid line). 3.3.1 Multipole expansion of the deformation energy The convergence of the multipole expansion (3.5) provides a check on the accuracy of our results. In Fig. 3.2, we show the energy difference, λ X Ediff (λ) = E[λ′ ] − E (3.26) λ′ =0 for 224 Ra at two values of the octupole deformation, β3 = 0.05 and 0.15. We see that at β3 = 0.15, the multipole components decrease exponentially with λ, with the monopole component off by about 150 MeV and the sum up to λ = 9 exhausted up to about 20 keV. At a small octupole deformation of β3 = 0.05, high-order contributions decrease. As expected, the octupole component brings now less energy as compared to the quadrupole one. The 28 results displayed in Fig. 3.2 convince us that cutting the multipole expansion of energy at λ = 9 provides sufficient accuracy. SLy4 UNEDF2 224 Ra ΔE[λ] (MeV) (a) (b) [0] 146 Ba [2] [1] [3] (c) (d) β3 Figure 3.3: Multipole components, ∆E[λ] (β3 ) = E[λ] (β3 ) − E[λ] (β3 = 0), of the total de- formation energy shown in Fig. 3.1, plotted for λ = 0 − 3 as functions of the octupole (0) deformation β3 at β2 . Upper (lower) panels show results for 224 Ra (146 Ba) obtained with the SLy4 (left) and UNEDF2 (right) EDFs. Figure 3.3 shows how the reflection-asymmetric deformation energy builds up. It presents the four leading multipole components ∆E[λ] (β3 ) = E[λ] (β3 ) − E[λ] (β3 = 0), for λ = 0 − 3, of the deformation energies shown in Fig. 3.1. We can see that the pattern of contributions of different multipolarities is fairly generic: it weakly depends on the choice of the nucleus or EDF. Figure 3.3 clearly demonstrates that the main driver of reflection-asymmetric shapes is a strong attractive octupole energy ∆E[3] . The attractive dipole energy ∆E[1] is much weaker. The monopole and quadrupole energies are repulsive along the trajectory of β3 (0) (with a fixed quadrupole deformation β2 ) and essentially cancel the octupole contribution. 29 Indeed, one can note that while individual multipole components can be of the order of tens of MeV, the total reflection-asymmetric deformation energy shown in Fig. 3.1 is an order of magnitude smaller. Therefore, the final reflection-asymmetric correlation results from a large cancellation between individual multipole components, and even a relatively small variation of any given component can significantly shift the net result. In addition, as discussed in Sec. 3.3.3 below, higher-order multipole components (λ > 3) can be important for the total energy balance. 3.3.2 Isospin and neutron-proton structure of the octupole defor- mation energy SLy4 UNEDF2 224 Ra ΔE[3] (MeV) (a) (b) 146 Ba (c) (d) β3 Figure 3.4: Similar to Fig. 3.3 but for different isospin and neutron-proton components of the octupole energy ∆E[3] . 30 To analyze the origin of the octupole energy ∆E[3] , in Fig. 3.4 we show its isospin and neutron-proton components as defined in Eqs. (3.21a) and (3.23). Again, a generic pattern t=0 . The emerges. In all cases, the octupole energy is almost equal to its isoscalar part ∆E[3] t=1 is indeed very small, even if the studied nuclei have a significant isovector energy ∆E[3] pair neutron excess. The contribution from the pairing energy ∆E[3] is also practically negligi- ble. In the neutron-proton scheme, the np component always clearly dominates the nn and t=0 ≈ ∆E np pp terms. The latter two are very small for UNEDF2 and hence ∆E[3] ≈ ∆E[3] [3] for this EDF. For SLy4, the nn and pp terms provide larger contributions to the octupole deformation energy, accompanied by a reduction of the np term. Regardless of these minor differences between the EDFs, we can safely conclude that it is the isoscalar octupole com- ponent (or the np octupole energy component) that plays the dominant role in building up the nuclear octupole deformation. 3.3.3 Reflection-asymmetric deformability along isotopic chains At this point, we are ready to study structural changes that dictate the appearance of nu- clear reflection-asymmetric deformations. The results shown in Figs. 3.3 and 3.4 tell us that a mutual cancellation of near-parabolic shapes of different components of the deforma- tion energy results in a clearly non-parabolic dependence of the total deformation energy, as seen in Fig. 3.1. Therefore, to trace back the positions and energies of the equilibrium reflection-asymmetric deformations to the properties of specific interaction components is not straightforward. To this end, we analyze the properties of reflection-asymmetric de- formabilities of nuclei, that is, we concentrate on the curvature of reflection-asymmetric deformation energies at β3 = 0. To investigate the variation of the reflection-asymmetric deformability with neutron number, we performed SLy4-HFB calculations for the isotopic 31 (a) (b) Ba Ra β2(β3=0) (c) U (d) Yb Neutron number (0) Figure 3.5: Equilibrium quadrupole deformations β2 as functions of N for the isotopic chains of (a) Ba, (b) Ra, (c) U, and (d) Yb computed with the SLy4 EDF. 32 chains of even-even 138−152 Ba, 214−232 Ra, and 216−234 U isotopes, which are in the region of reflection-asymmetric instability, as well as 166−180 Yb, which are expected to be reflection- (0) symmetric [21]. In Fig. 3.5 we show the baseline quadrupole deformations β2 . For the Ba, Ra, and U isotopic chains, spherical-to-deformed shape transitions are predicted slightly above the neutron magic numbers. The considered open-shell Yb isotopes are all predicted to be well deformed. Ba (a) Ra (b) ΔE (MeV) U (c) (d) Yb Neutron number Figure 3.6: Similar to Fig. 3.5 but for the deformation energy ∆E = E(β3 = 0.05) − E(β3 = 0). As a quantitative measure of the octupole deformability, we analyze the deformation energy ∆E = E(β3 = 0.05) − E(β3 = 0) calculated at a small octupole deformation of (0) β3 = 0.05, with the quadrupole deformation fixed at β2 . We have checked that for different energy components, curvatures ∆E/β32 are stable within about 1% up to β3 = 0.05, so values of ∆E taken at β3 = 0.05 constitute valid measures of the octupole stiffness. In Fig. 3.6 33 we show the values of ∆E calculated for the four studied isotopic chains. We see that the negative values of ∆E delineate regions of neutron numbers where reflection-asymmetric deformations set in in Ba, Ra, and U isotopes [21]. [0] Ba (a) Ra (b) [2] [1] ΔE[λ] (MeV) [3] U (c) Yb (d) Neutron number Figure 3.7: Similar to Fig. 3.5 but for the deformation energies ∆E[λ] = E[λ] (β3 = 0.05) − E[λ] (β3 = 0) for λ = 0 − 3. We now study ∆E[λ] , the multipole components of the total deformation energy, for the four isotopic chains considered to see whether they could provide insights into the neutron- number dependence of octupole deformations. Figure 3.7 shows that the answer is far from obvious. Indeed, we observe strong cancellations of contributions coming from different mul- tipole components of the reflection-asymmetric deformation energy. For example, both the repulsive monopole and attractive octupole components are an order of magnitude larger than the total deformation energies shown in Fig. 3.6. Therefore, we can expect that in order to understand the behavior of the deformation energies, higher-order multipole com- 34 ponents ∆E[λ] should be considered. Indeed, it has been early recognized that higher-order deformations can strongly influence the octupole collectivity of reflection-asymmetric nu- clei [109–116]. (a) (b) Ba Ra ΔE[0-λmax] (MeV) (c) (d) 0 5 U Yb 1 6 2 7 3 8 4 9 Neutron number Figure 3.8: Similar to Fig. 3.7 but for the deformation energies ∆E = E(β3 = 0.05)−E(β3 = 0) with multipole components summed up from λ = 0 to λmax . The values of λmax are listed in the legend. The regions of deformed isotopes exhibiting reflection-asymmetric instability in Fig. 3.6 are marked by shading. To better see accumulation effects with increasing multipolarity and subtle fluctuations at different orders, in Fig. 3.8 we plot multipole components of the octupole deformabil- ity summed up to λmax . Noting dramatically different scales of Figs. 3.6 and 3.8, we see that summations up to about λ = 5 or 7 are needed for the results to converge. Al- though the octupole component contributes by far the most to the creation of the reflection- asymmetric deformation energy, its effect is counterbalanced by a very large monopole com- ponent and, therefore, even higher multipole components are instrumental in determining 35 the total reflection-asymmetric deformability. This aspect is underlined in the results shown in Figs. 3.9 and 3.10, where we separately show analogous sums of only odd-λ (odd parity) and even-λ (even parity) components, respectively. It is clear that the octupole polarizability is a result of a subtle balance between positive (repulsive) effect of the even-parity multipoles and negative (attractive) effect of the odd-parity multipoles. (a) (b) Ba Ra ΔE[1,3,..λmax] (MeV) (c) (d) 1 3 5 7 Yb 9 U Neutron number Figure 3.9: Similar to Fig. 3.8 but for the cumulative sum involving odd-λ multipoles only. 3.3.4 Relation to shell structure To gain some insights into the shell effects behind the appearance of stable reflection- asymmetric nuclear shapes, Figs. 3.11 and 3.12 show, respectively, the s.p. level diagrams for 176 Yb and 224 Ra as functions of β2 . While such diagrams cannot predict symmetry breaking effects per se, they can often provide qualitative understanding. The well-deformed nucleus 176 Yb is characteristic of a stiff octupole vibrator. Indeed, 36 (a) 0 2 Ra ΔE[0,2,..λmax] (MeV) 4 6 8 Ba (b) U (c) (d) Yb Neutron number Figure 3.10: Similar to Fig. 3.8 but for the cumulative sum involving even-λ multipoles only. its nucleon numbers (Z = 70, N = 106) lie far from the “octupole-driving” numbers Noct . Due to the large deformed Z = 70 gap around β2 = 0.32, there are no s.p. states of opposite parity and the same projection Ω of the total s.p. angular momentum on the symmetry axis that could produce p-h excitations with appreciable λ = 3 strength across the Fermi level. As for the neutron s.p. levels, the low-Ω positive-parity states originating from the 1i13/2 shell lie below the Fermi level, which appreciably reduces the 1i13/2 ↔ 2f7/2 strength. Because of the large quadrupole deformations of Yb isotopes considered, the s.p. orbital angular momentum ℓ of normal-parity orbitals is fairly fragmented within the shell [117]. As seen in Figs. 3.10d and 3.9d, all multipole components of ∆E for 176 Yb vary very smoothly with neutron number. The Nilsson diagram shown in Fig. 3.12 is characteristic of transitional neutron-deficient 37 actinides in which the octupole instability is expected. The unique-parity shells, 1i13/2 pro- ton shell and 1j15/2 neutron shell, are of particle character, which results in an appearance of close-lying opposite-parity pairs of Nilsson levels with the same low Ω-values at intermediate quadrupole deformations. These levels can interact via the octupole field, with the dominant π1i13/2 ↔ π2f7/2 and ν1j15/2 ↔ ν2g9/2 couplings. As seen in Figs. 3.9 and 3.10, in the regions of octupole instability, the monopole and quadrupole deformation energies become locally reduced while the octupole and dotria- contapole (λ = 5) contributions to ∆E grow. According to our results, the effect of the dotriacontapole term is essential for lowering ∆E around Noct . This not surprising as the main contribution to the dotriacontapole coupling comes from the ∆ℓ = ∆j = 3 excita- tions [113, 115], i.e., the octupole and dotriacontapole correlations are driven by the same shell-model orbits. Interestingly, it is the attractive λ = 5 contribution to ∆E rather than the octupole term that exhibits the local enhancement in the regions of octupole instability. The shallow octupole minima predicted around 146 Ba result from an interplay between the odd-λ deformation energies, which gradually increase with N (see Fig. 3.9a) and the even-λ deformation energies, which gradually decrease with N (see Fig. 3.10b). Again, the dotriacontapole moment is absolutely essential for forming the octupole instability. 3.4 Summary: Origin of reflection-asymmetric deforma- tions In this chapter, we used the Skyrme-HFB approach to study the multipole expansion of interaction energies in both isospin and neutron-proton schemes in order to analyze their role in the appearance of reflection-asymmetric g.s. deformations. The main conclusions and 38 results of our study can be summarized as follows: 1. Based on the self-consistent HFB theory, reflection-asymmetric ground-state shapes of atomic nuclei are driven by the odd-multipolarity isoscalar (or, in neutron-proton scheme, np) part of the nuclear interaction energy. 2. The most favorable conditions for reflection-asymmetric shapes are in the regions of transitional nuclei with neutron and proton numbers just above magic numbers. For such systems, the unique-parity shell has a particle character, which creates favorable conditions for the enhanced ∆ℓ = ∆j = 3 octupole and dotriacontapole couplings. 3. The presence of high-multipolarity interaction components, especially λ = 5 are crucial for the emergence of stable reflection-asymmetric shapes. Microscopically, dotriacon- tapole couplings primarily come from the same ∆ℓ = ∆j = 3 p-h excitations that are responsible for octupole instability. According to our calculations, the attractive λ = 5 contribution to the octupole stiffness is locally enhanced in the regions of reflection- asymmetric g.s. shapes. In summary, stable pear-like g.s. shapes of atomic nuclei result from a dramatic can- cellation between even- and odd-multipolarity components of the nuclear binding energy. Small variations in these components, associated, e.g., with the s.p. shell structure, can thus be instrumental for tilting the final energy balance towards or away from the octupole instability. One has to bear in mind, however, that the shell effect responsible for the sponta- neous breaking of intrinsic parity is weak, as it is associated with the appearance of isolated ∆ℓ = ∆j = 3 pairs of levels (parity doublets) in the reflection-symmetric s.p. spectrum. In this respect, the breaking of the intrinsic spherical symmetry in atomic nuclei (presence 39 of ellipsoidal deformations) is very common as every spherical s.p. shell (except for those with j = 1/2) carries an intrinsic quadrupole moment that can contribute to the vibronic coupling. 40 (a) neutrons -3 126 -4 3p1/2 -5 2f 5/2 1i13/2 -6 3p 3/2 εF 104 -7 100 Single-particle energy (MeV) -8 1h9/2 2f7/2 -9 82 (b) protons -5 82 -6 -7 3s 1/2 70 1h11/2 εF -8 2d 3/2 -9 64 -10 2d5/2 -11 1g7/2 56 0.0 0.1 0.2 0.3 0.4 0.5 β2 Figure 3.11: Single-particle (canonical) neutron (top) and proton (bottom) SLy4-HFB levels as functions of β2 (β3 = 0) for 176 Yb. Solid (dashed) lines indicate positive- (negative-) parity levels. Fermi levels εF at N = 106 and Z = 70 are marked by dash-dotted lines. The equilibrium deformation of 176 Yb is indicated by a vertical dotted line. 41 2g7/2 3d3/2 -2 1j15/2 -3 -4 1i11/2 144 144 2g 140 142 -5 9/2 138 -6 130 132 134 136 Single-particle energy (MeV) -7 126 -8 (a) neutrons 2f5/2 -2 1i13/2 -3 2f7/2 -4 1h9/2 -5 88 130 132 134 138140 92 142 -6 136 144 -7 82 -8 (b) protons 0.0 0.1 0.2 0.3 β2 Figure 3.12: Similar to Fig. 3.11 but for 224 Ra. Fermi levels for even-even Ra isotopes with N = 130 − 144 are marked by circles. They have been shifted according to the position of the spherical 2g9/2 neutron and 1h9/2 proton shell. The equilibrium deformation of 224 Ra is indicated by a vertical dotted line. 42 Chapter 4 Nucleon localization function in rotating nuclei 4.1 Introduction As discussed in Sec. 1.4, nuclear rotation strongly depends on nuclear shell structure. Al- though rotation is essentially a time-dependent problem, the introduction of a rotating in- trinsic frame through the cranking approximation transforms the time-dependent problem into a time-independent one [118]. The cranking term added to the nuclear Hamiltonian can be interpreted as a constraint on the angular momentum, with the rotational frequency playing the role of the Lagrange multiplier. The spatial electron localization function (ELF) was originally introduced in the context of electronic HF studies to characterize shell structure in atoms and chemical bonds in molecules [119–124]. In nuclear structure research, the nucleon localization function (NLF) has been applied to identify clusters, fission fragments , and nuclear pasta [29–32, 34–38]. Compared to nucleonic densities that are fairly constant in the nuclear interior, the NLF’s characteristic oscillating pattern due to shell effects quantifies nuclear configurations much better. As a result, it is considered to be a good indicator of the interplay between collective 43 nuclear modes and s.p. motion. The NLF is defined as  !2 −1 τ (r) = 1 + Dqsµ (r) Cqsµ TF (r)  , (4.1) τqsµ 2 2 1 ∇ρqsµ j qsµ Dqsµ = τqsµ − − , (4.2) 4 ρqsµ ρqsµ TF is the Thomas- where τqsµ , ρqsµ and j qsµ are density dependent local densities, and τqsµ Fermi kinetic-energy density. More details and extensions of the NLF can be found in Ref. [39]. In the work [39], we used the NLF to study the nuclear response to rotation. We considered the case of superdeformed (SD) 152 Dy, a quintessential nuclear rotor that has been investigated in a number of self-consistent works [125–127]. My contribution to this project is to solve the cranked harmonic-oscillator (CHO) model and compare the patterns with the results of cranked HF calculations. 4.2 Cranked Harmonic Oscillator calculations In the previous study of the NLF, the HO model was used to provide an illustrative guidance [30]. In this work, we study the NLF patterns of the SD CHO model with frequencies ω⊥ = ωx = ωy = 2ωz . Since the HO potential is spin-independent, every s.p. HO level is doubly-degenerate. We assume that the rotation takes place around the y-axis. The s.p. Routhians and wave functions of the CHO can be obtained analytically [26,128,129]. We wish to emphasize that our CHO results were obtained without imposing the consistency relation between mean-field ellipsoidal deformation and the average density distribution [25, 129]. 44 [3, 0, 1] N shell = 7 carry the largest s.p. angular momentum. In Fig. 4.1 those are the [n1 ,n2 ,n3 ] = [0,0,7] (N = 7) and [0,0,6] (N = 6) Routhians. 4.3 Results and discussion 4.3.1 General considerations In a rotating system, the current density j characterizes the collective rotational behavior [127, 133–140]. Figure 4.2 shows how the current density builds up in the CHO model. As rotational frequency increases, a pattern of the vector field j resembling a rigid-body rotation gradually develops. At ω = 0.2ω0 , the lowest N = 7 Routhian [0,0,7] becomes occupied and the [3,0,0] level becomes empty, see Fig. 4.1. As the orbital [0,0,7] is strongly prolate-driving and it carries large s.p. angular momentum, and the Routhian [3,0,0] has large negative quadrupole moment (oblate), the associated configuration change (band crossing) results in a large increase in the angular momentum alignment and intrinsic deformation, see Fig. 4.2. This effect is also present in CHO calculations which consider the potential-density consistency relation [132]. 4.3.2 Simplified nucleon localization function An important consequence of the rigid-body flow is that the current density only contributes significantly to the NLF at the surface. This observation should be valid in most cases even if an irrotational flow exists (see examples in Refs. [136, 141]). The same argument is also 2 valid for the contribution to the NLF from the density gradient term ∇ρqs , which has a 46 ω = 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 4.2: Current density j in the x-z (y = 0) plane, calculated in the CHO model with 60 particles in a SD HO well for four values of rotational frequency ω (in units of ω0 ). The magnitude |j| (in fm−4 ) is shown by color and line thickness. surface character. Consequently, we define a simplified localization measure as:  !2 −1 τ (r) = 1 + τqsµ (r) Cqsµ TF (r)  , (4.3) τqsµ which does not include contributions from the current density and density gradient. Fig- ure 4.3 shows C, C τ , and their difference obtained in the CHO model; we indeed see that C τ exhibits the same pattern as C inside the nuclear volume. These two localization functions differ only in the surface region. At lower frequencies, this difference is even less pronounced. Therefore, the simplified localization function C τ is a useful tool to characterize intrinsic configurations in most cases, except perhaps for dynamic processes and high-energy modes where the current density and density gradient can become appreciable inside the nucleus. 47 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 4.3: C (top), C τ (middle), and their difference (bottom) in the x-z (y = 0) plane, calculated in the CHO model with 60 particles in a SD HO well for five values of rotational frequency ω (in units of ω0 ). 4.3.3 Angular-momentum alignment: Cranked harmonic-oscillator analysis In this section, we use the CHO model to illustrate some general features of the NLF and local densities. First, to show the usefulness of C τ when it comes to the visualization of nucleonic shell structure and angular-momentum alignment, we come back to Fig. 4.3. A characteristic regular pattern seen at ω = 0 gradually gets blurred with ω. At ω = 0.2ω0 , where the band crossing occurs, C τ rapidly changes. Namely, the number of maxima along 48 the z axis increases as the [0,0,7] orbit becomes occupies, and the number of maxima along the x axis decreases as the [3,0,0] state gets emptied. To clearly see the evolution of C τ with ω, we consider the indicator ∆C τ (r; ω) ≡ C τ (r; ω) − C τ (r; ω = 0). (4.4) This quantity is shown in Fig. 4.4 together with the corresponding variations ∆τ and ∆τ TF relative to the nonrotating case. 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 4.4: C τ (top), τ (in fm−5 , middle) and τ TF (in fm−5 , bottom) in the x-z (y = 0) plane, calculated in the CHO model with 60 particles in a SD HO well. The first column shows the reference plots at ω = 0 while the other columns show the rotational dependence relative to the ω = 0 reference as a function of ω (in units of ω0 ). One can notice that there is a clear correspondence between the peaks of ∆C τ and valleys (peaks) of ∆τ (∆τ TF ), which is consistent with Eq. (4.3). This observation suggests that 49 ∆τ and ∆τ TF are in antiphase, which results in a constructive interference when considering their ratio. a pattern in the last panel of Fig. 4.5, which is indicative of a change in τ due to rotation. Interestingly, this pattern is quite similar to that of Fig. 4.4 at ω = 0.15ω0 . We can thus conclude that, for a system that is strongly elongated along z axis, rotation-aligned s.p. states with large n3 leave a strong imprint on ∆τ and ∆C τ . 4.4 Summary: Nucleon localization function in rotating nuclei In this chapter, we solved the CHO model and used the NLF to interpret the results. The NLF involves various local densities, among which the current density j, density gradient ∇ρ, and spin-current tensor density J are appreciable only in the surface region. By neglecting surface effects, we defined a simplified localization measure C τ , which involves only the kinetic-energy density τ and the Thomas-Fermi kinetic energy density τ TF . We argue that C τ is amplified by the out-of-phase spatial oscillation of τ and τ TF attributed to the specific nodal structure of high-N s.p. states. The CHO results also demonstrate that C τ is an excellent indicator of the nodal structure for the rotating nuclei. These conclusions are consistent with the CHF analysis. 51 Chapter 5 Development of Three-dimensional Skyrme HFB solver 5.1 Introduction As we discussed in Sec. 1.3, the description of weakly bound exotic nuclei requires more com- putational effort and the pairing must be treated within the full HFB scheme. Furthermore, the corresponding self-consistent densities are usually very extended in space, necessitating the use of huge s.p. base or boxes. Two families of HFB solvers are briefly summarized in Sec. 2.4. The major difference be- tween basis-based and mesh-based methods is the treatment of one-quasiparticle continuum space [41, 43, 142, 143]. In the case of coordinate-space methods, the discretized continuum strongly depends on the geometry of the spatial box, the grid size, and the method em- ployed for the representation of derivatives. For large 3D boxes and dense grids, the size of the discretized continuum space quickly becomes intractable as the maximum allowed q.p. energy increases. As proposed in Refs. [64,65], a promising approach to solve the coordinate- space HFB problem is the canonical-basis HFB method. The one-body density matrix is diagonal in the canonical basis, and its eigenstates are spatially localized if the nucleus is particle-bound. Because of this localization, the s.p. continuum level density is significantly 52 reduced. Inspired by these features, we developed and benchmarked a 3D Skyrme-HFB solver HFBFFT in the coordinate-space representation. The theoretical framework of this project is described in Sec. 2.3.1. The new code is based on the published code Sky3D [77,78]. Sky3D has been well optimized for performance and parallelized with OpenMP and MPI [144]. In HFBFFT, we maintain a high-level parallelization, making it scalable on modern supercomputers. In order to overcome the pairing collapse problem mentioned in Ref. [65], we implemented the soft energy cutoff of pairing space and developed the annealing of pairing strengths to avoid pairing deadlock at an early stage. Furthermore, we introduced the sub-iteration method in the configuration space to stabilize and speed up the convergence. We also resolved the problem of Hermiticity violation in Sky3D brought by the incompatibility between the product rule and the Fourier- transform-based algorithm for derivatives. To test the validy and accuracy of HFBFFT, we studied several nuclear systems and benchmarked our results against HFBTHO and the coordinate-space HFB codes Sky1D and Sky2D, which solve the HFB problem in 1D (spherical) and 2D (axial) geometries, respectively. 5.2 Numerical representation 5.2.1 Numerical realization on a 3D coordinate-space grid In this section, we introduce the essentials of the numerical representation. For simplicity, our discretization strategy is explained here for one dimension as an example; the generalization to 3D is straightforward. All wave functions, densities and fields are defined on a three-dimensional equidistant 53 Cartesian grid. The grid points in the x direction are   Nx + 1 xν = − + ν δx, ν = 1, . . . , Nx , (5.1) 2 where Nx is the (even) number of grid points and δx is the grid spacing. Similar gridding applies to the y and z directions. The action of local operators on a coordinate-space grid is a simple multiplication of the local operator field and the wave function. The action of momentum operators, such as in the kinetic energy, requires first and second derivatives defined in Fourier space. The Fourier technique has been proved to be superior in precision and advantageous for large grids [145]. It is noteworthy that the direct Coulomb potential is also solved in Fourier space. The Coulomb solver has to fulfill the condition that the result in the box is the correct solution to Poisson’s equation with the boundary condition of zero potential at infinity. The algorithm to solve Poisson’s equation for an isolated charged distribution has been implemented in Sky3D. It follows the ideas of [146, 147] by doubling the 3D grid, folding the proton density with the 1/r Green’s function in momentum space and then restricting the final solution inside the original box. A subtle point is the setting of the Coulomb field at k = 0 in reciprocal space. We optimized it empirically to typical nuclear situations, for details see [77], which leaves a possible uncertainty of a few keV just at the edge of the last digit in the comparisons presented in this paper. The discrete grid points kn in Fourier space are related to the same number of grid 54 points xν in the coordinate space as:  Nx (n − 1)δk, n = 1, . . . ,   kn = 2 , (5.2a)  (n − Nx − 1)δk, Nx n= + 1, . . . , Nx  2 2π δk = . (5.2b) Nx δx Note that the coordinate-space grid (5.1) in combination with the conjugate momentum- space grid (5.2), imposes no spatial symmetry at all. But the particular examples considered for benchmarking in this study obey reflection symmetry in all three directions. A wave function ψ(xν ) in coordinate space is related to a wave function ψ(k e n ) in Fourier space by the discrete Fourier transform and its inverse XNx ψ(k e n) = exp (−ikn xν )ψ(xν ), (5.3a) ν=1 Nx 1 X ψ(xν ) = exp (ikn xν )ψ(k e n ). (5.3b) Nx n=1 Both can be efficiently computed via the FFT algorithm provided by the FFTW3 library [148]. This complex Fourier representation implies that the function ψ is periodic, i.e., ψ(x + Nx · δx) = ψ(x). The appropriate integration scheme that complies with the above summations is the trapezoidal rule Z Nx δx Nx 2 X dx f (x) ≈ f (xν )δx, (5.4) − N2x δx ν=1 where all terms are added up with equal weights. In Fourier space, the m-th derivative becomes a multiplication by (ikn )m . One proceeds 55 then in the following way: First, a forward transform (5.3a) is performed; then ψ(k e n ) is multiplied by (ikn )m ; and finally (ikn )m ψ(k e n ) is transformed back to the coordinate space by Eq. (5.3b). One should note that there is an arbitrariness about the choice of momentum kNx /2+1 : it can be taken as ± N2x δk. This arbitrariness does not alter the transforms (5.3a, 5.3b) but gives different results of the m-th derivative when m is odd (no impact when m is even). A natural choice is to equally split ψ(k e Nx /2+1 ) between the positive and negative momenta, making them cancel each other in the final result of an odd-order derivative. It is equivalent to setting ψ(k e Nx /2+1 ) = 0. This choice ensures that the derivative of a real-valued function is still real-valued; it also means that the second derivative is not equivalent to two consecutive first derivatives in this framework. The remaining problem is the Hermiticity breaking caused by the product rule; we will discuss it in Sec. 5.2.6. 5.2.2 Solution by accelerated gradient iteration The solution of the coupled HFB equations is obtained by interlaced iterations of the gap equation and the mean-field equation. The gap equation (2.31) can be solved in a closed form and it yields:   v  vα    u1 1 hαα − ϵF,qα =t ∓ q . (5.5) u 2 2 (h − ϵ 2 2  u  F,qα ) + ∆αα  α  αα The Fermi energy ϵF needs to be adjusted to fulfill the particle-number condition vα2 = Nq , X ϵF,q ←→ (5.6) α∈q 56 where Nq is the required particle number. Note that only the diagonal elements of the pairing potential and the HF Hamiltonian in the canonical basis enter (5.5); hence, no information about the non-diagonal elements is needed to determine the occupation amplitudes. The solution of the mean-field equation (2.32) is obtained by the damped gradient iteration [77, 145, 149, 150] interlaced with updating the matrix λ̂ for the orthonormality constraint. The steps are: 1. For given {ψα , vα , uα , α = 1, ..., Ω} compute the local densities, the HF Hamiltonian ĥ and the pairing potential h̃.ˆ 2. Compute the action of ĥ and h̃ ˆ on all ψ and store the result in work arrays Ψ and α α Ψe α , i.e., ĥψα −→ Ψα , (5.7a) ˆ h̃ψα −→ Ψα , (5.7b) e for α = 1, . . . , Ω. 3. Use Ψα and Ψ e α to compute and store the s.p. energies and pairing gaps hαα = ⟨ψα |Ψα ⟩, (5.8a) ∆αα = ⟨ψα |Ψe ⟩ . α (5.8b) 4. Evaluate and store the action of the generalized mean-field Hamiltonian (overwriting Ψα ) Hα ψα = vα2 Ψα + uα vα Ψe −→ Ψ . α α (5.9) 57 5. Apply the matrix of Lagrange multipliers on all ψα ; compute and store (again over- writing Ψα ) X Ψα − ψβ λβα −→ Ψα . (5.10) β 6. Apply the damping operation D̂ and orthonormalization Ô n o (new) ψα = Ô ψα − D̂Ψα , (5.11a) x0 D̂ = , (5.11b) vα2 (T̂ + E0 ) + 21 uα vα h̃0 where x0 , E0 , and h̃0 are adjustable numerical parameters. The empirical values h i x0 = 0.45, E0 = 100 MeV and h̃0 = max h̃n (r), h̃p (r) are used in our calculations. It is worth noting that the lower bound of uα vα and vα2 in D̂ is set to be 10−1 for numerical stability. 7. With the new hαα and ∆αα from step 3, compute new occupations vα and uα using Eq. (5.5). 8. Reevaluate the action of the generalized mean-field Hamiltonian on all ψα and compute the matrix of Lagrange multipliers ⟨ψβ |Ĥα |ψα ⟩ + ⟨ψα |Ĥβ |ψβ ⟩∗ λβα = . (5.12) 2 The above iteration usually starts from a number of HF+BCS steps, which are done in the same way as in Sky3D. The HF+BCS calculation is initialized by a 3D HO wave 58 function that can be triaxially deformed. To achieve better convergence, in step 1 the new densities are mixed linearly with the old ones: (n) κ(n) = (1 − γ)κ(n−1) + γκψ , κ = ρ, τ or ξ, (5.13) where n is the iteration number, subscript ψ denotes the density directly computed from the wave functions, and γ is the adjustable mixture ratio with a default value of 0.2. 5.2.3 Sub-iterations in configuration space The damped gradient scheme outlined in Sec. 5.2.2 converges, but requires more iterations in the HFB scheme as compared to the HF + BCS used in Sky3D. It also involves operations on the full 3D grid which can make computations cumbersome. The pairing part in the iterative steps works predominantly within the given space of canonical states. Thus one can reduce the total numerical expense by the sub-iteration method: switching between the full 3D step and a fast iterative solver in configuration space. To this end, we map the mean-field equations into configuration space with the expansion XΩ ψα = φn cnα , (5.14) n=1 where {φn } is a set of s.p. states acting as the expansion basis. For simplicity we choose (0) an expansion basis such that cnα = δnα at the beginning. Inserting (5.14) into the HFB mean-field equations (2.32) yields * + Ĥα − Ĥβ λ− c∗nβ X βα = φn φm cmα = 0. (5.15) mn 2 59 Eq. (5.15) is essentially the same as the symmetry condition (2.34a). It is solved by a simple damped gradient iteration:    (new)  δ X X  cnα = Ô cnα −  Hα,nm cmα − c λ nβ βα    hnn −h11 +E0 m β   (5.16)  δ  cnβ λ− X = Ô cnα − βα ,  hnn −h11 +E0  β where Hα,nm = ⟨φn |Ĥα |φm ⟩ and 1X ∗  λβα = cnβ Hα,nm + Hβ,nm cmα . (5.17) 2 mn The (interlaced) solution of the gap equations remains as before, but we do not update the local densities, the HF Hamiltonian ĥ and the pairing potential h̃ ˆ in configuration space. The convergence of the iteration is checked, again, by the symmetry condition (2.34b). The most efficient combination of the full 3D step with the iterations in configuration space is a matter of experience, see Sec. 5.3. 5.2.4 Soft cutoff on pairing-active space It is well known that the HFB equations with local interactions diverge when solved in infinite quasiparticle/canonical space [40]. To limit the pairing-active space, all local densities (2.5) are augmented by the cutoff factor wα , for instance the particle and pairing densities: wα vα2 |ψα (r, s)|2 , X X ρ(r) = (5.18a) α s |ψα (r, s)|2 . X X ξ(r) = wα uα vα (5.18b) α s 60 The same augment also applies to the kinetic-energy and the spin-orbit densities. A fixed number of states (realized by setting wα = 1 or 0) is dangerous for two reasons. First, it hinders the portability of the pairing functional between codes and nuclei, because the s.p. space depends on the basis representation. Second, level crossings near the hard cutoff can induce jumps of the pairing energy. These problems can be solved by pairing renormalization [42, 151] which, however, could be impractical in a full 3D treatment that involves huge canonical spaces. Therefore, a commonly used remedy is to use a soft pairing cutoff [152,153] 1 wα =  . (5.18c) hαα − ϵF − ∆ϵcut 1 + exp ∆ϵcut /10 The cutoff places a fixed band ∆ϵcut above the actual Fermi energy ϵF . We are going to use here ∆ϵcut = 15 MeV. It is important to note that the soft cutoff modifies the state-dependent Hamiltonian Ĥα : ˆ   Ĥα = wα vα2 ĥ + uα vα h̃ , (5.19) which defines all the ingredients entering the canonical HFB equations. 5.2.5 Strategies to avoid premature pairing breakdown The pairing comes along with a second-order superfluid-to-normal phase transition [62]. Below the critical pairing strength, the HFB pairing gap remains exactly zero. Above this critical strength, pairing becomes active and the gap starts to grow quickly. However, the onset of pairing is often delayed in a numerical calculation. The problem is that zero pairing remains a valid solution to the HFB (BCS) equations, but an unstable one. It can then take a very long time before the algorithm overcomes the instability and drives towards a 61 stable solution with nonzero pairing gap. As a consequence, an iteration scheme can easily be deadlocked due to a pairing breakdown. This is a well-known problem. Most algorithms incorporate recovery strategies, such as occasional kickoffs by giving the pairing gap an artificial value, small enough not to spoil the physics but large enough to revive the pairing mechanism. There is a more insidious problem with the state-dependent pairing gap ∆αα : it can happen that one canonical state logs out from the pairing scenario and gets stuck in its own pairing breakdown ∆αα → 0. To understand that, we inspect Eq. (2.32) and recall that ˆ . Far above the Fermi energy, we encounter states with u ≫ v such   Ĥα = vα vα ĥ + uα h̃ α α that Ĥα ≈ h̃ ˆ becomes a purely local operator. The solution to the mean field equation is ψ ∝ δ(r − r min ) where r min is the point h̃ ˆ has a minimum. In practice, this will be the representative of a δ-function on the grid, slightly mellowed by orthonormalization to other states. As a consequence, the state acquires a very high kinetic energy and a very high canonical s.p. energy, which drives the solution of the gap equation (5.5) even more toward vα → 0. This as such is a valid physical mechanism as long as the iterations curb down the occupations slowly from above. It becomes a problem if some vα gets stuck at zero at the very early stage of the iterative process. Once this has happened, the state α is locked out of the pairing space. In order to avoid this from happening, we adopt a strategy similar to simulated annealing [154] and start the iteration scheme with an enhanced effective pairing strength which gradually reduces to the physical strength as   (eff) max(Nenh − iter, 0) Vpair = Vpair ηenh +1 , (5.20) Nenh where iter is the iteration number. In practice, we use an enhancement factor ηenh = 2 62 and Nenh = 400. With this choice, the lock-in problem in the most critical early phase of iterations is avoided. 5.2.6 Hermiticity restoration According to Refs. [77, 78], the explicit expression of applying the Skyrme HF Hamiltonian ĥ on a wave function ψ can be written as: ĥψ = U (r)ψ − ∇ · [B(r)∇] ψ (5.21) i + [W · (σ × ∇) ψ + σ · ∇ × (W ψ)] . 2 This expression can be directly derived from the Skyrme EDF via Eq. (2.27), without invok- ing the product rule. In [78] it was noted that the product rule is not perfectly fulfilled when derivatives are evaluated via the discrete Fourier transform. Therefore, in Sky3D version 1.1 the commonly-adopted form of the spin-orbit term iW · (σ × ∇)ψ (5.22) was replaced by the one given in Eq. (5.21); with ∇ × W = 0, these two forms are connected by the product rule. However, the second term of ĥψ, which involves a position-varying differential operator, is still calculated through the product rule in Sky3D: X ∂B ∂ψ ∂ 2ψ ∇ · [B(r)∇] ψ = +B 2 . (5.23) ∂i ∂i ∂i i=x,y,z Unfortunately, evaluating Eq. (5.23) with the FFT-based differentiation breaks the Her- miticity of the operator [155]. This point is confirmed by computed results shown in Sec. 63 5.3.3. Instead of using Eq. (5.23), the simplest way to restore Hermiticity in the evaluation of ∇ · (B∇ψ) is to compute two consecutive first-order derivatives. But, as discussed in Sec. 5.2.1, this creates a problem with the second derivative that involves the Fourier component ψ(k e Nx /2+1 ). According to Ref. [155], one should keep the term ψ(kNx /2+1 ) in the two e first derivatives, and average the results of kNx /2+1 = ± N2x δk to maintain the symmetry in Fourier space. One can show that this “average” algorithm is equivalent to Algorithm 1 (Algorithm 3 in [155]), which is simpler to compute and thus implemented in HFBFFT. In Algorithm 1, one first computes an FFT-based first derivative, with ψ(k e Nx /2+1 ) saved and then zeroed before the inverse FFT is performed on ikn ψ(k e n ) (steps 1 through 3). Then one multiplies in coordinate space with the field B(x) involved, i.e., ϕ = Bψ ′ (step 4). Finally, one computes the derivative of the ϕ with ϕ(k e Nx /2+1 ) modified so that we can keep Hermiticity without losing the information of ψ(k e Nx /2+1 ) (steps 5 through 7). The position-varying differential operator also appears in many other physics equations, like the heat equation with varying diffusivity and Poisson’s equation with changing permittivity; hence, Algorithm 1 has a broad application range. h i d dψ Algorithm 1 Compute the one-dimensional position-varing differetiation dx B(x) dx . 1: Compute Fourier transform ψen = FFT[ψν ] with ψν = ψ(xν ). 2: Save ψ e′ e′ Nx /2+1 → Ψ, build ψ n = ikn ψn with ψ Nx /2+1 = 0. e e e 3: Compute inverse transform ψ ′ = FFT−1 [ψ ν e′ ]. n 4: Build ϕν = Bν ψν′ with Bν = B(xν ). 5: Compute Fourier transform ϕ en = FFT[ϕν ]. PNx 2 Bν Nx  ′ ′ 6: Build ϕ n = ikn ϕn and set ϕ Nx /2+1 = − ν=1 2 δk Ψ. e e e e h i Nx d B(x) dψ −1 e′ 7: Compute inverse transform dx dx ν = FFT [ϕ n ]. 64 5.2.7 Numerical realization in harmonic-oscillator basis The HFB solutions obtained with the HFBFFT code have been compared with the well- established code HFBTHO. This code has been extensively documented in several publica- tions [68,69]. The solver HFBTHO uses an expansion of the s.p. wave functions in the basis of axially symmetric HO (or transformed HO) states. The basis is given by the number of oscillator shells that defines the s.p. space size, as well as the oscillator length and deforma- tion that determine the HO wave functions. Local fields in HFBTHO are handled at the Gaussian integration points and the Gaussian integration rule is used to compute integrals. A major difference between the two codes lies in the way the HFB equations are solved. HFBFFT uses a representation in terms of the canonical basis, see Sec. 5.2, while HFBTHO works in a quasiparticle space. The results are fully equivalent if the same number of s.p. states is used. Differences appear in connection with the cutoff in the pairing space. HF- BFFT defines the cutoff in terms of the canonical s.p. energies, whereas HFBTHO does that in terms of the quasiparticle energies. This, taken together with the fact that the pairing strength has to depend on the size of the pairing space, means that the values of Vpair,q are not fully portable. This will affect in the benchmarking tests presented in Sec. 5.3. 5.3 Benchmarks In this section, we benchmark HFBFFT against HFBTHO, Sky1D, and Sky2D. These solvers have certain symmetry restrictions. Sky1D enforces spherical symmetry and can be used for magic nuclei. HFBTHO and Sky2D allow for axially symmetric shapes and cover all test cases here. Those codes can run with or without imposing reflection symmetry. First, we determine appropriate parameters to use in the calculations, including the box 65 size and grid spacing. Before making comparisons with other solvers, we quantify the effect brought by the Hermiticity restoration. In the next step, we compare some characteristic nuclei ranging from spherical doubly magic 132 Sn and 208 Pb, to spherical superfluid 120 Sn, to deformed superfluid 102,110 Zr, to superdeformed fission isomer in 240 Pu. In all these calculations, we use the Skyrme functional SLy4 [107] in the particle-hole channel and the mixed density-dependent δ interaction (ρ0,pair = 0.32 fm−3 in Eq. (2.11)) in the particle- particle channel. 5.3.1 Parameter determination To ensure the correct asymptotic behavior near the box boundary, we carry out calculations for 110 Zr to determine the appropriate box and grid sizes. The nucleus 110 Zr is chosen because it has a significant neutron excess and thus weakly bound canonical states. The calculated proton and neutron densities were inspected for different box lengths and different grids. Based on this analysis, we adopted a cubic box with a side length of 37.6 fm and 48 grid points in each dimension (spacing between two neighboring points is 0.8 fm). With the above settings, the proton and neutron densities are below 10−7 nucleons/fm3 at the boundary, which is small enough for our tests. For spherical nuclei such as 120 Sn, a smaller box is usually sufficient. We consider 176 neutron and 126 proton canonical states (Ωn = 176, Ωp = 126), and 15 MeV energy cutoff for the pairing window. This number of active states is determined by the tests for spherical 120 Sn and deformed 110 Zr nuclei. When we increase the number of active states to 200 neutron and 150 proton states, the total energy remains stable within 10 keV. In order to speed up the convergence, we perform 100 sub-iteration steps in the configuration space between two gradient iterations in the coordinate space, initialize with 66 30 HF+BCS steps, and employ the pairing enhancement factors defined in Sec. 5.2.5. For HFBTHO calculations, we take 25 HO shells for both protons and neutrons unless explicitly stated otherwise. An axially deformed HO basis with β2 = 0.2 is used in deformed ground-state calculations (102,110 Zr and 240 Pu) and β2 = 0.6 is used to calculate the 240 Pu fission isomer. For the spherical nuclei, we also compare HFBFFT results with the results of the 1D spherical HFB code Sky1D, which uses a radial coordinate-space mesh and the five-point finite difference formula for derivatives. The mesh spacing and the number of points we employ in Sky1D are 0.15 fm and 141, respectively. For the deformed nuclei, we compare HFBFFT results with the 2D axial HFB code Sky2D, which uses 31 points in both r- and z-directions with a mesh spacing of 0.7 fm. Since the nuclei considered in this study are all reflection-symmetric, the grid extends from z = 0 fm to z = 21 fm. 5.3.2 Pairing renormalization As mentioned in Sec. 5.2.7, pairing strengths are not portable between HFBFFT and HF- BTHO because of different descriptions of the pairing space and different structures of the one-quasiparticle continuum in these two solvers. Therefore, we need to renormalize the pairing strengths to compare results for open-shell nuclei for which pairing is essential. In- tuitively, there are several choices for pairing renormalization. For instance, one can tune pairing strengths to reproduce the pairing energies in different solvers. However, as discussed in [151, 156], the pairing energy density is divergent with respect to the cutoff energy. A better measure is the quantity q q q Ẽkin = Ekin + Epair (q = n or p), (5.24) 67 which is less sensitive to the pairing cutoff energy. As it will be shown in Sec. 5.3.4, the kinetic energy strongly depends on the basis size in HFBTHO. Therefore, in situations when the error related to the choice of the basis, or spatial grid, dominates, Ẽkin will be a poor renormalization measure. Another pairing measure is the spectral pairing gap [40, 42, 157] P 2 α∈q wα vα ∆αα ∆q ≡ P 2 (q = n or p). (5.25) α∈q wα vα This quantity has been used in numerous papers to adjust pairing strengths to observed odd- even mass differences and we shall use it in this study to renormalize the pairing channel for different solvers. 5.3.3 Energy shift by Hermiticity restoration As mentioned in Sec. 5.2.6, the product rule in the FFT-based differentiation violates the Hermiticity of the position-varying differential operator. To restore the Hermiticity, we implement Algorithm 1 in HFBFFT. The results are shown in Table 5.1 for several nuclei. The Hermiticity violation is demonstrated by a non-vanishing ∆S ∼ 10−7 MeV in the calculations of spherical nuclei 132 Sn and 208 Pb where the static pairing vanishes and hence the HFB calculation is reduced to HF. As for other open-shell nuclei with non-vanishing pairing, their ∆S values are similar before and after the Hermiticity restoration. These values of ∆S are characteristic of the accuracy typically achieved in HFBFFT and they are larger than the error due to the Hermiticity breaking. In terms of the total energy, the effect is of the order of a few keV, i.e., insignificant for many practical applications. Even so, Hermiticity breaking effects can affect some calculations if not remedied. For example, the small error brought by the Hermiticity breaking can accumulate step by step in a time- 68 dependent calculation. Hermiticity broken Hermiticity restored HFBFFT Etot ∆S Etot ∆S 132 Sn −1103.5429 3.44E-07 −1103.5423 1.60E-15 208 Pb −1635.6817 3.16E-07 −1635.6807 1.20E-15 120 Sn −1018.3310 3.44E-07 −1018.3305 4.01E-07 110 Zr −893.8578 4.59E-07 −893.8574 5.54E-07 102 Zr −859.4696 4.94E-07 −859.4692 3.93E-07 Table 5.1: Total energies Etot (in MeV) and ∆S (in MeV) for five nuclei calculated with HFBFFT without and with the Hermiticity restoration. The digits which do not coincide before and after the Hermiticity restoration are marked in bold. 132 208 5.3.4 Doubly magic nuclei: Sn and Pb We first calculate two doubly magic unpaired nuclei 132 Sn and 208 Pb. For these nuclei, the results of HFBFFT and Sky3D are identical. In Table 5.2, we list the ground-state energies as well as contributions from various functional terms, obtained from four solvers HFBFFT, HFBTHO, Sky1D and Sky2D for 132 Sn. Table 5.3 shows similar results for 208 Pb. When 132 Sn HFBTHO HFBFFT Sky1D Sky2D Etot −1103.49 −1103.54 −1103.57 −1103.56 n Ekin 1637.71 1637.97 1638.01 1638.02 p Ekin 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ρ∆ρ 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 Table 5.2: Energy contributions (in MeV) to the binding energy for 132 Sn computed with HFBTHO, HFBFFT, Sky1D, and Sky2D. The digits which do not coincide with HF- BFFT are marked in bold. we compare HFBFFT with Sky1D and Sky2D for 132 Sn , we find the energy differences do not usually exceed 40 keV. Such small differences can be traced back to different box 69 boundary conditions assumed in these codes. In HFBFFT, calculations are performed in a 3D rectangular box while the box is represented by a spherical shell in Sky1D and a cylindrical shape in Sky2D. For a well-bound nucleus and large spatial boxes, the results should be practically independent of the geometry of the box. As seen in Table 5.2 this indeed holds for 132 Sn. As we will see below, larger box-related errors are expected in superfluid and/or weakly bound nuclei. For nuclear matter and time-dependent calculations, the finite- size box errors can be appreciable; they can be greatly reduced by imposing twist-averaged boundary conditions [158]. 208 Pb N =15 N =20 N =25 N =30 HFBFFT Sky1D Etot −1634.25 −1635.16 −1635.46 −1635.62 −1635.68 −1635.70 n Ekin 2525.13 2527.80 2528.42 2528.83 2529.13 2529.16 p Ekin 1334.56 1336.34 1336.71 1336.91 1337.06 1337.07 Eρρ −7835.80 −7844.07 −7845.66 −7846.67 −7847.54 −7847.63 Eρτ 1327.84 1329.55 1329.79 1329.98 1330.20 1330.22 Eρ∆ρ 314.05 315.12 315.12 315.17 315.29 315.29 Eρ∇J −96.30 −96.44 −96.42 −96.43 −96.45 −96.45 ECoul 796.26 796.55 796.56 796.60 796.63 796.63 Table 5.3: Energies for 208 Pb from HFBTHO (computed with different numbers of HO shells N ), HFBFFT, and Sky1D. All energies are in MeV. The digits which do not coincide with HFBFFT are marked in bold. We find about 50 keV energy difference between HFBTHO and HFBFFT; this differ- ence can be primarily traced back to Ekin and Eρρ . As discussed in Refs. [159, 160], the kinetic energy converges slowly in the HO basis. To investigate this effect, we calculate 208 Pb using different numbers of HO shells in HFBTHO. We see in Table 5.3 that when we increase the number of HO shells to 30, the HFBTHO energies approach the HFBFFT values. It is also seen that Ekin and Eρρ exhibit the largest variations with N . In Refs. [159,161], the correction to the g.s. energy due to the finite number of HO shells 70 N has been derived: EL = E∞ + a0 e−2k∞ L , (5.26) p where L ≡ 2(N + 3/2 + 2)b, b is the oscillator length of our HO basis, and a0 , k∞ and E∞ are fit parameters. Then E∞ is the energy in the limit of infinitely large model space. The fit of Etot to Eq. (5.26) is presented in Fig. 5.1 and the resulting value of E∞ = −1635.786 MeV agrees fairly well with the HFBFFT and Sky1D values. Hence, obtaining an accurate kinetic as well as total energies in a HO basis-expansion solver requires a huge number of shells. In this context, the use of the coordinate-space representation is beneficial. -1634.0 -1634.5 Etot (MeV) -1635.0 -1635.5 E1 = ¡ 1635:786 MeV -1636.0 14 16 18 20 L (fm) Figure 5.1: Etot as a function of L for 208 Pb. The HFBTHO results are marked by red dots. The blue curve is fitted according to Eq. (5.26). 71 120 Sn HFBTHO HFBFFT Sky1D Sky2D HFBFFT Sky1D Sky2D n -renorm. Ẽkin ∆n -renorm. Etot −1018.77 −1018.34 −1018.45 −1018.37 −1018.78 −1018.92 −1018.74 ECoul 347.37 347.44 347.45 347.41 347.47 347.49 347.45 n Ekin 1340.51 1335.40 1335.18 1335.43 1339.17 1339.14 1338.72 p Ekin 830.75 830.97 831.01 830.01 831.25 831.31 831.28 n Epair −12.48 −7.37 −7.15 −7.40 −9.29 −9.14 −9.02 n Ẽkin 1328.03 1328.03 1328.03 1328.03 1329.88 1330.01 1329.70 ∆n 1.25 1.08 1.07 1.09 1.25 1.25 1.25 ϵF,n −8.02 −8.01 −8.01 −8.04 −8.00 −8.00 −8.04 Vpair,n −284.57 −342.70 −346.50 −354.90 −361.80 −367.30 −372.35 rrms 4.67 4.67 4.67 4.67 4.67 4.67 4.67 Table 5.4: Results of HFB + SLy4 calculations for 120 Sn using HFBTHO, HFBFFT, Sky1D, and Sky2D. Two neutron pairing renormalization variants are considered, by ad- justing the neutron pairing strengths in HFBFFT, Sky1D, and Sky2D to reproduce the HFBTHO values of Ẽkin n and ∆n . All energies are in MeV. The radius r rms is in fm. The digits which do not coincide with HFBFFT are marked in bold. 120 5.3.5 Spherical superfluid nucleus: Sn We now calculate 120 Sn which has a non-vanishing neutron pairing. The neutron pairing strength Vpair,n in HFBTHO is adjusted to the average experimental neutron pairing gap ∆n = 1.25 MeV. In HFBFFT, Sky1D and Sky2D, two pairing renormalizations are used. In the first variant, the neutron pairing strengths are adjusted to reproduce the HFBTHO n . In the second variant, the HFBTHO value of ∆n is matched. The results for value of Ẽkin both variants are displayed in Table 5.4. The neutron pairing strengths vary between the solvers, reflecting different structure of their quasiparticle pairing spaces, i.e., different pairing cutoff procedures and different structure of the discretized one-quasiparticle continuum. Although there are large discrepancies in Ekin n and E n pair between HFBFFT and HF- BTHO, in the first renormalization variant, the difference of the total energy, about 0.4 MeV, is quite reasonable considering the fact that the pairing space is treated differently and the HFBTHO results are affected by the basis truncation error. The difference in Etot 72 between the three coordinate-space solvers, less than 150 keV, reflects the dependence of the level density of the discretized quasiparticle continuum on the box boundary conditions assumed. In the pairing gap renormalization variant, the agreement of Etot is even better, with only 10-30 keV difference between HFBFFT, HFBTHO and Sky2D. In this variant, the magnitudes of the neutron pairing energy and kinetic energy are considerably larger as compared to the variant in which Ẽkin n is renormalized. Still, as seen in Table 5.4, both pairing renormalizations work reasonably well for 120 Sn. It is interesting to note that the total root-mean-square (rms) radii rrms are predicted very robustly in all renormalization variants. 102,110 5.3.6 Axially deformed nuclei: Zr The neutron-rich nuclei 102,110 Zr are suitable test cases, as they are known/expected to have large prolate deformations. In addition, 110 Zr is weakly bound, with the neutron chemical potential ϵF,n ≈ −3.5 MeV. The HFB proton pairing vanishes in this nucleus. In Table 5.5, we show results for 110 Zr with the two pairing renormalization schemes investigated in Sec. 5.3.5. It is seen that the HFBFFT results for various observables, i.e., total energy, quadrupole moments, and the rms radius, all agree well with those from HFBTHO in both pairing variants. In the case of 102 Zr, one also needs to consider proton pairing. In this case, we renor- malize both neutron and proton spectral pairing gaps by reproducing their values obtained from HFBTHO. In the calculation, 25 HO shells for both neutrons and protons are em- ployed in HFBTHO, which means that the s.p. proton and neutron spaces are the same. In HFBFFT, the canonical spaces are different as Ωn = 176, Ωp = 126. However, the actual 73 110 Zr HFBTHO HFBFFT Sky2D HFBFFT Sky2D n -renorm. Ẽkin n ∆ -renorm. Etot −893.97 −894.33 −894.32 −894.01 −894.01 ECoul 226.72 226.72 226.71 226.74 226.70 n Ekin 1368.08 1369.22 1368.98 1367.86 1367.13 p Ekin 632.03 632.05 632.13 632.16 632.05 n Epair −3.18 −4.31 −4.08 −2.30 −2.19 n Ẽkin 1364.90 1364.90 1364.90 1365.56 1364.94 ∆n 0.64 0.93 0.92 0.64 0.64 ϵF,n −3.55 −3.50 −3.52 −3.55 −3.57 Vpair,n −284.57 −409.80 −428.00 −371.00 −384.80 rrms 4.73 4.73 4.74 4.73 4.74 Qn20 789 794 795 791 796 p Q20 444 447 447 445 447 Table 5.5: Results of HFB + SLy4 calculations for 110 Zr with HFBTHO, HFBFFT and Sky2D. Two neutron pairing renormalization variants are considered, by adjusting the neu- tron pairing strengths in HFBFFT and Sky2D to reproduce the HFBTHO values of Ẽkin n p,n and ∆n . All energies are in MeV. The radius rrms is in fm and quadrupole moments Q20 are in fm2 . The HFB proton pairing vanishes in this nucleus. The digits which do not coincide with HFBFFT are marked in bold. pairing space is set by the soft-cutoff factor wα . It is seen in Table 5.6 that the benchmarking results following the pairing renormalization are very satisfactory. In particular, the results of HFBFFT, HFBTHO, and Sky2D are fairly close for the observables: Etot , rrms , and quadrupole moments. 240 5.3.7 Superdeformed heavy nucleus: Pu Compared with the HO basis, the coordinate-space representation can better capture strongly deformed configurations, such as the superdeformed fission isomer (f.i.) of 240 Pu. Indeed, very large configuration spaces are needed to guarantee the convergence of the HO expansion at large deformations [162, 163]. Given the large number of nucleons in 240 Pu, one needs to carefully consider the number of canonical states in HFBFFT and Sky2D calculations. To 74 102 Zr HFBTHO HFBFFT Sky2D Etot −859.65 −859.69 −859.67 ECoul 231.11 231.16 231.14 n Ekin 1202.02 1200.96 1201.97 p Ekin 651.25 651.22 651.27 n Epair −3.39 −2.50 −2.39 p Epair −1.97 −1.42 −1.38 n Ẽkin 1198.63 1199.53 1199.58 p Ẽkin 649.28 649.79 649.89 ∆n 0.69 0.69 0.69 ∆p 0.56 0.56 0.56 ϵF,n −5.43 −5.42 −5.44 ϵF,p −12.09 −12.09 −12.10 n Vpair −284.57 −367.00 −378.40 p Vpair −284.57 −372.00 −384.70 rrms 4.58 4.58 4.58 Qn20 639 639 640 p Q20 411 411 411 Table 5.6: Results of HFB + SLy4 calculations for 102 Zr using HFBTHO, HFBFFT and Sky2D. The pairing renormalization is carried out by adjusting the proton and neutron pairing strengths in HFBFFT and Sky2D to reproduce the HFBTHO values of ∆n and p,n ∆p . All energies are in MeV. The radius rrms is in fm and quadrupole moments Q20 are in fm2 . The digits which do not coincide with HFBFFT are marked in bold. this end, we performed a series of calculations by increasing the canonical space until the con- vergence had been reached. This has been done separately for the ground state (g.s.) and f.i. of 240 Pu. The final values are: (Ωn , Ωp ) = (300, 200) for the g.s. and (Ωn , Ωp ) = (400, 300) for the f.i. calculations. We renormalize the pairing strengths for 240 Pu to reproduce the g.s. ∆n and ∆p obtained in HFBTHO. The results are displayed in Table 5.7. In HFBTHO and Sky2D, the f.i. is found by performing quadrupole-moment con- strained calculations. The f.i. configuration in HFBFFT was computed by initializing the code with various HO deformations. As seen in Table 5.7, HFBTHO and HFBFFT results are very similar for g.s. energies, g.s. quadruple deformations, and radii. To test the functionality of HFBFFT for the f.i., we renormalize pairing strengths in 75 HFBFFT and Sky2D to the HFBTHO pairing gaps. Both coordinate-space solvers give very close results for the f.i., and they agree nicely with the HFBTHO results, see Table 5.7. Overall, the 240 Pu results obtained with HFBFFT for both the g.s. and f.i. show reasonable agreement with those from HFBTHO. ground state fission isomer 240 Pu HFBTHO HFBFFT HFBTHO HFBFFT Sky2D Etot −1802.11 −1802.43 −1797.00 −1797.35 −1797.35 ECoul 989.61 956.98 957.02 956.96 956.90 n Ekin 2938.92 2939.94 2922.56 2923.45 2923.43 p Ekin 1520.95 1521.46 1525.25 1525.52 1525.33 n Epair −3.11 −2.30 −3.52 −2.60 −2.48 p Epair −1.54 −1.22 −2.85 −2.19 −2.07 n Ẽkin 2935.81 2937.64 2919.03 2920.85 2920.55 p Ẽkin 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 ϵF,n −5.71 −5.70 −5.66 −5.65 −5.67 ϵF,p −5.69 −5.70 −5.76 −5.77 −5.79 rrms 5.93 5.93 6.40 6.40 6.40 Qn20 1784 1782 5063 5072 5071 p Q20 1166 1165 3336 3344 3343 n Vpair −284.57 −360.00 −284.57 −369.00 −384.60 p Vpair −284.57 −355.00 −284.57 −360.00 −375.80 s.p. space 25 shells (300, 200) 25 shells (400, 300) (400, 300) Table 5.7: Results of HFB + SLy4 calculations for 240 Pu ground state and fission isomer using HFBTHO, HFBFFT and Sky2D. The pairing strengths in HFBFFT and Sky2D were adjusted to reproduce the spectral pairing gaps obtained in HFBTHO for the g.s. and f.i. separately. The s.p. space for HFBFFT is defined by means of (Ωn , Ωp ). All energies are p,n in MeV, rrms is in fm, and Q20 are in fm2 . The digits which do not coincide with HFBFFT are marked in bold. 76 5.4 Summary: Development of Three-dimensional Skyrme HFB solver In this chapter, we developed a 3D Skyrme HFB solver HFBFFT in the coordinate-space representation using the canonical basis approach. The code is based on the well-optimized Sky3D solver. In HFBFFT, we implemented several new elements to facilitate calculations, namely (i) the sub-iteration method in configuration space to accelerate the convergence; (ii) the soft pairing cutoff and pairing annealing to avoid pairing breakdown; and (iii) a new algorithm to restore the Hermiticity of the HFB matrix. The new solver has been benchmarked for several spherical and deformed nuclei against HFBTHO, Sky2D, and (for spherical systems) Sky1D. The representation of the positive- energy continuum differs between HFB codes: In particular, it depends on the code’s geom- etry (spherical, cylindrical, Cartesian), the size of s.p. configuration space (number of HO shells, box size, grid size), and the effective pairing space. Consequently, even if the EDFs employed in two codes are identical, the pairing channel is usually described differently. This creates problems when comparing different HFB solvers as the perfect benchmarking is practically impossible [75]. In this work, we carried our careful inter-code comparisons by renormalizing pairing strengths to the spectral pairing gaps and/or the effective kinetic energy Ẽkin . While both methods give similar results, spectral pairing gaps are less sensitive to the s.p. space assumed. By carrying out calculations with different HFB solvers, we were able to assess the ranges of different uncertainties. For the total energy, the typical errors are: several keV due to the Hermiticity breaking; 10-80 keV due to different box boundary conditions assumed; 10-140 keV due to different quasiparticle continuum discretizations; and several hundred keV 77 due to the basis truncation in HO basis-expansion solvers. 78 Chapter 6 Conclusions and Outlook This dissertation consists of two parts: the applications of nuclear DFT to study nuclear reflection-asymmetric deformations and collective rotations, and the development of a 3D HFB solver. We first studied the microscopic origins of reflection-asymmetric deformations. In this project, we performed Skryme HFB calculations for octupole deformed Ba, Ra, and U iso- topes, as well as quadrupole deformed Yb isotopes. We applied the density multipole expan- sion method to decompose the total energy to different multipolarity contributions. Then, we analyzed their roles in the g.s. reflection-asymmetric deformation in both isospin and neutron-proton schemes. As a result, we concluded that the reflection-asymmetric shapes are driven by the odd-multipolarity isoscalar part of the nuclear interaction energy. To gain insights into the shell effects, we found that reflection-asymmetric shapes are favorable just above magic numbers where ∆ℓ = ∆j = 3 couplings are enhanced from plots of the single- particle levels. This study showed that the octupole deformation energy is subtle. Therefore, the dotriacontapole and higher-order components can be instrumental in the final octupole stability. The high-order effects could be explored in future studies. Next, we used the NLF to study the nuclear response to fast collective rotation. The CHO model was solved to give insight into the results of the self-consistent CHF calculations. In this work, we defined the simplified localization measure C τ by focusing on the interior 79 nodal structure and neglecting surface effects. We demonstrated that C τ is an excellent indicator of nuclear rotation. The concept of the NLF was generalized to anisotropic, spin- unsaturated, and spin-polarized systems. In the end, we suggested some applications of the NLF in the future, such as visualization of nuclear rotational and vibrational modes and time-dependent processes [33, 141, 164–169]. Finally, we built a reliable and efficient 3D Skyrme HFB solver HFBFFT in the coordinate-space representation using the canonical basis approach. We implemented sub- iteration, soft pairing cutoff, and Hermiticity restoration to speed up and stabilize the code. We benchmarked it against HFBTHO, Sky1D and Sky2D with different pairing renormal- izations. By comparing the results of several spherical and deformed systems, we confirmed the functionality of HFBFFT and assessed the ranges of different uncertainties to the total energy. To make HFBFFT more versatile, several enhancements are planned in the future. Most importantly, we intend to implement pairing regularization [44, 143, 151] to get rid of the dependence of pairing strengths on the cutoff energy. Another significant development is to compute potential energy surfaces defined by constraining one-body operators. This will enable us to use HFBFFT in the calculations of large-amplitude nuclear collective motion such as fission or fusion, for which the solvers based on the basis-expansion approach require the use of enormous configuration spaces. Finally, the performance of HFBFFT needs to be further optimized for modern supercomputer architectures. We hope HFBFFT become a powerful tool for exotic nuclei research. 80 BIBLIOGRAPHY 81 BIBLIOGRAPHY [1] Henri Becquerel. Sur les radiations émises par phosphorescence. Comptes rendus de 1’Academie des Sciences, Paris, 122:420–421, 1896. [2] E. Rutherford. The scattering of alpha and beta particles by matter and the structure of the atom. Phil. Mag. Ser. 6, 21:669–688, 1911. [3] James Chadwick. The existence of a neutron. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 136(830):692–708, 1932. [4] CF v Weizsäcker. Zur theorie der kernmassen. Zeitschrift für Physik, 96(7-8):431–458, 1935. [5] Petr Navrátil, Sofia Quaglioni, Guillaume Hupin, Carolina Romero-Redondo, and An- gelo Calci. Unified ab initio approaches to nuclear structure and reactions. Physica Scripta, 91(5):053002, apr 2016. [6] R. Machleidt and D.R. Entem. Chiral effective field theory and nuclear forces. Physics Reports, 503(1):1–75, 2011. [7] Bruce R. Barrett, Petr Navrátil, and James P. Vary. Ab initio no core shell model. Progress in Particle and Nuclear Physics, 69:131–181, 2013. [8] Maria Goeppert Mayer. On closed shells in nuclei. ii. Phys. Rev., 75:1969–1970, Jun 1949. [9] Otto Haxel, J. Hans D. Jensen, and Hans E. Suess. On the "magic numbers" in nuclear structure. Phys. Rev., 75:1766–1766, Jun 1949. [10] 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:427–488, Jun 2005. [11] Michael Bender, Paul-Henri Heenen, and Paul-Gerhard Reinhard. Self-consistent mean-field models for nuclear structure. Rev. Mod. Phys., 75:121–180, Jan 2003. [12] Hermann Arthur Jahn and Edward Teller. Stability of polyatomic molecules in de- generate electronic states-i—orbital degeneracy. Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences, 161(905):220–235, 1937. [13] P.-G. Reinhard and E. W. Otten. Transition to deformed shapes as a nuclear Jahn- Teller effect. Nucl. Phys. A, 420(2):173, 1984. 82 [14] H Schüler and Th Schmidt. Über abweichungen des atomkerns von der kugelsymmetrie. Zeitschrift für Physik, 94(7-8):457–468, 1935. [15] H. B. G. Casimir. On the Interaction between Atomic Nuclei and Electrons, pages 201–283. Springer Netherlands, Dordrecht, 1936. [16] W. Nazarewicz. Microscopic origin of nuclear deformations. Nucl. Phys. A, 574:27 – 49, 1994. [17] T.R. Werner, J. Dobaczewski, M.W. Guidry, W. Nazarewicz, and J.A. Sheikh. Micro- scopic aspects of nuclear deformation. Nucl. Phys. A, 578(1):1 – 30, 1994. [18] J. Dobaczewski, W. Nazarewicz, J. Skalski, and T. Werner. Nuclear deformation: A proton-neutron effect? Phys. Rev. Lett., 60:2254–2257, May 1988. [19] P. A. Butler and W. Nazarewicz. Intrinsic reflection asymmetry in atomic nuclei. Rev. Mod. Phys., 68(2):349–421, Apr 1996. [20] P. A. Butler. Pear-shaped atomic nuclei. Proc. R. Soc. A, 476(2239):20200202, 2020. [21] Yuchen Cao, S. E. Agbemava, A. V. Afanasjev, W. Nazarewicz, and E. Olsen. Land- scape of pear-shaped even-even nuclei. Phys. Rev. C, 102:024311, Aug 2020. [22] nuclear time-reversal violation and the schiff moment of 225 Ra. [23] Raj Kumar, J. A. Lay, and A. Vitturi. Nuclear fusion as a probe for octupole defor- mation in 224 Ra. Phys. Rev. C, 92:054604, Nov 2015. [24] Mengzhi Chen, Tong Li, Jacek Dobaczewski, and Witold Nazarewicz. Microscopic origin of reflection-asymmetric nuclear shapes. Phys. Rev. C, 103:034303, Mar 2021. [25] A. Bohr and B. R. Mottelson. Nuclear Structure, vol. II. W. A. Benjamin, Reading, 1975. [26] Zdzisĺaw Szymański. Fast Nuclear Rotation. Clarendon Press, January 1983. [27] Witold Nazarewicz. The nuclear collective motion. In José M. Arias and Manuel Lozano, editors, An Advanced Course in Modern Nuclear Physics, pages 102–140. Springer Berlin Heidelberg, Berlin, Heidelberg, 2001. [28] Stefan Frauendorf. Spontaneous symmetry breaking in rotating nuclei. Rev. Mod. Phys., 73(2):463–514, June 2001. [29] 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. 83 [30] Bastian Schuetrumpf, Chunli Zhang, and Witold Nazarewicz. Clustering and pasta phases in nuclear density functional theory. In W.-U. Schröder, editor, Nuclear Particle Correlations and Cluster Physics, pages 135–153. World Scientific, February 2017. [31] J-P Ebran, E Khan, T Nikšić, and D Vretenar. Localization and clustering in atomic nuclei. J. Phys. G, 44(10):103001, aug 2017. [32] B. Schuetrumpf and W. Nazarewicz. Cluster formation in precompound nuclei in the time-dependent framework. Phys. Rev. C, 96(6):064608, December 2017. [33] A. S. Umar, C. Simenel, and K. Godbey. Pauli energy contribution to the nucleus- nucleus interaction. Phys. Rev. C, 104:034619, Sep 2021. [34] C. L. Zhang, B. Schuetrumpf, and W. Nazarewicz. Nucleon localization and fragment formation in nuclear fission. Phys. Rev. C, 94(6):064323, December 2016. [35] Guillaume Scamps and Cédric Simenel. Impact of pear-shaped fission fragments on mass-asymmetric fission in actinides. Nature, 564(7736):382–385, 2018. [36] Guillaume Scamps and Cédric Simenel. Effect of shell structure on the fission of sub- lead nuclei. Phys. Rev. C, 100:041602, Oct 2019. [37] Zachary Matheson, Samuel A. Giuliani, Witold Nazarewicz, Jhilam Sadhukhan, and Nicolas Schunck. Cluster radioactivity of 294 118 Og176 . Phys. Rev. C, 99:041304, Apr 2019. [38] Jhilam Sadhukhan, Samuel A. Giuliani, Zachary Matheson, and Witold Nazarewicz. Efficient method for estimation of fission fragment yields of r-process nuclei. Phys. Rev. C, 101:065803, Jun 2020. [39] T. Li, M. Z. Chen, C. L. Zhang, W. Nazarewicz, and M. Kortelainen. Nucleon local- ization function in rotating nuclei. Phys. Rev. C, 102:044305, Oct 2020. [40] 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, 1984. [41] 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:783, 1987. [42] 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:2809–2840, Jun 1996. [43] J. Dobaczewski and W. Nazarewicz. Hartree-Fock-Bogoliubov solution of the pairing Hamiltonian in finite nuclei. In R. A. Broglia and V. Zelevinsky, editors, 50 Years of Nuclear BCS, page 40. World Scientific, 2013. 84 [44] A. Bulgac. Local density approximation for systems with pairing correlations. Phys. Rev. C, 65(5):051305, 2002. [45] 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. [46] Baishan Hu, Weiguang Jiang, Takayuki Miyagi, Zhonghao Sun, Andreas Ekström, Christian Forssén, Gaute Hagen, Jason D. Holt, Thomas Papenbrock, S. Ragnar Stroberg, and Ian Vernon. Ab initio predictions link the neutron skin of 208 pb to nuclear forces. 2021. [47] P. Hohenberg and W. Kohn. Inhomogeneous electron gas. Phys. Rev., 136:B864–B871, Nov 1964. [48] T.H.R. Skyrme. The effective nuclear potential. Nucl. Phys., 9(4):615 – 634, 1958. [49] JF Berger, M Girod, and D Gogny. Microscopic analysis of collective dynamics in low energy fission. Nucl. Phys. A, 428:23–36, 1984. [50] Stanislas Goriely, Stéphane Hilaire, Michel Girod, and S Péru. First Gogny-Hartree- Fock-Bogoliubov nuclear mass model. Phys. Rev. Lett, 102(24):242501, 2009. [51] F Chappert, M Girod, and S Hilaire. Towards a new Gogny force parameterization: Impact of the neutron matter equation of state. Phys. Lett. B, 668(5):420–424, 2008. [52] GA Lalazissis, Tamara Nikšić, Dario Vretenar, and Peter Ring. New relativistic mean-field interaction with density-dependent meson-nucleon couplings. Phys. Rev. C, 71(2):024312, 2005. [53] Tamara Nikšić, Dario Vretenar, and Peter Ring. Relativistic nuclear energy density functionals: Adjusting parameters to binding energies. Phys. Rev. C, 78(3):034318, 2008. [54] Peng-Wei Zhao, Zhi-Pan Li, Jiang-Ming Yao, J Meng, et al. New parametrization for the nuclear covariant energy density functional with a point-coupling interaction. Phys. Rev. C, 82(5):054319, 2010. [55] GA Lalazissis, S Karatzikos, R Fossion, D Pena Arteaga, AV Afanasjev, and P Ring. The effective force NL3 revisited. Phys. Lett. B, 671(1):36–41, 2009. [56] Nicolas Schunck, editor. Energy Density Functional Methods for Atomic Nuclei. 2053- 2563. IOP Publishing, 2019. [57] 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, Jan 2011. 85 [58] J. Erler, N. Birge, M. Kortelainen, W. Nazarewicz, E. Olsen, A.M. Perhac, and M. Stoitsov. The limits of the nuclear landscape. Nature, 486:509, 2012. [59] Jacek Dobaczewski. Current developments in nuclear density functional methods. J. Phys. Conf. Ser., 312(9):092002, sep 2011. [60] 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, 2001. [61] 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:034310, Mar 2009. [62] Peter Ring and Peter Schuck. The nuclear many-body problem. Springer-Verlag, Berlin, 1980. [63] Claude Bloch and Albert Messiah. The canonical form of an antisymmetric tensor and its application to the theory of superconductivity. Nuclear Physics, 39:95–106, 1962. [64] P-G Reinhard, Michaėl Bender, K Rutz, and Joachim Alexander Maruhn. An HFB scheme in natural orbitals. Z. Phys., A Hadrons nucl., 358(3):277–278, 1997. [65] Naoki Tajima. Canonical-basis solution of the Hartree-Fock-Bogoliubov equation on a three-dimensional Cartesian mesh. Phys. Rev. C, 69:034305, Mar 2004. [66] Peter Ring and Peter Schuck. The nuclear many-body problem. Springer-Verlag, Berlins, 2004. [67] Michael Bender, Rémi Bernard, George Bertsch, Satoshi Chiba, Jacek Dobaczewski, Noël Dubray, Samuel A Giuliani, Kouichi Hagino, Denis Lacroix, Zhipan Li, Piotr Magierski, Joachim Maruhn, Witold Nazarewicz, Junchen Pei, Sophie Péru, Nathalie Pillet, Jørgen Randrup, David Regnier, Paul-Gerhard Reinhard, Luis M Robledo, Wouter Ryssens, Jhilam Sadhukhan, Guillaume Scamps, Nicolas Schunck, Cédric Simenel, Janusz Skalski, Ionel Stetcu, Paul Stevenson, Sait Umar, Marc Verriere, Dario Vretenar, Michał Warda, and Sven Åberg. Future of nuclear fission theory. J. Phys. G, 47(11):113002, oct 2020. [68] 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. Comm., 167(1):43 – 63, 2005. [69] R. Navarro Perez, N. Schunck, R.-D. Lasseri, C. Zhang, and J. Sarich. Axially deformed solution of the Skyrme–Hartree–Fock–Bogolyubov equations using the transformed 86 harmonic oscillator basis (III) hfbtho (v3.00): A new version of the program. Comput. Phys. Commun., 220:363–375, 2017. [70] 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, 1997. [71] 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, 2017. [72] Jacek Jan Dobaczewski, Paweł Bączyk, Pierre Becker, Michael Bender, Karim Ben- naceur, Jeremy Bonnard, Yuan Gao, Andrea Idini, Maciej Konieczka, Markus Korte- lainen, and et al. Solution of universal nonrelativistic nuclear DFT equations in the Cartesian deformed harmonic-oscillator basis. (IX) HFODD (3.06h): a new version of the program. J. Phys. G, 2021. [73] M. Stoitsov, N. Michel, and K. Matsuyanagi. New efficient method for perform- ing Hartree-Fock-Bogoliubov calculations for weakly bound nuclei. Phys. Rev. C, 77:054301, May 2008. [74] K. Bennaceur and J. Dobaczewski. Coordinate-space solution of the Skyrme–Hartree–Fock–Bogolyubov equations within spherical symmetry. the program HFBRAD (v1.00). Comput. Phys. Commun., 168(2):96 – 122, 2005. [75] 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:064306, Dec 2008. [76] P.-G. Reinhard, B. Schuetrumpf, and J.A. Maruhn. The Axial Hartree–Fock + BCS Code SkyAx. Comput. Phys. Commun., 258:107603, 2021. [77] J.A. Maruhn, P.-G. Reinhard, P.D. Stevenson, and A.S. Umar. The TDHF code Sky3D. Comput. Phys. Commun., 185(7):2195 – 2216, 2014. [78] 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, 2018. [79] P.-G. Reinhard. Skyrme-Hartree-Fock calculations of the nuclear ground state. In K. Langanke, S.E. Koonin, and J.A. Maruhn, editors, Computational Nuclear Physics I - Nuclear Structure, page 28. Springer, Berlin, 1991. [80] P.-G. Reinhard, “HFB solvers Sky1D and Sky2D”, unpublished, 2021. 87 [81] P. Bonche, H. Flocard, and P.H. Heenen. Solution of the Skyrme HF+BCS equation on a 3D mesh. Comput. Phys. Commun., 171(1):49–62, 2005. [82] W. Ryssens, V. Hellemans, M. Bender, and P.-H. Heenen. Solution of the Skyrme–HF+BCS equation on a 3D mesh, II: a new version of the Ev8 code. Comput. Phys. Commun., 187:175–194, 2015. [83] Ryssens, W., Bender, M., and Heenen, P. -H. Iterative approaches to the self-consistent nuclear energy density functional problem - Heavy ball dynamics and potential pre- conditioning. Eur. Phys. J. A, 55(6):93, 2019. [84] Guillaume Scamps, Stephane Goriely, Erik Olsen Olsen, Michael Bender, and Wouter Ryssens. Skyrme-Hartree-Fock-Bogoliubov mass models on a 3D mesh: effect of triaxial shape. Eur. Phys. J. A, 57(12):333, Dec 2021. [85] J. C. Pei, G. I. Fann, R. J. Harrison, W. Nazarewicz, Yue Shi, and S. Thornton. Adaptive multi-resolution 3D Hartree-Fock-Bogoliubov solver for nuclear structure. Phys. Rev. C, 90:024317, Aug 2014. [86] Shi Jin, Kenneth J. Roche, Ionel Stetcu, Ibrahim Abdurrahman, and Aurel Bulgac. The LISE package: Solvers for static and time-dependent superfluid local density ap- proximation equations in three dimensions. Comput. Phys. Commun., 269:108130, 2021. [87] Shi Jin, Aurel Bulgac, Kenneth Roche, and Gabriel Wlazłowski. Coordinate-space solver for superfluid many-fermion systems with the shifted conjugate-orthogonal conjugate-gradient method. Phys. Rev. C, 95:044302, Apr 2017. [88] Yu Kashiwaba and Takashi Nakatsukasa. Coordinate-space solver for finite- temperature Hartree-Fock-Bogoliubov calculations using the shifted Krylov method. Phys. Rev. C, 101:045804, Apr 2020. [89] S. E. Agbemava, A. V. Afanasjev, and P. Ring. Octupole deformation in the ground states of even-even nuclei: A global analysis within the covariant density functional theory. Phys. Rev. C, 93:044304, Apr 2016. [90] S. E. Agbemava and A. V. Afanasjev. Octupole deformation in the ground states of even-even Z ∼ 96, N ∼ 196 actinides and superheavy nuclei. Phys. Rev. C, 96:024301, Aug 2017. [91] Zhong Xu and Zhi-Pan Li. Microscopic analysis of octupole shape transitions in neutron-rich actinides with relativistic energy density functional. Chin. Phys. C, 41(12):124107, 2017. 88 [92] William D. Myers and Wladyslaw J. Swiatecki. Nuclear masses and deformations. Nucl. Phys., 81(1):1 – 60, 1966. [93] P. Möller, R. Bengtsson, B. G. Carlsson, P. Olivius, T. Ichikawa, H. Sagawa, and A. Iwamoto. Axial and reflection asymmetry of the nuclear ground state. At. Data Nucl. Data Tables, 94(5):758–780, Sep 2008. [94] J.L. Egido and L.M. Robledo. Parity-projected calculations on octupole deformed nuclei. Nucl. Phys. A, 524(1):65 – 87, 1991. [95] L. M. Robledo. Enhancement of octupole strength in near spherical nuclei. Eur. Phys. J. A, 52(9):300, 2016. [96] S. Y. Xia, H. Tao, Y. Lu, Z. P. Li, T. Nikšić, and D. Vretenar. Spectroscopy of reflection-asymmetric nuclei with relativistic energy density functionals. Phys. Rev. C, 96:054303, Nov 2017. [97] Y. Fu, H. Wang, L.-J. Wang, and J. M. Yao. Odd-even parity splittings and octupole correlations in neutron-rich ba isotopes. Phys. Rev. C, 97:024338, Feb 2018. [98] V. M. Strutinsky. Remarks about pear-shaped nuclei. Physica, 22(6):1166–1167, Jan 1956. [99] Kiuck Lee and D. R. Inglis. Stability of pear-shaped nuclear deformations. Phys. Rev., 108:774–778, Nov 1957. [100] W. Nazarewicz, P. Olanders, I. Ragnarsson, J. Dudek, G.A. Leander, P. Möller, and E. Ruchowska. Analysis of octupole instability in medium-mass and heavy nuclei. Nucl. Phys. A, 429(2):269 – 295, 1984. [101] A. Bohr. The coupling of nuclear surface oscillations to the motion of individual nucleons. Kgl. Danske Videnskab. Selskab, Mat. fys. Medd, 26(14), 1952. [102] A. Bohr and B. R. Mottelson. Nuclear Structure, vol. II. W. A. Benjamin, Reading, 1975. [103] D.R. Bes, R.A. Broglia, and B.S. Nilsson. Microscopic description of isoscalar and isovector giant quadrupole resonances. Phys. Rep., 16(1):1 – 56, 1975. [104] D. Vautherin. Hartree-Fock calculations with Skyrme’s interaction. II. Axially de- formed nuclei. Phys. Rev. C, 7:296–316, Jan 1973. [105] 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:014316, Jan 2004. 89 [106] 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, 1975. [107] 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, 1998. [108] M. Kortelainen, J. McDonnell, W. Nazarewicz, E. Olsen, P.-G. Reinhard, J. Sarich, N. Schunck, S. M. Wild, D. Davesne, J. Erler, and A. Pastore. Nuclear energy density optimization: Shell structure. Phys. Rev. C, 89:054314, 2014. [109] P. Rozmej, S. Ćwiok, and A. Sobiczewski. Is octupole deformation sufficient to describe the properties of ‘octupolly’ unstable nuclei? Phys. Lett. B, 203(3):197 – 199, 1988. [110] Adam Sobiczewski, Zygmunt Patyk, Stefan Ćwiok, and Piotr Rozmej. Study of the potential energy of ‘octupole’-deformed nuclei in a multidimensional deformation space. Nucl. Phys. A, 485(1):16 – 30, 1988. [111] 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, 1989. [112] 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, 1990. [113] Stefan Ćwiok and Witold Nazarewicz. Ground-state shapes and spectroscopic proper- ties of Z ∼ 58, N ∼ 88 nuclei. Nucl. Phys. A, 496(2):367 – 384, 1989. [114] S. Ćwiok and W. Nazarewicz. Reflection-asymmetric shapes in transitional odd-A Th isotopes. Phys. Lett. B, 224(1):5 – 10, 1989. [115] S. Ćwiok and W. Nazarewicz. Reflection-asymmetric shapes in odd-A actinide nuclei. Nucl. Phys. A, 529(1):95 – 114, 1991. [116] W. Nazarewicz. Static multipole deformations in nuclei. Prog. Part. Nucl. Phys., 28:307 – 330, 1992. [117] R Bengtsson, J Dudek, W Nazarewicz, and P Olanders. A systematic comparison between the Nilsson and Woods-Saxon deformed shell model potentials. Phys. Scr., 39(2):196–220, feb 1989. [118] Takashi Nakatsukasa, Kenichi Matsuyanagi, Masayuki Matsuo, and Kazuhiro Yabana. Time-dependent density-functional description of nuclear dynamics. Rev. Mod. Phys., 88:045004, Nov 2016. 90 [119] 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. [120] Andreas Savin, Ove Jepsen, Jürgen Flad, Ole Krogh Andersen, Heinzwerner Preuss, and Hans Georg von Schnering. Electron Localization in Solid-State Structures of the Elements: the Diamond Structure. Angew. Chem. Int. Ed. Engl., 31(2):187–188, 1992. [121] Andreas Savin, Reinhard Nesper, Steffen Wengert, and Thomas F. Fässler. ELF: The Electron Localization Function. Angew. Chem. Int. Ed. Engl., 36(17):1808–1832, 1997. [122] T. Burnus, M. A. L. Marques, and E. K. U. Gross. Time-dependent electron localization function. Phys. Rev. A, 71(1):010501, January 2005. [123] Patricio Fuentealba, E. Chamorro, and Juan C. Santos. Chapter 5 Understanding and using the electron localization function. In Alejandro Toro-Labbé, editor, Comput. Theor. Chem., volume 19 of Theoretical Aspects of Chemical Reactivity, pages 57–85. Elsevier, January 2007. [124] Paul Jerabek, Bastian Schuetrumpf, Peter Schwerdtfeger, and Witold Nazarewicz. Electron and Nucleon Localization Functions of Oganesson: Approaching the Thomas- Fermi Limit. Phys. Rev. Lett., 120(5):053001, January 2018. [125] 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. [126] 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. [127] 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. [128] D. Glas, U. Mosel, and P. G. Zint. The cranked harmonic oscillator in coordinate space. Z Physik A, 285(1):83–87, March 1978. [129] T. Troudet and R. Arvieu. Shapes of nuclear configurations in a cranked harmonic oscillator model. Ann. Phys. (N.Y.), 134(1):1 – 44, 1981. [130] W. Nazarewicz and J. Dobaczewski. Dynamical symmetries, multiclustering, and oc- tupole susceptibility in superdeformed and hyperdeformed nuclei. Phys. Rev. Lett., 68:154–157, Jan 1992. [131] S. G. Nilsson and I. Ragnarsson. Shapes and Shells in Nuclear Structure. Cambridge University Press, Cambridge, 1995. 91 [132] Z. Szymański and W. Nazarewicz. Rotating pseudo-oscillator scheme: pseudo-spin symmetry and identical bands. Phys. Lett. B, 433(3):229 – 235, 1998. [133] Mark 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:1704– 1708, Nov 1976. [134] 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, 1978. [135] J. Kunz and U. Mosel. Flow patterns for collective rotations in heavy nuclei. Nucl. Phys. A, 323(2):271 – 284, 1979. [136] J. Fleckner, J. Kunz, U. Mosel, and E. Wuest. Current distributions for rotating nuclei. Nucl. Phys. A, 339(2):227 – 237, 1980. [137] M. Durand, P. Schuck, and J. Kunz. Semiclassical description of currents in normal and superfluid rotating nuclei. Nucl. Phys. A, 439(2):263 – 288, 1985. [138] H. Laftchiev, D. Samsœn, P. Quentin, and I. N. Mikhailov. Equivalence of pairing correlations and intrinsic vortical currents in rotating nuclei. Phys. Rev. C, 67:014307, Jan 2003. [139] 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:225 – 233, 2004. [140] 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:024329, Feb 2018. [141] Kristian Petrík and Markus Kortelainen. Thouless-Valatin rotational moment of inertia from linear response theory. Phys. Rev. C, 97(3):034321, March 2018. [142] N. Michel, K. Matsuyanagi, and M. Stoitsov. Gamow-Hartree-Fock-Bogoliubov method: Representation of quasiparticles with Berggren sets of wave functions. Phys. Rev. C, 78:044319, Oct 2008. [143] J. C. Pei, A. T. Kruppa, and W. Nazarewicz. Quasiparticle continuum and resonances in the Hartree-Fock-Bogoliubov theory. Phys. Rev. C, 84:024311, Aug 2011. [144] Md Afibuzzaman, Bastian Schuetrumpf, and Hasan Metin Aktulga. Scalable nuclear density functional theory with Sky3D. Comput. Phys. Commun., 223:34 – 44, 2018. 92 [145] 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, 1992. [146] R W Hockney. Potential calculation and some applications. Methods Comput. Phys., 9:135–211, 1 1970. [147] 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, 1979. [148] M. Frigo and S.G. Johnson. The design and implementation of FFTW3. Proc. IEEE, 93(2):216–231, 2005. [149] P.-G. Reinhard and R.Y. Cusson. A comparative study of Hartree-Fock iteration techniques. Nucl. Phys. A, 378(3):418 – 442, 1982. [150] C. Bottcher, M. R. Strayer, A. S. Umar, and P.-G. Reinhard. Damped relaxation techniques to calculate relativistic bound states. Phys. Rev. A, 40:4182–4189, Oct 1989. [151] 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:044319, Apr 2006. [152] P. Bonche, H. Flocard, P.H. 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, 1985. [153] S.J. Krieger, P. Bonche, H. Flocard, P. Quentin, and M.S. Weiss. An improved pair- ing interaction for mean field calculations using Skyrme potentials. Nucl. Phys. A, 517(2):275–284, 1990. [154] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes in C: The Art of Scientific Computing. Cambridge University Press, New York, 2 edition, 1992. [155] Steven G. Johnson. Notes on FFT-based differentiation. [156] T. Papenbrock and G. F. Bertsch. Pairing in low-density Fermi gases. Phys. Rev. C, 59:2052–2055, Apr 1999. [157] M. Bender, K. Rutz, P. Reinhard, and J. Maruhn. Pairing gaps from nuclear mean-field models. Eur. Phys. J. A, 8, 05 2000. 93 [158] B. Schuetrumpf, W. Nazarewicz, and P.-G. Reinhard. Time-dependent density func- tional theory with twist-averaged boundary conditions. Phys. Rev. C, 93:054304, May 2016. [159] R. J. Furnstahl, G. Hagen, and T. Papenbrock. Corrections to nuclear energies and radii in finite oscillator spaces. Phys. Rev. C, 86:031301, Sep 2012. [160] 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:044332, Apr 2016. [161] 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:044326, Apr 2013. [162] N. Nikolov, N. Schunck, W. Nazarewicz, M. Bender, and J. Pei. Surface symmetry energy of nuclear energy density functionals. Phys. Rev. C, 83:034305, Mar 2011. [163] N. Schunck. Density functional theory approach to nuclear fission. Acta Phys. Pol. B, 44(3):263, 2013. [164] T.S. Dumitrescu and Toru Suzuki. Spin distributions in the vibrational excitations of closed-shell nuclei. Nucl. Phys. A, 423(2):277 – 319, 1984. [165] Yoritaka Iwata and Joachim A. Maruhn. Enhanced spin-current tensor contribution in collision dynamics. Phys. Rev. C, 84:014616, Jul 2011. [166] Jun Xu and Bao-An Li. Probing in-medium spin–orbit interaction with intermediate- energy heavy-ion collisions. Phys. Lett. B, 724(4):346 – 350, 2013. [167] Yin Xia, Jun Xu, Bao-An Li, and Wen-Qing Shen. Spin-orbit coupling and the up- down differential transverse flow in intermediate-energy heavy-ion collisions. Phys. Rev. C, 89:064606, Jun 2014. [168] K. Godbey, A. S. Umar, and C. Simenel. Dependence of fusion on isospin dynamics. Phys. Rev. C, 95:011601, Jan 2017. [169] V. O. Nesterenko, A. Repko, J. Kvasil, and P.-G. Reinhard. Individual low-energy toroidal dipole state in 24 Mg. Phys. Rev. Lett., 120:182501, May 2018. 94