Non-perturbative collective inertias for fission: a comparative study

The non-perturbative method to compute Adiabatic Time Dependent Hartree Fock Bogoliubov (ATDHFB) collective inertias is extended to the Generator Coordinate Method (GCM) including the case of density dependent forces. The two inertias schemes are computed along the fission path of the $^{234}$U and compared with the perturbative results. We find that the non-perturbative schemes predict very similar collective inertias with a much richer structure than the one predicted by perturbative calculations. Moreover, the non-perturbative inertias show an extraordinary similitude with the exact GCM inertias computed numerically from the energy overlap. These results indicate that the non-perturbative inertias provide the right structure as a function of the collective variable and only a phenomenological factor is required to mock up the exact GCM inertia, bringing new soundness to the microscopic description of fission.


Introduction
Despite its discovery dates back almost 80 years, fission still remains a major challenge for nuclear theory [1]. The lack of a feasible full quantum formalism describing the evolution of the nucleus from the ground state to scission enforces the adoption of different approximations, which in turn provide a theoretical framework for the estimation of fission properties in nuclei. For instance, the starting point in any traditional energy density functional calculation is the original assumption that fission can be described using a reduced set of collective variables [2,3]. Within this approximation the fission probability is obtained as the probability of the nucleus to tunnel under the fission barrier, which is driven by the potential energy surface (PES) and the collective inertias felt by the nucleus in its way to scission [4,5]. Both quantities, together with the collective ground-state energy, enter in the collective action integral allowing for the calculation of the spontaneous fission lifetime by means of the semiclassical Wentzel-Kramers-Brillouin (WKB) approach.
A sound calculation of the PES, collective ground-state energy and collective inertias is thus essential for a proper estimation of fission lifetimes [1,6]. If the formalism that shall be used in the calculation of the first two quantities is well established, the same cannot be claimed for the collective inertias. Nowadays two theoretical frameworks allow for a derivation of a collective Schrödinger equation and its associated inertia: the Adiabatic Time Dependent Email addresses: giuliani@nscl.msu.edu (Samuel A. Giuliani), luis.robledo@uam.es (Luis M. Robledo) Hartree Fock Bogoliubov (ATDHFB) formalism and the Generator Coordinate Method (GCM) with the Gaussian overlap approximation (GOA) [1]. In both approaches the collective inertias can be written in terms of the collective momentum operators, which in turn can be related to the linear response matrix (LRM) and its inverse when acting on Hartree Fock Bogoliubov (HFB) wave functions. Given the whopping number of two-quasiparticle elementary excitations in realistic applications to fission, the dimensionality of the LRM is very high and therefore its inverse is difficult to evaluate. To avoid this bottleneck the assumption of diagonal dominance of the LRM is often used, leading to the traditional perturbative cranking formulas for the collective inertias involving denominators composed of two quasiparticle energies. A better approximation to the exact expression of the collective inertias was introduced in [7] (see also [8]), where the collective momentum operator is computed in terms of the derivatives of the density and pairing tensor with respect of the collective variables. This non-perturbative cranking calculation of the collective inertias, implemented in the AT-DHFB approach, showed that the numerical treatment of the derivatives gives rise to a less adiabatic behaviour of the collective inertias in comparison to the perturbative calculation.
It follows then that there are two different sources of uncertainty in the calculation of the collective inertias: one related to the choice of the theoretical framework (ATD-HFB vs GCM-GOA) and the other related to the approximations involved in the numerical evaluation of the inertias (exact vs non-perturbative vs perturbative). The purpose of this paper is then twofold: (i) to introduce for the first time the non-perturbative scheme in the GCM-GOA framework and (ii) compute the exact GCM-GOA collective inertias and use this result to study the suitability of both the ATDHFB and the GCM-GOA non-perturbative schemes. Using the actinide 234 U as a test case, we will show that, as the level of approximation improves, the results obtained in the ATDHFB and the GCM schemes naturally converge towards the same solution of the collective inertias, bringing new solidity to the theoretical description of fission. The present results represent a step forward in the microscopic description of fission providing the method with the credibility required to answer questions like the very existence of nuclei beyond oganesson [9].

