Improvements to the macroscopic-microscopic approach of nuclear fission
LLLNL-JRNL-813628
Improvements to the macroscopic-microscopic approach of nuclear fission
Marc Verriere ∗ Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USANuclear and Chemical Sciences Division, Lawrence Livermore National Laboratory, Livermore, California 94551, USA
Matthew Ryan Mumpower † Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA (Dated: August 18, 2020)The well established macroscopic-microscopic (mac-mic) description of nuclear fission enablesthe prediction of fission fragment yields for a broad range of fissioning systems. In this work,we present several key enhancements to this approach. We improve upon the microscopic sectorof nuclear potential energy surfaces by magnifying the Lipkin-Nogami equations’ resolution andstrengthening the Strutinsky procedure, thus reducing spurious effects from the continuum. Wefurther present a novel deterministic method for calculating fission dynamics under the assumptionof strongly damped nucleonic motion. Our technique utilizes the memoryless property of MarkovChains to produce fission yields that do not rely on the statistical accumulation of scission events.We show that our new technique is equivalent to the Metropolis random-walk pioneered over thepast decade by Randrup and colleagues. It further improves upon it, as we remove the need foraltering the nuclear landscape via a biased potential. With our final improvement, we calculatescission configurations using particle number projection, which affords the simultaneous calculationof both mass and charge yield distributions. Fission fragments are thus calculated from the quantummechanical A -body states of the potential energy surface rather than the collective mass asymmetryvariable ( α g ) of the Finite-Range Liquid-Drop Model (FRLDM) used in past work. We highlightthe success of our enhancements by predicting the odd-even staggering and the charge polarizationfor the neutron-induced fission of U and U. ∗ [email protected] † [email protected] a r X i v : . [ nu c l - t h ] A ug I. INTRODUCTION
The complexity of nuclear fission makes this reaction challenging to model theoretically from fundamental principles.Fission data, especially those based on fragment yields, are a key ingredient in many applications. For instance, theaccurate description of fission yields is influential for nuclear engineering and waste management [1, 2], the productionof radioactive isotopes via fragmentation [3–7], reactor neutrinos [8, 9], and in the pursuit of synthesizing superheavyelements [10–13]. Fission may also play an important role in the formation of the heavy elements in astrophysicalprocesses [14–19].There are many approaches to the theoretical description of fission. Fully microscopic models describe the fissionprocess by assuming the nuclear interaction between the nucleons only. The most frequently used implementationis the Energy Density Functional (EDF) theory, where the density-dependence of the energy is derived from aneffective interaction or directly parameterized. A quantum state is then extracted using a time-dependent mean-fieldapproximation, such as the Time-Dependent Hartree-Fock method (TDHF), with the possibility to include dissipativeeffects [20–28] or through the resolution of the Schr¨odinger equation in a collective space defined with a reduced numberof collective degrees of freedom (TDGCM), mainly the quadrupole and octupole multipolar moments [29–36]. Thelatter type of microscopic calculation can describe the fissioning system’s deformations before scission only recentlydue to lengthy computational costs. Statistical methodologies also provide a practical path to obtain a description offission and further enable large-scale calculations [37–39].Approaching the description of fission observables from macroscopic-microscopic theory provides yet another al-ternative. In this approach, the fissioning system is modeled to first approximation as a macroscopic system (e.g.,a liquid-drop or a droplet of nuclear matter) to obtain the smooth part of its energy [40, 41]. Corrections are thenapplied that account for the missing microscopic effects that may contribute to rapid variations in energy [42–44].The time-evolution of the fissioning system is then obtained statistically with the exact or approximate resolution ofthe Langevin equations [45–47]. Approximate methods that assume strong damping are modeled via a random walkon the potential energy surface and are found to produce a good agreement with known data [48, 49].The benefit of such semi-classical approaches over the microscopic approaches applied to the description of fissionis the inclusion of larger spaces for the collective shape degrees of freedom and the ability to model stochasticdynamics, whilst retaining reasonable calculation times. Thus, such models are applicable to a range of reactionsacross many heavy fissioning systems [50]. To date, these approaches have examined several starting configurationsfor the dynamics, including from the ground state, inner saddle region, and after the last saddle [47, 50]. Thosestudies which start after the last saddle explicitly ignore dynamics before this point, potentially neglecting importantpathways. Studies that start at the ground state may introduce a phenomenological tilting of the potential energysurface in an attempt to force the calculation to go over the fission barrier sufficiently quickly. The downside to thetilting is that it may introduce an unphysical bias to the results. Finally, the scission configurations of past work,e.g., [48], may be calculated with reference only to the collective asymmetry coordinate ( α g ), thereby limiting thepredictive power of the approach.In this work, we address several of the major shortcomings of past macroscopic-microscopic fission studies to describethe nascent fission fragment yields. In Sec. II A, we outline the calculation of the potential energy surface, which hasbeen improved by using an enhanced Strutinsky method that drastically reduces the continuum effects and allowsthe use of higher-quality particle basis. We present a new deterministic method to approximate fission dynamics inSec. II B and show that it has the capacity to reproduce the results of past random-walk methods. Our new methodhas comparable calculation time while avoiding the use of biasing the potential energy surface. Further, we improvethe estimation of the fission fragment mass and charge probability distribution based on a microscopic projectiontechnique [51, 52]. In Sec. III, we present the results we have obtained with our approach. II. THEORETICAL APPROACH
Our description of the fission process can be decomposed into three main steps. First, the shape-dependent potentialenergy surface (PES) of the corresponding fissioning system (typically target-plus-neutron) is determined using a semi-classical method based on the macroscopic Finite-Range Liquid-Drop Model (FRLDM) with microscopic corrections.An effective ground state and barrier height are then extracted from the PES, as presented in Sec. II A. In Sec. II B, wediscuss our new Random-Walk-based algorithm used to obtain the probability to populate each scission configurationof the PES. Finally, in Sec. II C, we present how we can deduce, for the first time in this type of approach, the fissionfragment probability distribution in charge and mass, Y ( Z, A ), based on a microscopic projection technique.
A. Potential energy landscape
Several parameterizations of the sharp macroscopic density have been developed, see e.g., Ref. [53], to specifythe relevant degrees of freedom associated with large deformations encountered in fission. In this work, we restrictourselves to the description of binary fission. Even though there is experimental evidence for ternary [54–57] andeven quaternary [58] fission, their contribution to the fragment probability distribution is smaller than binary fissionby orders of magnitude. Thus, we use the so-called Matched-Quadratic-Surface (MQS) parameterization that wasintroduced in Refs. [59–62] for the specification of our shape families. The MQS parameterization contains ninedegrees of freedom. Six degrees of freedom remain by ensuring a smooth junction between the bodies. Because one ofthe parameters corresponds to the center of mass, this may be set to the origin, reducing the number of parametersto five. We use the symbol q to denote a specific MQS shape. The lattice we use to calculate the PES is taken fromRef. [63]. We note that the nodes associated with α g = − .
02 are redundant due to parity-reversal symmetry, andtherefore are not explicitly calculated.For a nuclear system defined by Z protons, and N neutrons, the potential energy E ( q ) associated with a given setof MQS parameters can be written as E ( q ) = E mac ( q ) + ∆ E shell ( q ) + ∆ E pair . ( q ) , (1)where E mac ( q ) is the macroscopic energy, and the remaining terms define the microscopic corrections. It is obtainedassuming that the fissioning system is a nuclear drop of charged liquid. Implicitly, these terms also depend on Z and N as discussed in greater detail in Appendix A.In order to obtain the sharp contribution to the energy from shell effects, we determine a microscopic many-bodystate at the mean-field approximation. The effective averaged potential for the isospin τ is V ( τ ) ( q ) = V ( τ )1 ( q ) + V ( τ )C ( q ) + V ( τ )s . o . ( q ) , (2)where each term is taken from Ref. [62]. The first term, V ( τ )1 , corresponds to the parameterized mean-field associatedwith the central part of the nuclear interaction and is obtained assuming a Yukawa interaction between the nucleons, V ( τ )1 ( r ; q ) = − V τ πa (cid:90) V d r (cid:48) e −| r − r (cid:48) | /a pot | r − r (cid:48) | /a pot , (3)where V is the shape associated with the MQS parameters, q , scaled to have fixed volume, πAR ( R pot is definedby Eq. (81) of Ref. [64]). The potential depths V τ are given by V n = V s + V a ¯ δ (4) V p = V s − V a ¯ δ, (5)where V s and V a are parameters of the model and ¯ δ is given by Eq. (85) of Ref. [64]. The Coulomb term V ( τ )C ( r ; q ) isonly acting on protons and given by V (p)C ( r ; q ) = e Z πAr (cid:90) V d r | r − r (cid:48) | . (6)The spin-orbit term has the expression V ( τ )s . o . = λ τ (cid:18) (cid:126) m nuc c (cid:19) ∇ V ( τ )1 σ × p (cid:126) (7)where the interaction strength λ τ is taken for each isospin τ as a linear function of the mass A [64] such that λ τ = k τ A + l τ , (8)and k τ and l τ are parameters.Assuming that the particles of the compound system are independent, the many-body state can, therefore, beobtained as a Slater determinant of particles, the state of each particle being an eigenfunction of the Hamiltonianassociated with the energy e i of the particle. Differentiating ourselves from past FRLDM work, the shell correctionis calculated independently for each isospin by using the improved Strutinsky method presented in Ref. [65, 66].This procedure removes spurious contribution from the continuum, which happens when the number of shells inthe Harmonic Oscillator basis, N , is set too large, thus causing divergences in the calculation of the energy. Theavoidance of such spurious contributions is the reason why the size of the particle basis was limited to N = 12 inpast work. In this work, we use N = 20, which is sufficiently larger than N = 12 while remaining computationallymanageable. We have tested our Strutinksy procedure up to N = 30 and found no anomalies.Pairing correlations are obtained using the Lipkin-Nogami approach with the seniority-pairing approximation onthe Slater-determinant of particles used to estimate shell effects [67, 68]. The relevant equations read N pair = L min + L max (cid:88) k = L min v k (9)2 G = L max (cid:88) k = L min (cid:112) ( (cid:15) k − λ ) + ∆ (10) v k = 12 (cid:34) − (cid:15) k − λ (cid:112) ( (cid:15) k − λ ) + ∆ (cid:35) (11) (cid:15) k = e k + (4 λ − G ) v k (12) λ = G L max (cid:88) k = L min u k v k L max (cid:88) l = L min l (cid:54) = k u l v lL max (cid:88) k = L min u k v k L max (cid:88) l = L min l (cid:54) = k u l v l . (13)This is a nonlinear system of 2 N v + 3 equations, where N v = L max − L min + 1 is the number of pairs in the valencespace. The unknowns of this systems are the pairing gap ∆, the Fermi energy λ , the number-fluctuation constant λ and for k = L min . . . L max , v k are the occupation amplitudes and (cid:15) k are the shifted single-particle energies. Thissystem of equations is completely determined by the single-particle energies e k , the number of paired levels N pair , thefirst and last levels 0 ≤ L min ≤ L max in the valence space and the seniority-pairing strength, G . The latter is obtainedusing a method based on Ref. [69]. Specifically, to obtain an expression for G , we assume that the spacing betweenthe energy levels is constant (cid:15) k − ˜ λ = k − N pair ˜ ρ . (14)In this expression, ˜ λ is the smooth Fermi energy of the smoothed single-particle energy and ˜ g is the smooth leveldensity obtained with the Strutinsky method. The pairing gap ∆ in (11) is approximated by an effective pairing gap¯∆ = (cid:40) r mic B s N / for neutrons r mic B s Z / for protons . (15)We then substitute the sum of (11) by an integral L max (cid:88) k = L min f ( (cid:15) k − λ ) ≈ ˜ ρ (cid:90) y y f ( x ) d x (16)to obtain the following expression for GG = 2˜ ρ (cid:104) arsinh (cid:16) y ¯∆ (cid:17) − arsinh (cid:16) y ¯∆ (cid:17)(cid:105) − , (17)where y = L min − N pair − ˜ ρ (18) y = L max − N pair + ˜ ρ . (19)Note that these expressions differ from the one obtained in [41, 64], being shifted by 1 / (2˜ ρ ). The expression of theaverage pairing correlation plus quasi-particle energy ˜ E p . c . is then obtained by inserting these quantities into Eq. (110)of Ref. [64].We solve the Lipkin-Nogami equations (9)-(13) using a new method based on the analytical calculation of the fullLipkin-Nogami Jacobian and the action of its inverse on any vector coupled with a fifth-order numerical scheme. Weset multiple starting points to avoid local minima and find that our method greatly enhances the success rate of theresolution of the Lipkin-Nogami equations. Thus, the time required to solve the Lipkin-Nogami equations using ourmethod is faster than in older work. This means the calculation time for our procedure is now negligible comparedto the calculation and diagonalization time of the Hamiltonian, thus allowing for larger shape families or potentialenergy surfaces to be explored in the future. This new method is described in detail in Appendix B. B. Novel approach to strongly damped nuclear motion
Many different methods have been developed to calculate nascent fission fragment yields. Most of these methodsseek to simplify the time-evolution of the complex nuclear motion subject to various assumptions, see Ref. [50]for a recent overview. One of the most successful methods pursued over the past ten years is the assumption ofstrongly damped nuclear dynamics in which the Smoluchowski equations reduce to a Metropolis random walk (RW)[48, 70]. This method has been used in a large range of fission reactions using Markov Chain Monte Carlo (MCMC)sampling [49, 50, 71–73]. The primary drawback to this method is the long calculation time required for near-barrierfission. In what follows, we briefly review the construction of the finite-temperature PES and the MCMC samplingmethod. We then present our new Deterministic-Probabilistic Algorithm (DPA) for attaining scission configurationsunder the assumption of strongly damped motion.The finite-temperature potential energy U ( q ) for each points of the PES is obtained through the insertion of asuppression factor S [ E ∗ ( q )] in Eq. (1), U ( q ) = E mac ( q ) + S [ E ∗ ( q )]∆ E s+p ( q ) (20)∆ E s+p ( q ) = ∆ E shell ( q ) + ∆ E pair . ( q ) (21) S [ E ∗ ] = 1 + exp( − E /E )1 + exp(( E ∗ − E ) /E ) , (22)where E ∗ ( q ) = E ∗ − E ( q ) is the local excitation energy as in Ref. [49]. A discrete or continuous Random-Walk (RW)is then used on the finite-temperature PES U ( q ) to ascertain the scission configurations for the given incident energy, E ∗ .The MCMC RW procedure is illustrated in Fig. 1a. A single step can be decomposed as follows:1. start at an initial point q (state P );2. randomly choose one of its D neighbors q k ( k = 1 , . . . , D ) using a uniform distribution (transition from P toone of the I k );3. if(a) U ( q k ) ≤ U ( q ), then, P acc . (cid:1) k = 1 P rej . (cid:1) k = 0 (23)where P acc . (cid:1) k is the acceptance probability and P rej . (cid:1) k is the rejection probability. We accept q k as a the nextpoint (transition from state I k to another state P k );(b) U ( q k ) > U ( q ), then, we sample a random variable r ∼ U ([0 , r < P rej . (cid:1) k = 1 − P acc . (cid:1) k , where P acc . (cid:1) k = exp (cid:18) − U ( q k ) − U ( q ) T (cid:19) , (24)then go back to q (transition from I k to P ). Otherwise, accept q k as the next point (transition from I k to P k ).Results using this approach are always obtained through the generation of paths by explicitly repeating these stepsfrom an initial point up to reaching a scission shape. A scission shape is defined such that the radius of its neck, if P P D I D P I P I (a) Graph of the Markov-Chain transitions associated withthe Random-Walk method. P P D P P (b) Loop-less graph of the Markov Chain transitions of theMarkov-Chain associated with the Random-Walk method. present, is lower than a given parameter r neck . Many paths have to be calculated to reduce the statistical uncertaintyand obtain a reasonable estimation of the probability density function (PDF) for the scission configurations.We propose here to directly determine the evolution of the PDF according to the number of steps in a path.This transformation is analogous to the determination of the Fokker-Plank equations associated with the Langevinequations. The first step is to unravel the loops P (cid:1) I k (cid:1) P appearing in Fig. 1a. This can be visualized by thetransition graph of Fig. 1b.The probabilities P (cid:1) k to reach P k from P is the sum of the probabilities of all possible paths between P and P k .Such a path can be decomposed into two parts: the first one is the self-looping from P to P and the second oneis P (cid:1) I k (cid:1) P k and corresponds to the last two transitions of the path. Using the Markov property of memoryless transitions, the probability P (cid:1) k to reach P k from P through any path can be decomposed as the product of theprobability C to reach P starting from P which is independent of k and the probability to reach P k from P withoutgoing around any loop P (cid:1) k = C × P acc . (cid:1) k D . (25)and D is the enumeration of states directly connected with P . Using again the memoryless transitions property, itcan be shown that the probability to reach P after 2 N iterations is P (2 N )0 (cid:1) = (cid:34) D (cid:88) k =1 P rej . (cid:1) k D (cid:35) N , (26)which goes to zero as N goes to infinity. Therefore, we have D (cid:88) k =1 P (cid:1) k = C × D (cid:88) k =1 P acc . (cid:1) k D = 1 . (27)This relation enables the calculation of C without resorting to its definition as an infinite sum to obtain C = D (cid:80) Dk =1 P acc . (cid:1) k , (28)leading to P (cid:1) k = P acc . (cid:1) k (cid:80) Dk =1 P acc . (cid:1) k . (29)Note that these formulas can also be derived by directly using the definition of CC ≡ ∞ (cid:88) N =0 P (2 N )0 (cid:1) = 1 (cid:80) Dk =1 P acc . (cid:1) k D , (30)where Eq. (26), the Taylor series of the function x (cid:55)→ / (1 − x ) and the relation P rej . (cid:1) k = 1 − P acc . (cid:1) k have been used. Thecorresponding transition probabilities from a point q to q k are thus proportional to P acc . (cid:1) k . Our derivation shows thatin the standard RW method, there exists a parameter that correspond to an energy threshold, ∆ U thresh . = 0 MeV, onthe energy difference, U ( q k ) − U ( q ), below which a transition from I k to P k is certain to happen (with probability1).With the probability of an individual step well defined, we now seek to calculate the probability distribution aftera fixed number of RW steps. This quantity, p ( n +1) R ( q f ), is defined as the probability distribution associated with therandom variable R n +1 to reach a point, q f , after n + 1 RW steps and can be obtained recursively from the distribution p ( n ) R after n steps using p ( n +1) R ( q f ) = (cid:88) q P acc . q (cid:1) q f p ( n ) R ( q ) . (31)From this definition it is clear that the probability distribution after n steps is thus completely determined by thechoice of the initial probability distribution p (0) R .Rather than using Eq. (31) to directly compute our scission PDF, we generalize our procedure by allowing adistribution of ending configurations. This amounts to adding an absorption mechanism on top of the nuclear PES,analogous to the idea introduced in Ref. [30]. At each step n of our calculation we add an absorption field A ( q ) suchthat ¯ a ( n ) R ( q ) = ¯ a ( n − R + A ( q ) p ( n ) R ( q ) (32) a ( n ) R ( q ) = ¯ a ( n ) R ( q ) (cid:80) q (cid:48) ¯ a ( n ) R ( q (cid:48) ) (33)where q and q (cid:48) are two distinct shape configurations, ¯ a ( n ) R is the probability to reach a given shape at a given iteration, a ( n ) R ( q ) is the probability to reach a shape q after n iterations given that q is a scission shape. We let ¯ a (0) R ( q ) = 0 andthe values of the absorption field 0 ≤ A ( q ) ≤ q . The Eq. (31) is then replaced by p ( n +1) R ( q f ) = (cid:88) q P acc . q (cid:1) q f ¯ p ( n ) R ( q ) . (34)where the initial probability distribution at the next iteration ¯ p ( n ) R ( q ) is then taken to be¯ p ( n ) R ( q ) = (1 − A ( q )) p ( n ) R ( q ) . (35)In this first work, we set the absorption field, A ( q ), to 100% for scission configurations and zero otherwise. This allowsus to directly compare to previous work in the Section III. The absorption field can in principle be configured in anynumber of ways, including a distribution for various neck sizes, instead of the implicitly assumed sharp function ofthis work. We plan to study such possibilities in future work.Our novel method also affords the calculation of a convergence criteria. We can use this quantity as an estimateof the statistical uncertainty on calculated fission yields. To this end, we assume that the 1-distance (the distanceassociated with the 1-norm || . || ) between two distributions a ( n ) R and a ( n + j ) R , defined as∆ ( n )( j ) = (cid:88) q (cid:12)(cid:12)(cid:12) a ( n + j ) R ( q ) − a ( n ) R ( q ) (cid:12)(cid:12)(cid:12) , (36)is an inverse-quadratic function ∆ ( n )( j ) ≈ c ( j ) n + d ( j )] , (37)where j is an integer, c ( j ) and d ( j ) are two real parameters that are ultimately obtained using a fit procedure. Thevalidity of this assumption is discussed in III for j = 1000. The convergence error at a step n can be introduced asthe distance between the distribution at step n and the distribution at infinity ε ( n ) = ∆ ( n )( ∞ ) . (38)Using the subadditivity of the 1-distance (commonly referred as the triangle inequality), we can determine the con-vergence criteria, ε ( n ) ≤ ∞ (cid:88) k =0 ∆ ( n + kj )( j ) . (39)Inserting (37) into (39) gives ε ( n ) ≤ ∞ (cid:88) k =0 jc ( j )) k + ( c ( j ) n + d ( j ))] . (40)A closed form of the right-hand side of this expression can be obtained by resorting to the trigamma function definedas φ (1) ( z ) = d d z ln Γ( z ) , (41)that satisfies φ (1) ( z ) = ∞ (cid:88) k =0 k + z ) . (42)Employing these properties of the trigamma function, we finally obtain ε ( n ) ≤ ε ( n,j )eff . , (43)where ε ( n,j )eff . = 1[ jc ( j )] φ (1) (cid:18) c ( j ) n + d ( j ) jc ( j ) (cid:19) . (44)We take ε ( n,j )eff . as our convergence criteria as it is an upper limit on the error ε ( n ) after n iterations. C. Mass and charge yields of fission fragments
Most of the current mac-mic models used to estimate the probability distribution associated with the fission fragmentproperties before prompt emissions resort to obtaining the mass, Y ( A f ), or charge, Y ( Z f ), yields separately using therelevant macroscopic shape parameter, e.g., the procedure of Ref. [72]. While this method has been highly successful,see e.g. [48, 70], it does not provide a means to calculate the full fragment yield, Y ( Z f , A f ) or equivalently, Y ( Z f , N f ).In the past, the full mass and charge yields have obtained through the direct analysis of systematics on knownexperimental data [74], thus with low predictive power albeit high-quality data, or using the Wahl systematics [75]e.g., in Ref. [50], which introduces a free parameter σ Z that controls the dispersion in charge of the isobaric yields,and presumes the Unchanged Charge Distribution (UCD) assumption that relies on the ratio η ≡ Z f A f = ZA . Anothermethod, presented in Ref. [76], aims at obtaining the full mass and charge yields but relies on the addition of asixth macroscopic shape parameter in the PES. Without adding specific parameters, none of these past approachescan predict the fission fragments’ charge polarization, which is the experimentally observed deviation from the UCDassumption.In the following, we construct an approach to predict the full probability distribution of the fission fragmentmass and charge Y ( Z f , A f ) directly from the quantum mechanical wavefunctions. As it will be shown in Sec. III A,our projection technique is able to reproduce charge polarization and odd-even staggering in the fragments withoutadditional parameters.We only consider here the nascent fragments in binary fission. Therefore, Y ( Z f , A f ) can be decomposed accordingto the probability distribution, Y L ( Z f , A f ), (normalized to 1) associated with the mass and charge of only the leftfragment, Y ( Z f , A f ) = Y L ( Z f , A f ) + Y L ( Z − Z f , A − A f ) . (45)The law of total probability enables the decomposition of the probability distribution of the number of particles inthe left fragment Y L ( Z f , A f ) before prompt particle emission as follows, Y L ( Z f , A f ) = (cid:90) P (( Z, A ) L = ( Z f , A f ) | R = q ) a R ( q ) ( ∞ ) d q , (46)where the integral iterates over all the scission shapes parameterized by q , a ( ∞ ) R ( q ) is the limit for n → ∞ of thedistributions introduced in Eq. (33) and P (( Z, A ) L = ( Z f , A f ) | R = q ) is the probability associated with a left fragmentof mass A f and charge Z f when the fissioning system is in the shape q . Note that to use the law of total probability,we had to implicitly assume that there exists a set of random variables R , . . . , R ∞ , as in the previous section. Itis not always the case, such as in the Time-Dependent Generator Coordinate Method formalism [77]. In this case,coupling terms have to be added in the decomposition of Y L ( Z f , A f ).The probability P (( Z, A ) L = ( Z f , A f ) | R = q ) is extracted from the microscopic state calculated to estimate theshell+pairing correction for each coordinate q of the PES after projection on the good mass and charge of the totalfissioning system. Because all these states preserve the isospin, we can further decompose P (( Z, A ) L = ( Z f , A f ) | R = q ) = P ( N L = N f | R = q ) P ( Z L = Z f | R = q ) , (47)where N f is the number of neutrons in the left fragment. Both factors in the right-hand side of Eq. (47) are calculatedthrough the particle-number projection-based technique on the fragments mass and charge developed first in Ref. [51]in the context of time-dependent mean-field calculations for transfer reactions. This technique was first applied tofission in Ref. [78] and adapted to the case of static mean-field calculations in Ref. [52]. When applied to fission, thistechnique gives the probabilities associated with the number of nascent fragment neutrons ( X = N , X f = N f ) andprotons ( X = Z , X f = Z f ) using, P ( X L = X f | R = q ) = (cid:10) Φ( q ) (cid:12)(cid:12) ˆ P ( L ) X f ˆ P X (cid:12)(cid:12) Φ( q ) (cid:11)(cid:10) Φ( q ) (cid:12)(cid:12) ˆ P X (cid:12)(cid:12) Φ( q ) (cid:11) , (48)or equivalently P ( X L = X f | R = q ) = (cid:10) Φ( q ) (cid:12)(cid:12) ˆ P ( L ) X f ˆ P ( R ) X (cid:12)(cid:12) Φ( q ) (cid:11)(cid:10) Φ( q ) (cid:12)(cid:12) ˆ P X (cid:12)(cid:12) Φ( q ) (cid:11) . (49)A double projection is required, where ˆ P X = N,Z is the operator restoring the good number of particle in the totalsystem while ˆ P ( L ) X f is an operator projecting on X f particles in the left fragment. The definition of the latter relies onthe position of the neck along the symmetry axis. We define this quantity, in the standard way, as the position of theminimum of the local one-body density of the microscopic state in q between the two pre-fragments [34, 79, 80]. Theprojection-based method to calculate the fragment distribution is already known to describe the odd-even staggeringof the charge distribution of the fragments in case of time-dependent mean-field methods [78] and is able to givenon-zero probability for the existence of fragments with odd-number of particles [52]. At first sight, it could seem tobe a paradox in the case of static mean-field calculations in the case of even-even systems since:1. the state describing the fissioning system is time-even;2. the operators ˆ P X and ˆ P ( L ) X f are time-even;3. a state describing odd-number fragments cannot be time-even.0However, this paradox is only apparent and can be solved by noticing that the projection operators are both acting onthe full X -body wavefunction (for each isospin). Even in the case where X f is odd, the state ˆ P ( L ) X f ˆ P X | Φ( q ) (cid:105) containsan even number of particles: X f in the left fragment and X − X f in the right one. It is thus time-even. It can easilybe seen from a simple example of a time-even state having only two particles | Φ (cid:105) = ˆ a † i ˆ a † ¯ i | (cid:105) , (50)where ˆ a † i is the creation operator of a particle in state i , ¯ i is the time-reveral state of i . Both states can be decomposedin a similar way as in [51, 52] as ˆ a † i = α ( L ) i ˆ a ( L ) † i + α ( R ) i ˆ a ( R ) † i (51)ˆ a † ¯ i = α ( L ) (cid:63)i ˆ a ( L ) † ¯ i + α ( R ) (cid:63)i ˆ a ( R ) † ¯ i , (52)where ˆ a ( L ) † k and ˆ a ( R ) † k are respectively the left and right parts of ˆ a † k for k = i, ¯ i . The creation operators on the leftcommute with the ones on the right due to the complete separation of their spatial domain, and each of them commutewith their time-reversal. Therefore, by injecting (51) and (52) into (50) and developing the resulting expression, weobtain | Φ (cid:105) = α ( L ) i α ( L ) (cid:63)i ˆ a ( L ) † i ˆ a ( L ) † ¯ i | (cid:105) + (cid:0) α ( L ) i α ( R ) (cid:63)i ˆ a ( L ) † i ˆ a ( R ) † ¯ i + α ( R ) i α ( L ) (cid:63)i ˆ a ( R ) † i ˆ a ( L ) † ¯ i (cid:1) | (cid:105) + α ( R ) i α ( R ) (cid:63)i ˆ a ( R ) † i ˆ a ( R ) † ¯ i | (cid:105) . (53)The three terms are orthogonal to each other and are all time-even. The first one corresponds to a state with twoparticles in the left side and zero in the right one, the second one corresponds to one particle on each side and the lastone corresponds to two particles in the right side and zero in the left one. Therefore, even though | Φ (cid:105) is time-evenand ˆ P ( L )1 is time-even, we have a non-zero probability to have odd-number fragments when the state i spreads onboth the left and right domains P ( X L = X f ) = (cid:10) Φ (cid:12)(cid:12) ˆ P ( L ) X f (cid:12)(cid:12) Φ (cid:11) (cid:104) Φ | Φ (cid:105) = 2 | α ( L ) i α ( R ) i | . (54)However, as shown and explained in Ref. [52] in the case of static time-even Bogoliubov states with an even-numberof particles, the probability associated with odd-number fragments collapse to zero as soon as the fragments areseparated enough and don’t interact anymore. This is a direct consequence of the finite-range character of the nuclearinteraction and the minimization of the energy: if two subsystems S and S of a system S do not interact with eachother, the energy of the total system is the sum of the energies of both subsystems and thus, the state that minimizethe energy of S is the product of the states minimizing each subsystems. III. APPLICATION TO 233,235-U(N,F)
We illustrate our model improvements in what follows by showcasing the well-known neutron-induced fission of twoisotopes of Uranium,
U(n,f). We also prove that our implementation can reproduce the results of the discreterandom-walk method used in past work.In Refs. [41, 63, 64, 71–73, 81], the authors implemented shell-plus-pairing corrections through the resolution of theSchr¨odinger equation using an Axial Harmonic Oscillator (AHO) basis with only N sh = 12 shells. The limitation ofthis previous approach artificially introduces spurious contributions of the continuum [66]. The use of the improvedStrutinsky method [65] allows us to remove these contributions and thus use instead N sh = 20 shells without anyenergy truncation. In addition, we optimize the oscillator scaling factor b and deformation q for each point of thePES using a variational principle such that the mean-field-plus-pairing energy of the microscopic state is minimized.The order of the Strutinsky method is p = 8, and the corresponding range is γ = C sr C cur A / B s , (55)where the relative surface energy B s ( q ) is the ratio of the nucleus surface at shape q with the surface of the samenucleus at spherical shape. We have used the same parameters as in Ref. [64], listed in the following table.1 Parameter Value Unit C sr C cur
41 MeV
TABLE I: Microscopic parameters associated with the Strutinsky correction.The pairing correction is obtained using the Lipkin-Nogami method. We have solved the Lipkin-Nogami (LN) foreach point of the PES using a pairing window of ± r mic = 3.2 MeV. The full LN equations are often numerically solvedby splitting them into two or more subsets of equations, solved separately at each iteration. It adds overhead in theresolution time and can also lead to spurious divergences. Instead, we have developed a new method to solve theseequations. Our method is based on the iterative fifth-order method presented in Ref. [82] that only requires first-orderderivatives (i.e., the Jacobian matrix). The main steps of our method are presented in Appendix B.The fragment probabilities at each q are calculated using Eq. (49). The double-projection on the numerator iscalculated using the Pfaffian technique presented in Refs. [83, 84]. The determination of the integrals over the gaugeangles are determined through a Fomenko discretization method [85]. For each scission shape, we find the optimalnumber of integration point N Fom . by assuming that the width of the fragment distribution associated with bothisospin is N Fom . . We check if this assumption is correct by checking that P ( X L = X f | R = q ) < . × − (56)for all X f = (cid:98) X mean − N Fom . / (cid:99) + ∆ X and X f = (cid:100) X mean + N Fom . / (cid:101) − ∆ X for ∆ X = 0 , , N Fom . that satisfiesit, starting at N Fom . = 30.The PES at a given excitation energy is obtained through the finite temperature method of Ref. [49] where thedamping of the shell-plus-pairing correction as in Eq. (20) invokes the damping parameter S [ E ∗ ] defined in Eq. (22).The two parameters we have taken to define the damping coefficient are E = 20 MeV and E = 15 MeV. Theexcitation energy dependency of the temperature is taken at the Thomas-Fermi approximation to be E ∗ ( q ) = aT ( q ) , (57)with the nuclear level density a = A/
8. To obtain an implementation equivalent to the state-of-the-art random-walk,we have set ∆ U thresh . = 0 MeV. Recall that this threshold does not directly appear in the standard formalism ofthe random-walk and corresponds to the energy difference U ( q k ) − U ( q ) between an initial point q and one of itsneighbor q k below which a transition from the states I k to P k is certain in the Markov Chain Fig. 1a. The potentialenergy surface is calculated on a regular grid following the work of Ref. [63]. Two points q and q are neighbors ifall the integer coordinates on the lattice differ by at most one unit. Such a definition in five dimensions leads to amaximum of 3 − q loc ;3. we determine the minimum energy E sad . required to reach scission configurations from q loc ;4. we calculate the set C of all the configurations accessible from q loc with an energy lower than E sad . , and definethe effective ground-state as the node q g . s . ∈ C associated with the lowest energy E g . s . .By using this procedure, we also obtain the saddle energy E sad , as well as an effective barrier height E B = E sad − E g . s . .A site of the lattice is a scission configuration if its corresponding sharp macroscopic density has a neck radius r neck < r sciss . . We have calculated the fragment probability distribution for r sciss = 1 . , . , .
75 fm (58)and three different excitation energies E ∗ such that x = E ∗ − E B = 0 . , . , . Z0510152025 Y ( Z ) ( % ) r neck = 1.75fm r neck = 2.25fm r neck = 2.75fmQuade 1988Wehring 198330 32 34 36 38 40 42Z05101520 Y ( Z ) ( % ) x = 0.1MeV x = 2.0MeV x = 4.0MeVQuade 1988Wehring 1983 (a) Fragment charge yields obtained with our approach forthe reaction U( n th . ,f) at x = 0 . r neck = 2 .
25 fm at differentexcitation energies (bottom panel). Our results are comparedwith experimental data with thermal incident neutronenergies ( E neut . = 0 . Z0102030 Y ( Z ) ( % ) r neck = 1.75fm r neck = 2.25fm r neck = 2.75fmLang 198032 34 36 38 40 42Z05101520 Y ( Z ) ( % ) x = 0.1MeV x = 2.0MeV x = 4.0MeVLang 1980 (b) Fragment charge yields obtained with our approach forthe reaction U( n th . ,f) at x = 0 . r neck = 2 .
25 fm at differentexcitation energies (bottom panel). Our results are comparedwith experimental data with thermal incident neutronenergies ( E neut . = 0 . A. Results
We have used our approach to calculate the fission fragment charge and mass probability distribution before promptemission for the reactions
U(n,f) and
U(n,f). In Figs. 2a and 2b we highlight the charge yields of these reactionsrespectively. The upper panels of Fig. 2a and Fig. 2b show an important variation of the yields according to thedefinition of the scission line for both reactions. This behavior is due to the minimization of the energy to obtainour microscopic states [52]. The lower panels of Fig. 2a and Fig. 2b exhibit very little variation of our charge yieldsaccording to the excitation energy in the range 0 . − . r neck = 2 .
25 fm, our chargeyields present a very good quantitative agreement with experimental data for nearly all proton numbers.As previously mentioned, our approach enables the determination of the mass and charge yields Y ( Z, A ). Ourcharge yields for different fixed fragment masses are presented in Fig. 3 for the reaction
U(n,f) and
U(n,f).Our independent yields (solid orange lines) are obtained by assuming that the prompt neutron emission multiplicitydistribution, P n , depends on the mass of the fragments only. To estimate it, we first fit the parameter p of theprobability mass function of the binomial random variable B ( N = 5 , p ) on experimental data from Ref. [88] assuggested in [89]. We then shift the distribution for each fragment mass such that the expected value is equal to¯ ν ( A ) from Refs. [90, 91] to obtain P n ( A ). Finally, the pre-neutron yields are corrected according to P n ( A ) to obtainour estimation of independent yields. This method is easy to implement and fast, albeit simple; providing a meansfor comparison with experimental independent yield data. The presented theoretical yields correspond to x = 0 . U( n th . ,f) and U( n th . ,f), where the neck condition is r neck = 2 .
25 fm. The main feature in these resultsis the reproduction of the charge polarization of the fragments. To estimate the charge polarization of the fragments,we have first calculated the average number of protons for each mass. Then, we have linearly fit it for the light- andheavy-fragments on all masses associated with a probability greater than 1%. The result of the fit is the upper bluelines for the light fragment and the lower ones for the heavy fragment. In each case, the expression of the linear curveis Z Lf = η L A + s L (60) Z Hf = η H A + s H , (61)3 A = 86(4.2%)Quade 1988Wehring 1983 This work This work (ind.) A = 87(5.2%)050 A = 88(6.1%) A = 89(6.8%)050 A = 90(7.5%) A = 91(8.1%)050 A = 92(8.5%) A = 93(8.4%)050 Y ( Z , A ) ( % ) A = 94(8.0%) A = 95(7.2%)050 A = 96(6.0%) A = 97(4.6%)050 A = 98(3.2%) A = 99(2.1%)30 32 34 36 38 40 42Z050 A = 100(1.2%) 30 32 34 36 38 40 42Z A = 101(0.6%) 050 A = 87(3.7%)Lang 1980 This work This work (ind.) A = 88(4.6%)050 A = 89(5.5%) A = 90(6.4%)050 A = 91(7.2%) A = 92(7.9%)050 A = 93(8.4%) A = 94(8.6%)050 Y ( Z , A ) ( % ) A = 95(8.3%) A = 96(7.7%)050 A = 97(6.7%) A = 98(5.4%)050 A = 99(4.0%) A = 100(2.6%)32 34 36 38 40 42Z050 A = 101(1.6%) 32 34 36 38 40 42Z A = 102(0.9%) FIG. 3: Isobaric fragment charge yields (dashed blue lines) and the independent charge yields (solid orange lines),after prompt neutron emission, obtained with our approach for the reactions U( n th . ,f) (left panel) and U( n th . ,f) (right panel) with a neck condition z neck = 2 .
25 fm and an excitation energy above the barrier of x = 0 . A to 100% and the fragment probability associated with each ofthe masses is listed in parentheses. We compare our results with experimental independent yields (dashed blackcurve) from Refs. [86, 87, 92].where Z L / Hf is the charge of the light/heavy fragment. Since we only consider binary fission, η L = η H ≡ η . Ifthe UCD assumption is satisfied (recall η = ZA ) then η ≈ . U(n,f) and η ≈ . U(n,f). However, the values of s L and s H are equal only if the UCD assumption is correct and the mass and chargedistributions of the light- and heavy-fragments are aligned. Therefore, we introduce the quantity ∆ Z = s L − s H ,which measures the charge alignment and vanishes if the UCD is verified. The ∆ Z values, along with the one of η , s L and s H for each excitation energy x are reported in Table II for the reaction U( n th . ,f) and Table II for thereaction U( n th . ,f). For these two reactions, our approach indicates the light and heavy fragments distributions aremisaligned by ∆ Z > . B. Validity of the convergence criteria
To analyze the convergence properties of our new algorithm, we calculate the probability distributions of the scissionconfigurations a n S ( q ) for all steps, n ≤ n max = 150 , U(n,f). From this calculations, we can
U(n,f)
U(n,f) x (MeV) η s L s H ∆ Z η s L s H ∆ Z TABLE II: Characteristics of the charge alignment between the light- and heavy- fragment distributions for thereactions
U(n,f) and
U(n,f) represented by the blue curves of Fig. 4. If the UCD is verified, η ≈ . U(n,f), η ≈ . U(n,f), and ∆ Z = 0.4 Z x = 0.1MeV30405060 Z x = 2.0MeV80 90 100 110 120 130 140 150 160A30405060 Z x = 4.0MeV 10 F r a g m e n t p r o b a b ili t y ( % ) Z x = 0.1MeV30405060 Z x = 2.0MeV80 90 100 110 120 130 140 150 160A30405060 Z x = 4.0MeV 10 F r a g m e n t p r o b a b ili t y ( % ) FIG. 4: Fragment mass and charge yields, before prompt neutron emission, obtained with our approach for thereactions U( n th . ,f) (left panel) and U( n th . ,f) (right panel) with a neck condition z neck = 2 .
25 fm. Each panelshows the yield with different excitation energies above the barrier. (a) Evolution of the scission probability according to thenumber of DPS iterations for the reaction
U(n,f). (b) Inverse-square-root of the evolution of the distancebetween successive probability distributions of the scissionconfigurations for the reaction
U(n,f). To keep a readablegraph, we have displayed only one point over four. A linearfit of the data is presented for comparison purposes. extract the probability to reach scission after n iterations P sciss .n = (cid:88) q ¯ a ( n ) R ( q ) . (62)The evolution of this quantity with increasing n is shown in Fig. 5a. As expected, the lower excitation energy (redcurve) is significantly below the run with higher excitation energy (green curve). For the same range of iterations,the evolution of the inverse-square-root of ∆ ( n )(1000) is presented in Fig. 5b. Recall that ∆ ( n ) j is the convergence errorgiven by Eq. (36).5The inverse-quadratic regime for ∆ ( n )(1000) is reached around n ≈ ,
000 iterations where the ratio of paths that arereaching scission is 3 . × − % for x = 0 . . × − % for x = 2 . ε ( n, . , we have estimated it through a fit on the data applicable inthe range of iterations n , . . . , n . We compare our convergence criteria against data through the calculation of therelative error of the convergence criteria using j = 1000 D n ,n = (cid:20) ε ( n ,j )eff . (cid:12)(cid:12)(cid:12) n n − ε ( n max +1 ,j )eff . (cid:12)(cid:12)(cid:12) n n (cid:21) − n max − n j (cid:88) k =0 ∆ ( n + jk )( j ) ε ( n ,j )eff . (cid:12)(cid:12)(cid:12) n n − ε ( n max +1 ,j )eff . (cid:12)(cid:12)(cid:12) n n , (63)where ε ( n,j )eff . (cid:12)(cid:12)(cid:12) n n corresponds to the estimation of our convergence criteria using a fit between n and n . The evolutionof D n ,n according to the range used for the fit is presented Fig. 6 for which n > ,
000 and n − n ≥ , D n ,n associated with our convergence criteria between the iterations n and n is below 2% after only10,000 iterations in both cases. D n ,n is greater than 0 almost everywhere, which means that our criteria slightlyoverestimates the convergence error.FIG. 6: Relative error D n ,n , defined by Eq. (63), between our convergence criteria and the exact convergenceerror, for each interval n , . . . , n of fitted data, for the reaction U(n,f) at an excitation energy above the barrierof x = 0 . x = 2 . C. Comparison with the Metropolis implementation
We compare our new calculations to past work in order to show that we can reproduce these efforts within thecontext of our more general methodology. We use the implementation of a discrete random walk (DRW) as inRef. [50] as the baseline FRLDM mass yield calculations. We perform these calculations for the reaction
U(n,f) atan excitation energy of x ≈ .
25 fm. We accumulate 100,000 scission configurationsfor this comparison. With these parameters, we have nearly identical inputs as our new results shown in Sec. III A.While the older mass yields rely only on the mass asymmetry coordinate, α g , it is not sufficient to compare only thisvariable at scission as there could be changes in the distribution of other coordinates. Figure 7a shows the distributionof the scission configurations in the full collective coordinate lattice space, ( i, j, k, l, n ), between the standard randomwalk method and the method presented in this work. Despite the differences of the two approaches, exceedingly good6 (a) Probability distributions to obtain a scissionconfiguration associated with each integer indices ( i, j, k, l, n )as defined in [63] estimated with state-of-the-art random walkwith the code DRW (blue) and with this work (orange). (b) Absolute error between the probability distributionspresented in Fig. 7a. agreement is found between the two algorithms. The absolute error between the two approaches is shown in Fig. 7b.We find that the statistical nature of the DRW algorithm leads to a maximum of ∼
2% error in the distribution of thescission neck radius, while the statistical error in α g is the lowest of all the coordinates, on the order of .1%. Thesetwo figures show that we are successfully able to reproduce past work with our new technique and that older worksindeed have quantitatively very good estimates of the mass yields within the context of FRLDM so long as a largenumber of scission events are calculated. IV. CONCLUSION
With this work, we have improved the quality and predictive power of the mac-mic method in several areas. First,we have enhanced the quality of the nuclear PES by removing spurious continuum effects in our five-dimensionalfinite-range model. Further, our new resolution procedure of the Lipkin-Nogami equations enables the descriptionof pairing effects with very high accuracy. Second, our new Deterministic-Probabilistic Algorithm (DPA) completelyremoves statistical uncertainties when computing the fission fragment distribution of a particular nucleus. DPA enablesthe starting point of our calculation to be located at the ground state (easily identifiable for all nuclei) without therequirement of including a biased potential that artificially tilts the PES. We have defined a high-accuracy convergencecriterion associated with our algorithm that affords ability to monitor the error associated with the obtained results.Last, but not least, we have generalized the particle number projection technique introduced for independent quasi-particle states in Ref. [51, 78] and calculated scission configurations with it. This projection technique allows for thecalculation of the coupled fragment charge and mass yield, Y ( Z, A ). We refer to these improvements colloquially asthe ‘Enhanced Finite-Range Liquid-Drop Model’ or eFRLDM for short.Our first eFRLDM results are presented for the pre-neutron fission fragments probability distributions of thereactions
U(n,f) at different excitation energies. Our method can reproduce the odd-even staggering in thecharge yields as well as the charge polarization of the fragments without any additional free parameters in the model.We find that a charge misalignment exists between the light and heavy fragments for these two reactions on the orderof ∆ Z > . ACKNOWLEDGMENTS
The authors would like to thank Nicolas Schunck for his advice, his careful reading of the manuscript and his helpfulcomments, as well as Patrick Talou, Jørgen Randrup, Toshihiko Kawano, Ionel Stetcu, Arnie Sierk, and Peter M¨ollerfor valuable discussions over the years.M.V. and M.R.M. were supported by the US Department of Energy through the Los Alamos National Laboratory.Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear SecurityAdministration of U.S. Department of Energy (Contract No. 89233218CNA000001). M.V. and M.R.M. were partlysupported by the Fission In R-process Elements (FIRE) topical collaboration in nuclear theory, funded by the U.S.Department of Energy and through the Los Alamos Laboratory Directed Research & Development ExploratoryResearch project entitled “Low-energy fission dynamics”. This work was partly performed under the auspices of theU.S Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.
Appendix A: Finite-Range Liquid Drop Model (FRLDM)
We group here the main formulas and parameters that define our model that we use to calculate the nuclearpotential energy surface. The macroscopic liquid-drop energy for even-even nuclei is E mac ( q ) = M H Z + M n N − a v (1 − κ v I ) A + a s (1 − κ s I ) B ( q ) A / + a A B W ( q )+ c Z A / B ( q ) − c Z / A / + f ( k f r p ) Z A − c a ( N − Z )+ W | I | B W ( q ) − a el Z . . (A1)In this expression, A , Z , N are respectively the number of nucleons, protons and neutrons and I is the relative neutronexcess, I = N − ZA . (A2)We note the pairing term is zero for even-even nuclei, and thus does not appear in the formula.We shift globally the energy the PES such that the energy of the spherical point is zero. In this case, only theshape-dependent terms contribute to the PES. There are four such terms in our approach:81. the surface energy, calculated assuming a finite-range Yukawa-plus-exponential nuclear interaction with nofolding [93],2. the Coulomb term, defined with a Yukawa folding of the sharp macroscopic density [94],3. the A energy from Ref. [41],4. the Wigner term from Ref. [41].By setting σ a = | r − r (cid:48) | a (A3)and letting V represent the sharp-macroscopic density, the shape-dependent energies are the relative surface energyassociated with a Yukawa-plus-exponential finite-range two-body interaction, B ( q ) = − A − / π r a (cid:90) (cid:90) V [ σ a − e − σ a σ a d r d r (cid:48) , (A4)the relative Coulomb energy of a folded-Yukawa macroscopic density, B ( q ) = 15 A − / π a den r (cid:90) (cid:90) V d r d r (cid:48) σ a den (cid:34) − (cid:16) σ a den (cid:17) e − σ a den (cid:35) , (A5)and the shape-dependency of the A and Wigner terms is defined as B W ( q ) = (cid:16) − S S (cid:17) a d + 1 if there is a neck1 otherwise. (A6)In the last expression, S S q at the neck location. The definition of B W is slightly different than in Refs. [41, 50]where the condition is only relative to the MQS parameter, σ . However, σ < z -axis. Defining l , l and l as the respective centers of the left, middle and right bodies of the MQS shape, in the case where l ≤ l or l ≤ l , the shape cannot exhibit a neck whatever the sign of σ . In this situation, B W defined as in older work is notcontinuous at σ = 0. Our new definition prevents this situation.The model parameters we have used for the calculation of the PES are often referred to as FRLDM2002 whichcorrespond to the model parameters introduced in Ref. [64] with additional corrections of Ref. [81]. We present hereonly the parameters having an influence on the shape-dependent terms of the PES. The values of the fundamentalconstants we have used are Table IV references the parameters associated with the macroscopic part of the energy,
Parameter Value Unit e m amu m nuc m n m p TABLE III: Truncation of the fundamental constants.taken from [64]. Table V reports the macroscopic parameters taken from [81].The parameters associated with the potential V ( r ; q ) defined in (2) are Appendix B: Resolution of the Lipkin-Nogami equations
When using a seniority-pairing interaction, the Lipkin-Nogami equations associated with a valence space of N v energy levels are the set of 2 N v + 3 nonlinear equations (9)-(13) with the same number of unknowns v k , (cid:15) k , ∆, λ Parameter Value Unit r a a den W TABLE IV: Part of the macroscopic parameters used in our approach to describe the fission process from [64].
Parameter Value Unit a d a s κ s a TABLE V: Part of the macroscopic parameters used in our approach to describe the fission process from [81].and λ . Some of these equations are associated with high-derivatives. To reduce the amplitude and the number ofnon-zeros derivatives, we substitute u k , v k and (cid:15) k by the variable x k and θ k according to u k = cos( θ i ) (B1) v k = sin( θ i ) (B2) x k = (cid:15) k − λ. (B3)The Lipkin-Nogami equations can then be rewritten and reorganized as F ( p ) = 0 , (B4)where, setting ¯ k = k − L min , F N v ( p ) = L min − N pair + L max (cid:88) k = L min sin ( θ k ) (B5) F N v +2 ( p ) = (cid:34) L max (cid:88) k = L min √ x k + ∆ (cid:35) − G (B6) F k +1 ( p ) = 12 (cid:20) − x k √ x k + ∆ (cid:21) − sin( θ k ) (B7) F k ( p ) = (4 λ − G ) sin( θ k ) + e k − x k − λ (B8) F N v +1 ( p ) = A ( θ ) λ − G B ( θ ) , (B9)where the notations A ( θ ) = (cid:34) L max (cid:88) k = L min cos( θ k ) sin( θ k ) (cid:35) − L max (cid:88) k = L min cos( θ k ) sin( θ k ) (B10) B ( θ ) = (cid:34) L max (cid:88) k = L min cos( θ k ) sin( θ k ) (cid:35) × (cid:34) L max (cid:88) k = L min cos( θ k ) sin( θ k ) (cid:35) − L max (cid:88) k = L min cos( θ k ) sin( θ k ) , (B11)and p = (cid:16) θ L min , x L min , . . . , θ L max , x L max , λ, λ , ∆ (cid:17) . (B12)The analysis of the dependencies of each equations leads to a maximum of 10 N v + 2 non-vanishing elements in theJacobian matrix J F ( p ). Also, J F p is block-arrowhead , which means that J F = (cid:18) A BC D (cid:19) , (B13)0 Parameter Value Unit V s V a A den B den a pot k p l p k n l n a a J L Q K TABLE VI: Microscopic parameters associated with the potential (2).where A is a block-diagonal matrix. In our case, the blocks of A are 2-dimensional matrices A ¯ k = (cid:32) ∂F k ∂θ k ∂F k ∂x k ∂F k +1 ∂θ k ∂F k +1 ∂x k (cid:33) . (B14)The block-column matrix B and row-column matrix C are respectively associated with the following 2 × × B ¯ k = (cid:32) ∂F k ∂λ ∂F k ∂λ
00 0 ∂F k +1 ∂ ∆ (cid:33) (B15) C ¯ k = ∂F N v ∂θ k ∂F N v+1 ∂θ k ∂F N v+2 ∂x k , (B16)In the following, we propose a method to solve the Lipkin-Nogami equations and more generally any system ofequations associated with a block-arrowhead Jacobian matrix at each point p based on generalizations of the iterativeNewton method. In our case, we use the cubic and the fifth-order iterative methods developed respectively byHomeier [95] and by Sharma & Gupta [82]. The idea of these methods is to improve the convergence properties ofthe Newton scheme by evaluating the Jacobian at different p . For example, in the Sharpa & Gupta scheme, one stepfrom iteration i to i + 1 is x ( i ) (cid:0) p ( i ) − J F ( p ( i ) ) − F ( p ( i ) ) (B17) y ( i ) (cid:0) p ( i ) − J F ( x ( i ) ) − F ( p ( i ) ) (B18) p ( k +1) (cid:0) y ( i ) − (cid:104) aJ F ( x ( i ) ) − + bJ F ( y ( i ) ) − (cid:105) F ( y ( i ) ) , (B19)where the fifth-order convergence is obtained when a = 2 and b = −
1. The Newton scheme is recovered by doing onlythe first step and p ( i +1) (cid:0) x ( i ) , while the two first steps are present in the Homeier scheme and p ( i +1) (cid:0) y ( i ) .The blockwise inversion theorem gives the inverse of the Jacobian matrix as J F ( p ) − = (cid:18) A − + IM − IT − − M T − (cid:19) , (B20)where I = A − B (B21) T = D − CI (B22) M = T − CA − . (B23)1Note that it is assumed that A and T are invertible. when it is not the case, we slightly perturb the diagonal elementsof the non-invertible matrix. This method requires inversion of N v A ¯ k and one 3-dimensionalmatrix T . However, the procedure gives a dense matrix. Instead, we directly calculate the four vectors a ( i ) = J F ( p ( i ) ) − F ( p ( i ) ) (B24) b ( i ) = J F ( x ( i ) ) − F ( p ( i ) ) (B25) c ( i ) = J F ( x ( i ) ) − F ( y ( i ) ) (B26) d ( i ) = J F ( y ( i ) ) − F ( y ( i ) ) , (B27)such that x ( i ) (cid:0) p ( i ) − a ( i ) (B28) y ( i ) (cid:0) p ( i ) − b ( i ) (B29) p ( k +1) (cid:0) y ( i ) − (cid:104) a c ( i ) + b d ( i ) (cid:105) . (B30)In the following, we note e = a ( i ) , b ( i ) , c ( i ) or d ( i ) and f = F ( p ( i ) ) or F ( p ( i ) ) according to the equation (B24)-(B27)considered. By injecting (B20) in (B24)-(B27) we obtain e = A − f + IM f − IT − f (B31) e = − M f + T − f , (B32)where e = (cid:18) e e (cid:19) f = (cid:18) f f (cid:19) , (B33)and e and f are vectors of dimension 2 N v and e and f are vectors of dimension 3. The expression of e appearsin the expression of e . Therefore, once e is obtained by using (B32), e can be obtained through the expression e = A − f − I e . (B34)Lastly, we note our convergence criteria is set to ε LN ( p ) < − , where, using the aforementioned functions F k , ε LN ( p ) = max k | F k ( p ) | . (B35) [1] N. Tsoulfanidis and R. G. Cochran, Nuclear Technology , 263 (1991).[2] “Neutron chain fission reactors,” in Nuclear Reactor Physics (John Wiley & Sons, Ltd, 2018) Chap. 2, pp. 33–42.[3] R. C. Pardo, G. Savard, S. Baker, C. Davids, E. F. Moore, R. Vondrasek, and G. Zinkann, Nuclear Instruments andMethods in Physics Research B , 965 (2007).[4] J. Carlson, M. P. Carpenter, R. Casten, C. Elster, P. Fallon, A. Gade, C. Gross, G. Hagen, A. C. Hayes, D. W. Higinbotham,C. R. Howell, C. J. Horowitz, K. L. Jones, F. G. Kondev, S. Lapi, A. Macchiavelli, E. A. McCutchen, J. Natowitz,W. Nazarewicz, T. Papenbrock, S. Reddy, M. A. Riley, M. J. Savage, G. Savard, B. M. Sherrill, L. G. Sobotka, M. A.Stoyer, M. B. Tsang], K. Vetter, I. Wiedenhoever, A. H. Wuosmaa, and S. Yennello, Progress in Particle and NuclearPhysics , 68 (2017).[5] B. M. Sherrill, in European Physical Journal Web of Conferences , European Physical Journal Web of Conferences, Vol.178 (2018) p. 01001.[6] R. Surman and M. Mumpower, EPJ Web of Conferences (2018), 10.1051/epjconf/201817804002.[7] E. P. Abel, M. Avilov, V. Ayres, E. Birnbaum, G. Bollen, G. Bonito, T. Bredeweg, H. Clause, A. Couture, J. DeVore,M. Dietrich, P. Ellison, J. Engle, R. Ferrieri, J. Fitzsimmons, M. Friedman, D. Georgobiani, S. Graves, J. Greene, S. Lapi,C. S. Loveless, T. Mastren, C. Martinez-Gomez, S. McGuinness, W. Mittig, D. Morrissey, G. Peaslee, F. Pellemoine, J. D.Robertson, N. Scielzo, M. Scott, G. Severin, D. Shaughnessy, J. Shusterman, J. Singh, M. Stoyer, L. Sutherlin, A. Visser,and J. Wilkinson, Journal of Physics G: Nuclear and Particle Physics , 100501 (2019).[8] G. Mention, M. Fechner, T. Lasserre, T. A. Mueller, D. Lhuillier, M. Cribier, and A. Letourneau, Phys. Rev. D , 073006(2011).[9] A. C. Hayes, J. L. Friar, G. T. Garvey, G. Jungman, and G. Jonkmans, Phys. Rev. Lett. , 202501 (2014). [10] M. Itkis, J. Aysto, S. Beghini, A. Bogachev, L. Corradi, O. Dorvaux, A. Gadea, G. Giardina, F. Hanappe, I. Itkis, M. Jan-del, J. Kliman, S. Khlebnikov, G. Kniajeva, N. Kondratiev, E. Kozulin, L. Krupa, A. Latina, T. Materna, G. Montagnoli,Y. Oganessian, I. Pokrovsky, E. Prokhorova, N. Rowley, V. Rubchenya, A. Rusanov, R. Sagaidak, F. Scarlassara, A. Ste-fanini, L. Stuttge, S. Szilner, M. Trotta, W. Trzaska, D. Vakhtin, A. Vinodkumar, V. Voskressenski, and V. Zagrebaev,Nuclear Physics A , 136 (2004).[11] Y. Oganessian and V. Utyonkov, Nuclear Physics A , 62 (2015), special Issue on Superheavy Elements.[12] K. Godbey and A. S. Umar, Frontiers in Physics , 40 (2020).[13] C. Simenel, A. Wakhle, B. Avez, D. J. Hinde, R. du Rietz, M. Dasgupta, M. Evers, C. J. Lin, and D. H. Luong, in European Physical Journal Web of Conferences , European Physical Journal Web of Conferences, Vol. 38 (2012) p. 09001.[14] G. Mart´ınez-Pinedo, D. Mocelj, N. T. Zinner, A. Keli´c, K. Langanke, I. Panov, B. Pfeiffer, T. Rauscher, K.-H. Schmidt,and F.-K. Thielemann, Progress in Particle and Nuclear Physics , 199 (2007).[15] S. Goriely, J.-L. Sida, J.-F. Lemaˆıtre, S. Panebianco, N. Dubray, S. Hilaire, A. Bauswein, and H.-T. Janka, Phys. Rev.Lett. , 242502 (2013).[16] M. Eichler, A. Arcones, A. Kelic, O. Korobkin, K. Langanke, T. Marketin, G. Martinez-Pinedo, I. Panov, T. Rauscher,S. Rosswog, et al. , The Astrophysical Journal , 30 (2015).[17] M. Mumpower, T. Kawano, T. Sprouse, N. Vassh, E. Holmbeck, R. Surman, and P. M¨oller, The Astrophysical Journal , 14 (2018).[18] N. Vassh, R. Vogt, R. Surman, J. Randrup, T. M. Sprouse, M. R. Mumpower, P. Jaffke, D. Shaw, E. M. Holmbeck, Y. Zhu, et al. , Journal of Physics G: Nuclear and Particle Physics , 065202 (2019).[19] N. Vassh, M. R. Mumpower, G. C. McLaughlin, T. M. Sprouse, and R. Surman, Astrophys. J. , 28 (2020).[20] P. Bonche, S. Koonin, and J. W. Negele, Phys. Rev. C , 1226 (1976).[21] Y. Engel, D. Brink, K. Goeke, S. Krieger, and D. Vautherin, Nuclear Physics A , 215 (1975).[22] R. Y. Cusson and H. W. Meldner, Phys. Rev. Lett. , 694 (1979).[23] S. Levit, J. W. Negele, and Z. Paltiel, Phys. Rev. C , 1979 (1980).[24] A. Bulgac, P. Magierski, K. J. Roche, and I. Stetcu, Phys. Rev. Lett. , 122504 (2016).[25] C. Simenel and A. S. Umar, Phys. Rev. C , 031601(R) (2014).[26] P. Goddard, P. Stevenson, and A. Rios, Phys. Rev. C , 054610 (2015).[27] C. Simenel and A. Umar, Progress in Particle and Nuclear Physics , 19 (2018).[28] Y. Tanimura, D. Lacroix, and S. Ayik, Phys. Rev. Lett. , 152501 (2017).[29] N. Schunck and L. M. Robledo, Reports on Progress in Physics , 116301 (2016).[30] J. Berger, M. Girod, and D. Gogny, Computer Physics Communications , 365 (1991).[31] H. Goutte, J. F. Berger, P. Casoli, and D. Gogny, Phys. Rev. C , 024316 (2005).[32] D. Regnier, N. Dubray, N. Schunck, and M. Verri`ere, Phys. Rev. C , 054611 (2016).[33] D. Regnier, N. Dubray, and N. Schunck, Phys. Rev. C , 024611 (2019).[34] W. Younes and D. Gogny, Phys. Rev. C , 054313 (2009).[35] W. Younes, D. M. Gogny, and J.-F. Berger, A microscopic theory of fission dynamics based on the generator coordinatemethod , Vol. 950 (Springer, 2019).[36] J. Zhao, T. Nikˇsi´c, D. Vretenar, and S.-G. Zhou, Phys. Rev. C , 014618 (2019).[37] K.-H. Schmidt, B. Jurado, C. Amouroux, and C. Schmitt, Nuclear Data Sheets , 107 (2016), special Issue on NuclearReaction Data.[38] J.-F. Lemaˆıtre, S. Panebianco, J.-L. Sida, S. Hilaire, and S. Heinrich, Phys. Rev. C , 034617 (2015).[39] J.-F. Lemaˆıtre, S. Goriely, S. Hilaire, and J.-L. Sida, Phys. Rev. C , 034612 (2019).[40] P. M¨oller, W. D. Myers, H. Sagawa, and S. Yoshida, Phys. Rev. Lett. , 052501 (2012).[41] P. Mller, A. Sierk, T. Ichikawa, and H. Sagawa, Atomic Data and Nuclear Data Tables , 1 (2016).[42] V. M. Strutinsky, Nuclear Physics A , 420 (1967).[43] V. M. Strutinsky, Nuclear Physics A , 1 (1968).[44] M. BRACK, J. DAMGAARD, A. S. JENSEN, H. C. PAULI, V. M. STRUTINSKY, and C. Y. WONG, Rev. Mod. Phys. , 320 (1972).[45] T. Wada, N. Carjan, and Y. Abe, Nuclear Physics A , 283 (1992).[46] P. N. Nadtochy, A. Keli´c, and K.-H. Schmidt, Phys. Rev. C , 064614 (2007).[47] A. J. Sierk, Phys. Rev. C , 034603 (2017).[48] J. Randrup and P. M¨oller, Phys. Rev. Lett. , 132503 (2011).[49] J. Randrup and P. M¨oller, Phys. Rev. C , 064606 (2013).[50] M. R. Mumpower, P. Jaffke, M. Verriere, and J. Randrup, Phys. Rev. C , 054607 (2020).[51] C. Simenel, Phys. Rev. Lett. , 192701 (2010).[52] M. Verriere, N. Schunck, and T. Kawano, Phys. Rev. C , 024612 (2019).[53] R. W. Hasse and W. D. Myers, Geometrical relationships of macroscopic nuclear physics (Springer Science & BusinessMedia, 1988).[54] A. V. Ramayya, J. H. Hamilton, J. K. Hwang, L. K. Peker, J. Kormicki, B. R. S. Babu, T. N. Ginter, A. Sandulescu,A. Florescu, F. Carstoiu, W. Greiner, G. M. Ter-Akopian, Y. T. Oganessian, A. V. Daniel, W. C. Ma, P. G. Varmette, J. O.Rasmussen, S. J. Asztalos, S. Y. Chu, K. E. Gregorich, A. O. Macchiavelli, R. W. Macleod, J. D. Cole, R. Aryaeinejad,K. Butler-Moore, M. W. Drigert, M. A. Stoyer, L. A. Bernstein, R. W. Lougheed, K. J. Moody, S. G. Prussin, S. J. Zhu,H. C. Griffin, and R. Donangelo, Phys. Rev. C , 2370 (1998).[55] A. V. Ramayya, J. K. Hwang, J. H. Hamilton, A. Sandulescu, A. Florescu, G. M. Ter-Akopian, A. V. Daniel, Y. T. Oganessian, G. S. Popeko, W. Greiner, J. D. Cole, and GANDSCollaboration (GANDS95 Collaboration), Phys. Rev.Lett. , 947 (1998).[56] S. Vermote, C. Wagemans, O. Serot, J. Heyse, J. Van Gils, T. Soldner, and P. Geltenbort, Nuclear Physics A , 1(2008).[57] S. Vermote, C. Wagemans, O. Serot, J. Heyse, J. Van Gils, T. Soldner, P. Geltenbort, I. AlMahamid, G. Tian, and L. Rao,Nuclear Physics A , 176 (2010).[58] F. Go et al. , Nuclear Physics A , 213 (2004).[59] J. Nix, Further studies in the liquid-drop theory of nuclear fission , Tech. Rep. LBNL Report , 241 (1969).[61] J. R. Nix, Annual Review of Nuclear Science , 65 (1972).[62] M. Bolsterli, E. O. Fiset, J. R. Nix, and J. L. Norton, Phys. Rev. C , 1050 (1972).[63] P. M¨oller, A. J. Sierk, T. Ichikawa, A. Iwamoto, R. Bengtsson, H. Uhrenholt, and S. ˚Aberg, Phys. Rev. C , 064304(2009).[64] P. Moeller, J. R. Nix, W. D. Myers, and W. J. Swiatecki, Atomic Data and Nuclear Data Tables (1995),10.1006/adnd.1995.1002.[65] A. Kruppa, Physics Letters B , 237 (1998).[66] N. Tajima, Y. R. Shimizu, and S. Takahara, Phys. Rev. C , 034316 (2010).[67] H. J. Lipkin, Annals of Physics , 272 (1960).[68] Y. Nogami, Phys. Rev. , B313 (1964).[69] D. G. Madland and J. Nix, Nuclear Physics A , 1 (1988).[70] J. Randrup, P. M¨oller, and A. J. Sierk, Phys. Rev. C , 034613 (2011).[71] P. M¨oller, J. Randrup, and A. J. Sierk, Phys. Rev. C , 024306 (2012).[72] P. M¨oller, J. Randrup, A. Iwamoto, and T. Ichikawa, Phys. Rev. C , 014601 (2014).[73] P. M¨oller and J. Randrup, Phys. Rev. C , 044316 (2015).[74] M. Caama˜no, F. Rejmund, and K.-H. Schmidt, Journal of Physics G: Nuclear and Particle Physics , 035101 (2011).[75] A. C. Wahl, Los Alamos Report (2002), LA-13928.[76] P. M¨oller and T. Ichikawa, The European Physical Journal A , 1 (2015).[77] M. Verriere and D. Regnier, Frontiers in Physics , 233 (2020).[78] G. Scamps, C. Simenel, and D. Lacroix, Phys. Rev. C , 011602(R) (2015).[79] D. Regnier, N. Dubray, M. Verriere, and N. Schunck, Computer Physics Communications , 180 (2018).[80] R. N. Perez, N. Schunck, R.-D. Lasseri, C. Zhang, and J. Sarich, Computer Physics Communications , 363 (2017).[81] P. M¨oller, A. J. Sierk, and A. Iwamoto, Phys. Rev. Lett. , 072501 (2004).[82] J. R. Sharma and P. Gupta, Computers & Mathematics with Applications , 591 (2014).[83] G. F. Bertsch and L. M. Robledo, Phys. Rev. Lett. , 042505 (2012).[84] L. M. Robledo, Phys. Rev. C , 021302(R) (2009).[85] V. Fomenko, Journal of Physics A: General Physics , 8 (1970).[86] U. Quade, K. Rudolph, S. Skorka, P. Armbruster, H.-G. Clerc, W. Lang, M. Mutterer, C. Schmitt, J. Theobald, F. Gn-nenwein, J. Pannicke, H. Schrader, G. Siegert, and D. Engelhardt, Nuclear Physics A , 1 (1988).[87] B. W. Wehring, S. Lee, R. Strittmatter, and G. Swift, Trans. Am. Nucl. Soc.; (United States) (1980).[88] R. Gwin, R. Spencer, and R. Ingle, Nuclear Science and Engineering , 381 (1984).[89] B. C. Diven, H. C. Martin, R. F. Taschek, and J. Terrell, Phys. Rev. , 1012 (1956).[90] K. Nishio, M. Nakashima, I. Kimura, and Y. Nakagome, Journal of nuclear science and technology , 631 (1998).[91] K. Nishio, Y. Nakagome, H. Yamamoto, and I. Kimura, Nuclear Physics A , 540 (1998).[92] W. Lang, H.-G. Clerc, H. Wohlfarth, H. Schrader, and K.-H. Schmidt, Nuclear Physics A , 34 (1980).[93] H. J. Krappe, J. R. Nix, and A. J. Sierk, Phys. Rev. C , 992 (1979).[94] K. T. R. Davies and J. R. Nix, Phys. Rev. C , 1977 (1976).[95] H. Homeier, Journal of Computational and Applied Mathematics169