Methodology
This section is devoted to the derivation of the different expressions used for the calculation of the collective inertias. The key element is the momentum operatorP q associated to the collective variable q which is derived in the quasiparticle representation in section 2.1. This quantity is then used to obtain the non-perturbative expression of the GCM-GOA mass in section 2.2 while the extension to density dependent forces is presented in section 2.3. In section 2.4 we discuss how to compute the exact GCM-GOA mass using the numerical derivatives of the exact Hamiltonian and norm kernels. In section 2.5 we will briefly review the derivation of the ATDHFB non-perturbative formula. Section 2.6 is devoted to obtaining the explicit expression of the perturbative masses, both in the GCM-GOA and ATDHFB framework. Finally, the connection between collective inertias and moments of inertia is presented in section 2.7.

Momentum operator
Given a collective variable q like, for instance, the quadrupole moment, its associated collective momentP q can be defined through the relation (we use = 1 in the following) Both |Φ(q + δq) and |Φ(q) are HFB wave functions satisfying the HFB equation with the corresponding constraints Φ(q+δq)|Q 20 |Φ(q+δq) = q+δq and Φ(q)|Q 20 |Φ(q) = q, respectively. To evaluate |Φ(q + δq) in terms of |Φ(q) we use linear response theory in the quasiparticle representation. We notice that |Φ(q + δq) and |Φ(q) are related by a Thouless transformation that, for infinitesimal δq, can be written as On the other hand, |Φ(q + δq) satisfies the HFB equation with constraints where the quasiparticles creation operators above are defined at deformation q + δq Expanding in powers of δq we obtain at zero order which is an identity because |Φ(q) is the constrained HFB solution at deformation q. At first order in δq the following identity has to be satisfied as well as its complex conjugated. In the above expression ∆Ô =Ô − Ô , (Q j ) 20 * µν = β ν β µQj is the 20 part of the operatorQ j , andĤ =Ĥ − j λ jQj . Introducing the matrices with the properties A µ ν µν = A * µνµ ν and B µ ν µν = B µ ν µν , Eq. (7) becomes µ <ν µ <ν To simplify the notation it is convenient to introduce indexes ρ and σ corresponding to the pair of indexes µ and ν with the restriction µ < ν. The ordering of the correspondence is irrelevant in what follows. With the new indexes, Z µν becomes the vector Z ρ and the four index quantities A µ ν µν become the matrix elements of a hermitian matrix A ρ ρ . The same applies to the B µ ν µν that become the matrix elements of a symmetric matrix B ρ ρ . In terms of the new indexes the previous equation becomes Introducing the linear response matrix (LRM) L which is closely related to the matrix appearing in the Random Phase Approximation (RPA), it is easy to express Z in terms of the partial derivatives of the chemical potentials The partial derivatives of the chemical potentials are determined by plugging the above result in the definition of the constraints As a consequence of this requirement we get where the quantity M (−1) ij is given by Collecting together all the partial results we finally obtain The evaluation of the momentum matrix elements requires the inversion of L which is in general a tremendous task, given the typical number of two quasiparticle excitations involved in a realistic calculation. An alternative (and useful) expression for the momentum operator (or Z) can be obtained by evaluating the derivatives of the densities (both normal and abnormal) with respect to the constraints (see Appendix).

The Generator Coordinate Method inertia
The GCM does not directly provide an expression of the collective inertia. It is only after introducing some local approximation that the Hill-Wheeler equation can be reduced to a collective Schrödinger equation and yield the associated inertia [10]. Traditionally, the GOA is the approximation of choice to make this connection.
Assuming that the width of the Gaussian does not depend on q the GCM-GOA mass is given by [1,6,10] 1 where as required to preserve particle number on the average also for the GCM wave functions [11,12]. If the constant width is not assumed [13,11,12] the above expression remains valid, but one has to replace the partial derivatives by covariant ones that include in their definition the affine connection or Christoffel symbols of differential geometry. We will use in the following the constant width formula to preserve the traditional connection with the momentum operator defined above. Assuming time reversal invariant states |φ(q) such that φ(q)| ∂ ∂q |φ(q ) |q =q = 0 and computing second derivatives of the HFB states as [14] are omitted) we finally obtain that leads to the compact expression Please note that the A and B matrices above are not exactly the same as those of Eqs. (8) which are defined in terms ofĤ instead ofĤ eff . The differences, associated with the collective constrains, are zero for the ground state and very small elsewhere as we have checked in our example below. In the following we will assume them to be the same. Using the definition of Z, L and introducing the matrix η = 1 0 0 −1 we finally obtain that is written in a way that can be easily generalized to the multidimensional case. TheM (−1) is given bȳ The width γ can be obtained in a similar manner: where M (−2) is defined in analogy with Eq. (15) but replacing L −1 by L −2 . In the non-perturbative cranking approach we use in Eq. (19) the Z obtained from the partial derivatives of the density matrix and pairing tensor (see Eq. (A.8)). Additionally, we use the cranking approximation for L where B = 0 and A is replaced by its diagonal Inserting this approximation into the general equation we arrive to that is the expression used in this paper. We do not use BCS like approximations like the one discussed in Ref [7].

Density dependent forces
For density dependent forces like Gogny or Skyrme the above formulation has to be slightly modified. In Eq. (4) the HamiltonianĤ(q + δq) has to be replaced byĤ(q + δq) + ∂Γ(q + δq) where the one body rearrangement term is given by ∂Γ(q + δq) = ij ∂Γ(q + δq) ij c † i c j with matrix elements When expressing those quantities in terms of the corresponding ones at deformation q, derivatives with respect to q of bothĤ(q) and ∂Γ(q) have to be considered. The expressions obtained are rather involved but straightforward to derive and they will not be given here. In addition, those derivatives only enter the LRM and therefore they are not required in the non-perturbative case except for the definition of the one-quasiparticle energies that must be computed with the Hamiltonian including the rearrangement termĤ(q) + ∂Γ(q).

The GCM-GOA inertia
To compute the GCM-GOA inertia without using the cranking approximation we use Eq. (17) evaluating the derivatives numerically. The required Hamiltonian overlap in Eq. (18) is evaluated using the expressions of the generalized Wick theorem [15]. For the phenomenological density dependent part of the Gogny force we use the mixed density prescription as discussed in Refs [16,17]. First order finite difference formulas are used for the sec- with a value of h conveniently chosen according to the collective variable used (see below). The width γ is computed numerically in the same way from the norm overlap.

The Adiabatic Time Dependent HFB inertia
The ATDHFB inertia [18] can be evaluated using the same framework as above, but imposing additional constraints on the momentum operators [19]. We are not going to provide the details here, but following the same steps as above in the quasiparticle picture one gets Introducing the matrix which is very similar in structure toM (−1) lm of Eq. (20), we obtain A method for the exact evaluation of the ATDHFB inertia has been formulated in [20] and applied to very simple cases. A more recent attempt based on a direct evaluation of the LRM and its numerical inversion seems to be rather impractical due to the enormous computational cost [21].
Recently, it has been suggested [22] that the finite amplitude method [23] could be useful for this task, but so far, it has only been applied to the evaluation of the Thouless-Valatin moment of inertia.

The perturbative approximation
In the perturbative cranking approximation the L matrix is approximated by its diagonal also in the definition of the momentum operator Eq. (16). The approximate L commutes with η and the quantity defined in Eq. (20) becomes the M (−1) defined in Eq. (15) which in turn becomes M PE (−1) with the momenta or order n defined as

Connection with rotational band moments of inertia and inequalities
We would like to mention the similarity between the non-perturbative inertias and the Inglis-Belyavev and approximate Yoccoz moments of inertia [10]. In this case, the "momentum operator" is dictated by symmetry considerations (it is the J x operator) and therefore the approximate expressions used in the literature to compute moments of inertia fall into the non-perturbative cranking category discussed here [10]. Concerning the exact moments of inertia, the Thouless-Valatin moment of inertia is obtained in a similar framework to the ATDHFB case, whereas the Yoccoz moment of inertia corresponds to the GCM-GOA inertia. A comparison of the two moments of inertia [14,24] reveals that the Thouless-Valatin is typically a factor 1.4 larger than the Yoccoz one both when computed exactly and when computed in the nonperturbative cranking spirit.
In [25] it was shown, using the Schwarz inequality, that at the minima of the potential energy surface the exact inertias satisfy M ATDHFB ≥ M GCM . The same inequality was proved true for the whole fission path in the case of the perturbative inertias. In the case of the non-perturbative inertias and following the same arguments as in [25] it is straightforward to show that the same inequality also applies.

Results
In this section we compare the numerical results obtained for the quadrupole collective inertia in the typical case of the fission of the actinide 234 U.
In Fig 1 the results for 234 U are shown as a function of the mass quadrupole moment Q 20 expressed in barns. In panel a) the potential energy surface, given by the HFB energy, is depicted: the characteristic normal deformed minimum at Q 20 = 12 b (β 2 = 0.25) is obtained, followed by a fission isomer (or super-deformed state) at Q 20 = 40 b (β 2 = 0.71) with an excitation energy of 4.2 MeV. A very broad and high (9.5 MeV) second fission barrier is found next. In panel b) the particle-particle energy E pp defined as 1 2 tr(∆ * κ) is given separately for protons (full) and neutrons (dashed). A rather intricate behaviour is observed with the quadrupole moment, due to different level crossings that increase or decrease the level density around the Two main conclusions regarding the collective inertias can be extracted from panels c) and d). The first one is that the ATDHFB np and GCM-GOA np cranking inertias have exactly the same peak structure. As expected (see the discussion in Sec. 2.7), and in both the NP and PE cases, the ATDHFB mass is larger than the GCM-GOA one, as shown in panel e). What is remarkable is that in the two numerical schemes (NP and PE) the ATDHFB and GCM-GOA inertias differ by a factor around 1.5 and rather constant over the whole deformation range up to the region corresponding to two separate fragments (Q 20 > 126 b). After this point on the quadrupole moment has a geometric origin (being proportional to the square of the distance between the two fragments) and the ATDHFB np /GCM-GOA np ratio remains very close to one. We also point out that the 1.5 value of the ratio is consistent with other studies concerning the values of the Thouless-Valatin and Yoccoz moments of inertia (see section 2.7). The second finding is that both the ATDHFB np and GCM-GOA np inertias have more pronounced peaks compared to the perturbative calculations. Looking at the variations of the particle-particle energy E pp it is possible to relate these peaks with a larger sensitivity of the non-perturbative inertias to the presence of level crossings. This "lack of adiabaticity" of the non-perturbative inertias is consistent with the ATDHFB results of Baran et al. [7]. Finally, the fluctuations of the GCM-GOA np /GCM-GOA pe ratio depicted in panel e) (dashed line) suggest that the recipe of multiplying the perturbative inertia by a constant factor in order to simulate the non-perturbative masses is not a reasonable assumption [26]. The agreement between the non-perturbative inertias diminishes the uncertainties arising from the ambiguity in the choice of the theoretical scheme (ATDHFB vs GCM-GOA), but still the suitability of this numerical approximation has to be proved. In order to address this point we computed the exact GCM-GOA collective inertias (GCM-GOA) and compared the results with the perturbative and non-perturbative calculations. The lower panel of figure 2 shows the different inertias computed in this work and the upper panel represents the ratio of the perturbative and non-perturbative calculations to the exact GCM-GOA collective inertias. Surprisingly the GCM-GOA inertias have the same peak structure and evolution with quadrupole deformation of the non-perturbative calculations, with a ATDHFB np /GCM-GOA ratio very close to one (or equivalently GCM-GOA np /GCM-GOA∼ 0.6) and virtually independent of Q 20 . This result brings consistency to the nuclear fission theory, indicating that the GCM-GOA inertia can be obtained either by using the ATDHFB np scheme or the GCM-GOA np mass multiplied by a constant factor around 1.5. On the other hand, the upper panel of figure 2 shows that the ratio of perturbative cranking to GCM-GOA masses depends on the quadrupole deformation, with discrepancies in some cases as large as a factor of 5. This comparison confirms the results found in the non-perturbative study and indicating the inadequacy of multiplying the perturbative inertias by a phenomenological factor to grasp the structure of the exact GCM-GOA collective inertia [26].

Conclusion
In summary, this work provides a solution to the uncertainties arising from the ambiguity in the choice of both the theoretical framework and the numerical approximations involved in the calculation of collective inertias. Taking 234 U as a benchmark, the non-perturbative cranking and exact collective inertias were calculated for the first time using the Generator Coordinate Method (GCM) in the Gaussian overlap approximation (GOA) and compared with the Adiabatic Time Dependent Hartree Fock Bogolyubov (ATDHFB) non-perturbative and perturbative cranking inertias.
The ATDHFB np , GCM-GOA np and GCM-GOA inertias present the same peak structure along the whole fission path, being the GCM-GOA np calculations smaller by a roughly constant factor around 1.5. These inertias show a much richer structure compared to the perturbative calculations indicating a stronger sensitivity to level crossings. These results are not only important for fission but also for approximate models used to describe collective dynamics within the Bohr Hamiltonian or the Collective Schrödinger equation. The use of the non-perturbative inertias along with a phenomenological stretching factor can be a good substitute for the more elaborated beyond mean field calculations with the GCM-GOA. It would be highly desirable to extend this comparison to the exact ATDHFB inertia, but the complications arising from the calculation of the inverse of the linear response matrix L prevent this comparison for the moment. Work aimed to compute exact ATDHFB inertias is under way and will be reported elsewhere.

Acknowledgement
We are grateful to G. Martinez-Pinedo for his continuous encouragement in all the stages of this work as well as for many suggestions and improvements. SAG acknowledges support from the U.S. Department of Energy under Award Number DOE-DE-NA0002847 (NNSA, the Stewardship Science Academic Alliances program), the Deutsche Forschungsgemeinschaft through contract SFB 1245, and the BMBF-Verbundforschungsprojekt number 05P15RDFN1. The work of LMR has been supported in part by Spanish grant Nos FIS2015-63770 MINECO and FPA2015-65929 MINECO.
Appendix A. Derivatives of the density matrix The relationship between the momentum operator matrix elements and the derivative of the densities with respect to the constraints is established. Let us consider the density ρ ij (q) and the pairing tensor κ ij (q) corresponding to a set of values of the constraints q. By shifting one of the constraints q k by an infinitesimal δq we have to consider ρ ij (q k + δq) = φ(q k + δq)|c † j c i |φ(q k + δq) φ(q k + δq)|φ(q k + δq) . (A.1) With Eqs. (2,3) and the contractions φ(q)|c † j β † µ |φ(q) = V jµ and φ(q)|c j β † µ |φ(q) = U jµ we easily arrive to and In order to solve for Z k in the above expressions it is far more convenient to work with the unitary block matrix W of Bogoliubov amplitudes W = U V * V U * (A.4) and the associated generalized density matrix Using them we can write the expressions for the derivatives in a compact way as that can straightforwardly be solved for Z leading to (A.8)