Microscopic theory of nuclear fission: a review

In the MM approach, the energy (1) is a function of all the parameters needed to describe the nuclear shape. In other words, the nuclear surface must be parametrized explicitly. Numerous parametrizations have been introduced over the years; see [55] for a comprehensive review. In this section, we wish to recall some of the most common parametrizations, especially those that have been introduced to describe extremely deformed nuclear shapes near scission.

The simplest and most common parametrization is based on the multipole expansion of the nuclear radius [56],

$$\begin{eqnarray}&&R(\theta,\varphi )={{R}_{0}}c(\alpha )\left[1+\underset{\lambda \geqslant 2}{\sum}\,\underset{\mu =-\lambda}{\overset{+\lambda}{\sum}}\,{{\alpha}_{\lambda \mu}}{{Y}_{\lambda \mu}}(\theta,\varphi )\right],\end{eqnarray} \tag{ 2 }$$

where ${{R}_{0}}={{r}_{0}}{{A}^{1/3}}$ is a parametrization of the nuclear radius for a spherical nucleus of mass A (${{r}_{0}}\approx 1.2$ fm); $c(\alpha )$ is a factor accounting for the conservation of nuclear volume with deformation; the ${{\alpha}_{\lambda \mu}}$ are the deformation parameters; and ${{Y}_{\lambda \mu}}(\theta,\varphi )$ are the usual spherical harmonics. This parametrization of the nuclear surface is ideal for small to moderate deformations. However, a well-known problem of this formulation discussed in [57, 58] is that it does not provide a unique representation of the nuclear surface. In addition, the number of active deformation parameters needed to characterize very elongated shapes can become large. These difficulties also manifest themselves in the DFT framework since, as we will discuss in section 2.3.1, the collective variables most commonly used are the mass multipole moments of the nucleus, which are closely related to the spherical harmonics of (2). For illustration, we show in figure 2 a cross-section of the shapes obtained with the expansion (2) for the nucleus ²⁴⁰Pu with $\lambda =2,3,4$ and $\mu =0$.

Zoom In Zoom Out Reset image size

Because of the limitations of expansion (2), alternative parametrizations of the nuclear surface have been advocated over the years. One of the most popular is the Funny Hills shapes of [9]. Assuming axial symmetry along the z-axis of the intrinsic reference frame, a point on the surface is characterized by the usual cylindrical coordinates $(\rho,z,\varphi )$. In this case, the distance ρ from the nuclear surface to the axis of symmetry is given by

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\rho}^{2}}=R_{0}^{2}{{c}^{2}}\left(1-{{\xi}^{2}}\right)\left[A+\alpha \xi +B{{\xi}^{2}}\right], & B\geqslant 0, \end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\rho}^{2}}=R_{0}^{2}{{c}^{2}}\left(1-{{\xi}^{2}}\right)[A+\alpha \xi ]{{\text{e}}^{B{{c}^{3}}{{\xi}^{2}}}}, & B<0, \end{array}\end{eqnarray} \tag{ 4 }$$

where A, B, and α are the parameters characterizing the shape; c is related to the elongation of the system, and defines the dimensionless variable $\xi =z/\left(c{{R}_{0}}\right)$; in practice, A is obtained from the volume conservation condition, and the parameter B is substituted by another, $h=\frac{1}{2}B-\frac{1}{4}(c-1)$, which is related to the thickness of the neck. Therefore, the Funny Hills parametrization is most commonly characterized by the set $(c,h,\alpha )$. The figure 3 illustrates the shapes obtained for $(h,\alpha )=(0,0)$ (top panel) and $(h,\alpha )=(0.2,0.3)$ (bottom panel), each case with c varying between 1.0 and 2.2 by step of 0.3.

Zoom In Zoom Out Reset image size

An alternative to the Funny Hills parametrization is the one proposed by the Los Alamos group in [7]. It also applies to axially-symmetric shapes. In cylindrical coordinates, the nuclear surface is described by three quadratic surfaces of revolution

$$\begin{eqnarray}{{\rho}^{2}}=\left\{\begin{array}{*{35}{l}} a_{1}^{2}-\frac{a_{1}^{2}}{c_{1}^{2}}{{\left(z-l_{1}^{2}\right)}^{2}}, & {{l}_{1}}-{{c}_{1}}\leqslant z\leqslant {{z}_{1}}, \\ a_{2}^{2}-\frac{a_{2}^{2}}{c_{2}^{2}}{{\left(z-l_{2}^{2}\right)}^{2}}, & {{z}_{2}}\leqslant z\leqslant {{l}_{2}}+{{c}_{2}}, \\ a_{3}^{2}-\frac{a_{3}^{2}}{c_{3}^{2}}{{\left(z-l_{3}^{2}\right)}^{2}}, & {{z}_{1}}\leqslant z\leqslant {{z}_{2}}. \end{array}\right. \end{eqnarray} \tag{ 5 }$$

The nine parameters ${{\left({{a}_{i}},{{l}_{i}},{{c}_{i}}\right)}_{i=1,2,3}}$ are not all independent: the condition of volume conservation eliminates one parameter; the condition of continuity and differentiability of the surface at the two interfaces z₁ and z₂ eliminate two more parameters; fixing the centre of mass of the shape can also eliminate an additional parameter, so that only five independent parameters are ultimately needed. They are denoted $\left({{Q}_{2}},{{\alpha}_{g}},{{\varepsilon}_{f1}},{{\varepsilon}_{f2}},d\right)$ and correspond respectively to the elongation of the whole nucleus, the mass asymmetry of the two nascent fragments, the axial quadrupole deformation of the left and right fragment, and the thickness of the neck; see [49] for details and derivations.

The last category of nuclear shape parametrization originally introduced in [59] is based on Cassini ovals. The general expression for Cassini ovals is

$$\begin{eqnarray}&&{{\left({{a}^{2}}+{{z}^{2}}+{{\rho}^{2}}\right)}^{2}}=4{{a}^{2}}{{z}^{2}}+{{b}^{4}},\end{eqnarray} \tag{ 6 }$$

with a and b the only parameters. If a  =  0, then the shape reduces to a sphere of radius b; if a  >  b (6) represents two separate fragments. When using this parametrization to represent the nuclear shape, the volume conservation condition allows eliminating one of the two parameters. Therefore, Cassini ovals represent a single-parameter parametrization of the nucleus that can reproduce shapes from the spherical point to two distinct fragments separated by an arbitrary distance. To account for mass asymmetry, Cassini ovals can be distorted by substituting ${{a}^{2}}\to {{a}_{2}}-{{\epsilon}^{2}}$ in the left-hand side of (6). The figure 4 illustrates some of the typical shapes obtained in the symmetric case of (6).

Zoom In Zoom Out Reset image size

Once a parametrization of the nuclear surface has been established, the macroscopic part of the energy ${{E}_{\text{mac}}}\left(\boldsymbol{q}\right)$ can be computed. In the most advanced parametrization of the macroscopic energy, the surface, surface-symmetry, Coulomb and Wigner terms depend on the deformation via geometrical form factors that can be computed in a straightforward manner when the nuclear shape is known. The other terms are deformation-independent and adjusted to global nuclear properties [50–52].

In addition to the macroscopic energy, the MM approach relies on various quantum corrections. Here, the adjective 'quantum' recalls the fact that these corrections can only be computed by solving the Schrödinger equation of quantum mechanics.

The shell correction is the most important correction. It was introduced in [60] in order to account for the single particle shell structure in the calculation of the total energy. The starting point of the method is a phenomenological nuclear potential $U\left(\boldsymbol{r}\right)$ including a central part such as e.g. the Nilsson, Woods–Saxon or folded-Yukawa potential, a spin–orbit potential and the Coulomb potential for protons. The geometry of these potentials should be consistent with that of the nuclear surface: in the case of the Woods–Saxon or folded-Yukawa potential, this implies that the various terms depend on the distance of the current point to the nuclear surface, parametrized as discussed in the previous section. Given the potential and its deformation, one solves the one-body Schrödinger equation in order to obtain the (deformed) single particle energies e_n and wave functions ${{\varphi}_{n}}\left(\boldsymbol{r}\right)$. The last step is to approximate the discrete sequence of states ${{\left\{{{e}_{n}}\right\}}_{n}}$ by a continuous distribution $\tilde{g}(e)$ following the Strutinsky averaging procedure. The shell correction is then defined as

$$\begin{eqnarray}&&{{E}_{\text{shell}}}=\underset{n}{\sum}\,{{e}_{n}}-{\int}^{}\tilde{g}(e)e\text{d}e.\end{eqnarray} \tag{ 7 }$$

In the last expression, the shell correction should be computed for both protons and neutrons; the summation extends either over the Z and N lowest orbitals (ground-state) or over a subset of orbitals (excited states). In addition to providing a more realistic estimate of binding energies, the introduction of a one-body potential has given birth to a very powerful phenomenology based on single particle orbitals as basic degrees of freedom of nuclei; see [11] for a reference textbook on this topic.

Pairing correlations are incorporated into the MM approach in the form of an average pairing energy in the macroscopic energy ${{E}_{\text{mac}}}$ supplemented by a pairing energy correction. The average pairing energy does not contribute to fission (it is independent of the deformation), while the pairing energy correction is treated in a very similar manner to shell effects. The solution to the BCS equations define the pairing correlation energy ${{E}_{\text{pair}}}$. Assuming again a continuous level distribution, one can define a smooth pairing energy ${{\tilde{E}}_{\text{pair}}}$. The difference $\delta {{E}_{\text{pair}}}={{E}_{\text{pair}}}-{{\tilde{E}}_{\text{pair}}}$ defines the pairing correction energy. In most applications to fission, pairing is described within the BCS theory with schematic interactions such as a constant seniority pairing force. Particle number conservation can be accounted for exactly by projection techniques as described in [61] or approximately through the Lipkin–Nogami prescription, see application in [8, 43]. Such a scheme was applied extensively in the fission studies by the LANL/LBNL and Warsaw groups, see, e.g. [48, 49, 62–66] and references therein for some recent work.

The HFB approximation to the energy is centred on the Bogoliubov transformation defining quasiparticle creation and annihilation operators in terms of a given single particle basis $\left({{c}_{i}},c_{i}^{\dagger}\right)$ of a (restricted) Fock space

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\beta}_{\mu}} & = & \underset{i}{\sum}\,U_{i\mu}^{\ast}{{c}_{i}}+\underset{i}{\sum}\,V_{i\mu}^{\ast}c_{i}^{\dagger}, \end{array}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \beta _{\mu}^{\dagger} & = & \underset{i}{\sum}\,{{U}_{i\mu}}c_{i}^{\dagger}+\underset{i}{\sum}\,{{V}_{i\mu}}{{c}_{i}}. \end{array}\end{eqnarray} \tag{ 9 }$$

While the Hilbert space of single particle wave functions is infinite, in practice summations are truncated up to a maximum basis state N. It is convenient to introduce a block matrix notation and matrices of double dimension by writing the equation above as

$$\begin{eqnarray}&&\left(\begin{array}{c} \beta \\ {{\beta}^{\dagger}} \end{array}\right)=\left(\begin{array}{c} {{U}^{+}} & {{V}^{+}} \\ {{V}^{T}} & {{U}^{T}} \end{array}\right)\left(\begin{array}{c} c \\ {{c}^{\dagger}} \end{array}\right)={{W}^{+}}\left(\begin{array}{c} c \\ {{c}^{\dagger}} \end{array}\right),\end{eqnarray} \tag{ 10 }$$

which defines the matrix W of the Bogoliubov transformation

$$\begin{eqnarray}W=\left(\begin{array}{c} U & {{V}^{\ast}} \\ V & {{U}^{\ast}} \end{array}\right).\end{eqnarray} \tag{ 11 }$$

The quasiparticle operators must satisfy canonical fermion commutation relations, which can be summarized by the condition

$$\begin{eqnarray}{{W}^{+}}\sigma W=\sigma,~~~\sigma =\left(\begin{array}{c} 0 & {{I}_{N}} \\ {{I}_{N}} & 0 \end{array}\right),\end{eqnarray} \tag{ 12 }$$

with I_N the N-dimensional identity matrix. The HFB wave function $| \Phi \rangle$ is defined as the vacuum of the quasiparticle annihilation operators, ${{\beta}_{\mu}}| \Phi \rangle =0$ for all μ. This leads to writing $| \Phi \rangle ={\prod}_{\mu}{{\beta}_{\mu}}|0\rangle$ where $|0\rangle$ is the particle vacuum of Fock space and the product runs over all the μ indexes that render $| \Phi \rangle$ non zero.

Densities.

If the ground-state wave function takes the form of a Bogoliubov vacuum, the Wick theorem guarantees that the basic degrees of freedom are the one-body density matrix ρ, the pairing tensor κ and its complex conjugate ${{\kappa}^{\ast}}$ [61, 81]. Each of these objects can be expressed as a function of the Bogoliubov transformation,

$$\begin{eqnarray}&&{{\rho}_{ij}}=\langle \Phi |c_{j}^{\dagger}{{c}_{i}}| \Phi \rangle ={{\left({{V}^{\ast}}{{V}^{T}}\right)}_{ij}},~~~~{{\kappa}_{ij}}=\langle \Phi |{{c}_{j}}{{c}_{i}}| \Phi \rangle ={{\left({{V}^{\ast}}{{U}^{T}}\right)}_{ij}}.\end{eqnarray} \tag{ 13 }$$

The notation can be further condensed by introducing the generalized density matrix,

$$\begin{eqnarray}\mathcal{R}=\left(\begin{array}{c} \rho & \kappa \\ -{{\kappa}^{\ast}} & 1-{{\rho}^{\ast}} \end{array}\right),\end{eqnarray} \tag{ 14 }$$

which is given in terms of the W matrix of the Bogoliubov transformation and the matrix of quasiparticle contractions

$$\begin{eqnarray}\mathbb{R}=\left(\begin{array}{c} \langle \beta _{\mu}^{\dagger}{{\beta}_{\nu}}\rangle & \langle {{\beta}_{\mu}}{{\beta}_{\nu}}\rangle \\ \langle \beta _{\mu}^{\dagger}\beta _{\nu}^{\dagger}\rangle & \langle {{\beta}_{\mu}}\beta _{\nu}^{\dagger}\rangle \end{array}\right)=\left(\begin{array}{c} 0 & 0 \\ 0 & {{I}_{N}} \end{array}\right)\end{eqnarray} \tag{ 15 }$$

as $\mathcal{R}=W\mathbb{R}{{W}^{\dagger}}$. The generalized density matrix encapsulates in a single matrix all degrees of freedom of the theory and simplifies the formal manipulations required in the application of the variational least energy principle.

Variational principle.

Since ρ, κ and ${{\kappa}^{\ast}}$ are the only degrees of freedom, the energy can be expressed most generally as the functional $E\left[\rho,\kappa,{{\kappa}^{\ast}}\right]\equiv E\left[\mathcal{R}\right]$. In practice, one often distinguishes between the particle–hole and particle–particle channel,

$$\begin{eqnarray}&&E\left[\rho,\kappa,{{\kappa}^{\ast}}\right]={{E}^{\text{ph}}}[\rho ]+{{E}^{\text{pp}}}\left[\rho,\kappa,{{\kappa}^{\ast}}\right]\end{eqnarray} \tag{ 16 }$$

although the distinction is a bit artificial since ρ and κ are related through ${{\rho}^{2}}-\rho =-\kappa {{\kappa}^{\dagger}}$. Please also note that pairing interactions depending on the density ρ are very popular in nuclear physics, which adds an extra dependence in ${{E}^{\text{pp}}}$. The actual matrix $\mathcal{R}$ is obtained by minimizing the energy under the condition that the HFB solution remains a quasiparticle vacuum, which is equivalent to ${{\mathcal{R}}^{2}}=\mathcal{R}$. We thus have to minimize

$$\begin{eqnarray}&&\mathcal{E}=E-\text{tr}\left[\Lambda \left({{\mathcal{R}}^{2}}-\mathcal{R}\right)\right],\end{eqnarray} \tag{ 17 }$$

where $\Lambda$ is a matrix of (undetermined) constraints. In this expression and the rest of this paper, 'tr' refers to the trace in the single particle basis (possibly doubled). The variations of the energy are given by

$$\begin{eqnarray}&&\delta \mathcal{E}=\underset{ab}{\sum}\,\frac{\partial E}{\partial {{\mathcal{R}}_{ab}}}\delta {{\mathcal{R}}_{ab}}-\underset{ab}{\sum}\,\frac{\partial}{\partial {{\mathcal{R}}_{ab}}}\text{tr}\left[\Lambda \left({{\mathcal{R}}^{2}}-\mathcal{R}\right)\right]\delta {{\mathcal{R}}_{ab}}.\end{eqnarray} \tag{ 18 }$$

Since $\partial {{\mathcal{R}}_{dc}}/\partial {{\mathcal{R}}_{ab}}={{\delta}_{da}}{{\delta}_{cb}}$, we find after some simple algebra,

$$\begin{eqnarray}&&\underset{ab}{\sum}\,\frac{\partial}{\partial {{\mathcal{R}}_{ab}}}\text{tr}\left[\Lambda \left({{\mathcal{R}}^{2}}-\mathcal{R}\right)\right]=\underset{ab}{\sum}\,{{\left(\mathcal{R} \Lambda + \Lambda \mathcal{R}- \Lambda \right)}_{ba}}.\end{eqnarray} \tag{ 19 }$$

Let us denote $\frac{1}{2}{{\mathcal{H}}_{ba}}=\partial E/\partial {{\mathcal{R}}_{ab}}$. In matrix form, the condition that the variations of the energy should be zero is expressed by the equation

$$\begin{eqnarray}&&\frac{1}{2}\mathcal{H}-\left(\mathcal{R} \Lambda + \Lambda \mathcal{R}- \Lambda \right)=0.\end{eqnarray} \tag{ 20 }$$

Pre- and post-multiplying this equation by $\mathcal{R}$ and subtracting the results after noticing that ${{\mathcal{R}}^{2}}=\mathcal{R}$, we find the final form of the HFB equation

$$\begin{eqnarray}&&\left[\mathcal{H},\mathcal{R}\right]=0.\end{eqnarray} \tag{ 21 }$$

Since $\frac{1}{2}{{\mathcal{H}}_{ba}}=\partial E/\partial {{\mathcal{R}}_{ab}}$, we have by definition $\delta E=$ $\frac{1}{2}\rm{tr}{({\mathcal{H}{\delta}{\mathcal{R}}})}$. Considering the form (14) of the generalized density, we find that the matrix of $\mathcal{H}$ reads

$$\begin{eqnarray}\mathcal{H}=\left(\begin{array}{c} h & \Delta \\ -{{ \Delta }^{\ast}} & -{{h}^{\ast}} \end{array}\right),\end{eqnarray} \tag{ 22 }$$

with the mean field h and pairing field $\Delta$ defined by

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{h}_{ij}}=\frac{\partial E}{\partial {{\rho}_{ji}}}, & {} & {{ \Delta }_{ij}}=\frac{\partial E}{\partial \kappa _{ij}^{\ast}}, \\ \Delta _{ij}^{\ast}=\frac{\partial E}{\partial {{\kappa}_{ij}}}, & {} & h_{ij}^{\ast}=\frac{\partial E}{\partial \rho _{ji}^{\ast}}. \end{array}\end{eqnarray} \tag{ 23 }$$

The form of the HFB equation (21) is not well suited for its practical solution and, therefore, it is often reinterpreted by considering that $\left[\mathcal{H},\mathcal{R}\right]=0$ implies that $\mathcal{H}$ and $\mathcal{R}$ can be diagonalized by the same Bogoliubov transformation. Since the generalized density matrix is diagonalized by W as shown in (15), the same must hold for $\mathcal{H}$,

$$\begin{eqnarray}&&\mathcal{H}W=W\mathcal{E}.\end{eqnarray} \tag{ 24 }$$

This represents a non-linear diagonalization problem (since $\mathcal{H}$ depends upon W through the densities) with eigenvalues $\mathcal{E}$ and eigenvectors W. It turns out that the matrix of eigenvalues can be written schematically

$$\begin{eqnarray}\mathcal{E}=\left(\begin{array}{c} {{E}_{\mu}} & 0 \\ 0 & -{{E}_{\mu}} \end{array}\right),\end{eqnarray} \tag{ 25 }$$

that is, for each positive eigenvalue ${{E}_{\mu}}$, there is another one of opposite sign $-{{E}_{\mu}}$. The eigenvalues of the HFB matrix are referred to as 'quasiparticle energies'. In section 4.1 p 33, we discuss the technical aspects of solving such a non-linear eigenvalue problem. The quasiparticle energies ${{E}_{\mu}}$ defined above should not be confused with the single particle energies typical of the phenomenological mean field potentials of the MM method. The equivalent to the single particle energies in the HFB case are the eigenvalues of the Hartree–Fock (HF) Hamiltonian h of (23). The study of those quantities allows obtaining the same degree of insight that can be obtained by the study of the Nilsson orbitals.

Energy and fields.

In nuclear physics, the energy functional $E\left[\rho,\kappa,{{\kappa}^{\ast}}\right]$ is typically composed of two parts: one is derived from an effective two-body Hamiltonian,

$$\begin{eqnarray}&&\hat{H}=\underset{ij}{\sum}\,{{t}_{ij}}c_{i}^{\dagger}{{c}_{j}}+\frac{1}{4}\underset{ijkl}{\sum}\,{{\bar{v}}_{ijkl}}c_{i}^{\dagger}c_{j}^{\dagger}{{c}_{l}}{{c}_{k}},\end{eqnarray} \tag{ 26 }$$

where t_ij refers to the matrix elements of the one-body kinetic energy operator and ${{\bar{v}}_{ijkl}}$ to the anti-symmetrized matrix elements of the two-body potential. The other part of the energy functional comes in the form of various phenomenological density dependent terms that have to be handled with some care as they introduce additional terms in the HF Hamiltonian and pairing fields that come in the form of rearrangement potentials (see below). In the case of such an effective two-body Hamiltonian, the HFB energy is simply given by $E=\langle \Phi |\hat{H}| \Phi \rangle$, with $| \Phi \rangle ={\prod}_{\mu}{{\beta}_{\mu}}|0\rangle$ the quasiparticle vacuum introduced earlier. Considering the definitions (13) for the density matrix and pairing tensor, the application of the Wick theorem gives

$$\begin{eqnarray}&&E=\text{tr}(t\rho )+\frac{1}{2}\text{tr}(\Gamma \rho )-\frac{1}{2}\text{tr}\left(\Delta {{\kappa}^{\ast}}\right).\end{eqnarray} \tag{ 27 }$$

In this particular case, the mean field potential and the pairing field are given as a function of the two-body matrix elements by

$$\begin{eqnarray}&&{{ \Gamma }_{ik}}=\underset{jl}{\sum}\,{{\bar{v}}_{ijkl}}\;{{\rho}_{lj}},~~~~~~~{{ \Delta }_{ij}}=\frac{1}{2}\underset{kl}{\sum}\,{{\bar{v}}_{ijkl}}\;{{\kappa}_{kl}}.\end{eqnarray} \tag{ 28 }$$

and the HF Hamiltonian reads $h=t+ \Gamma$. The notation can be further condensed by defining the generalized kinetic energy and mean field matrices by

$$\begin{eqnarray}\mathcal{T}=\left(\begin{array}{c} t & 0 \\ 0 & -{{t}^{\ast}} \end{array}\right),~~~~~~\mathcal{K}=\left(\begin{array}{c} \Gamma & \Delta \\ -{{ \Delta }^{\ast}} & -{{ \Gamma }^{\ast}} \end{array}\right).\end{eqnarray} \tag{ 29 }$$

With these notations, we note that $\mathcal{H}\left[\mathcal{R}\right]=\mathcal{T}+\mathcal{K}\left[\mathcal{R}\right]$. The total energy can then be written in the very compact form

$$\begin{eqnarray}&&E=\frac{1}{4}\text{tr}\left[\left(\mathcal{H}+\mathcal{T}\,\,\right)\mathcal{S}\right],\end{eqnarray} \tag{ 30 }$$

with

$$\begin{eqnarray}\mathcal{S}=\left(\begin{array}{c} \rho & \kappa \\ -{{\kappa}^{\ast}} & -{{\rho}^{\ast}} \end{array}\right)=\mathcal{R}-\left(\begin{array}{c} 0 & 0 \\ 0 & {{I}_{N}} \end{array}\right).\end{eqnarray} \tag{ 31 }$$

Density-dependent interactions.

In the practical implementation of the HFB theory in nuclear structure, special attention should be paid to the very common case of two-body, density-dependent effective 'forces' ${{\hat{v}}^{\text{DD}}}\left(\boldsymbol{r},\boldsymbol{r}^{{\prime}}\right)={{\hat{v}}^{\text{DD}}}\left[\rho \left(\boldsymbol{R}\right)\right]\delta \left(\boldsymbol{r}-\boldsymbol{r}^{{\prime}}\right)$ with $\boldsymbol{R}=\left(\boldsymbol{r}+\boldsymbol{r}^{{\prime}}\right)/2$. Note that in this case, the potential ${{\hat{v}}^{\text{DD}}}$ cannot be put into strict second quantized form—as emphasized in [82]. However, one can still define an energy functional $E\left[\rho,\kappa,{{\kappa}^{\ast}}\right]$ by taking the expectation value of ${{\hat{v}}^{\text{DD}}}$ in coordinate space on a quasiparticle vacuum. Given the energy functional, the mean field and pairing field can then be defined from (23). When computing these derivatives with respect to the density matrix, additional contributions related to the matrix elements of $\partial {{\hat{v}}^{\text{DD}}}/\partial \rho$ will arise (rearrangement terms). They have the generic form

$$\begin{eqnarray}&&{{ \Gamma }_{ik}}\to {{ \Gamma }_{ik}}+\underset{jl}{\sum}\,\partial v_{ijkl}^{\text{DD}}{{\rho}_{lj}},\end{eqnarray} \tag{ 32 }$$

with

$$\begin{eqnarray}&&\partial v_{ijkl}^{\text{DD}}={\int}^{}{{\text{d}}^{3}}\boldsymbol{r}{\int}^{}{{\text{d}}^{3}}\boldsymbol{r}^{{\prime}}\;\varphi _{i}^{\ast}(\boldsymbol{r})\varphi _{j}^{\ast}(\boldsymbol{r}^{{\prime}})\frac{\partial {{{\hat{v}}}^{\text{DD}}}}{\partial \rho}(\boldsymbol{R})\delta (\boldsymbol{r}-\boldsymbol{r}^{{\prime}}){{\varphi}_{k}}(\boldsymbol{r}){{\varphi}_{l}}(\boldsymbol{r}^{{\prime}}).\end{eqnarray} \tag{ 33 }$$

Note that one could also consider pairing-tensor-dependent potentials in the same way. Recently, an extension of the above scheme so as to consider finite range density-dependent forces has been proposed [83]. In this case, the simple formulas above have to be adapted but the same conceptual procedure remains valid.

Representation based on the Thouless theorem.

There is an alternative way to obtain and represent the HFB equation, which is based on the Thouless theorem of [84]. We recall that the Thouless theorem presented, e.g. in [61, 81], establishes that two non-orthogonal quasiparticle vacua $|{{ \Phi }_{0}}\rangle$ and $|{{ \Phi }_{1}}\rangle$ (corresponding to two different Bogoliubov transformations W₀ and W₁) are related through

$$\begin{eqnarray}&&|{{ \Phi }_{1}}\rangle =\exp \left(\text{i}\hat{Z}\right)|{{ \Phi }_{0}}\rangle,\end{eqnarray} \tag{ 34 }$$

where $\hat{Z}$ is an hermitian one-body operator written in the quasiparticle basis corresponding to $|{{ \Phi }_{0}}\rangle$ as

$$\begin{eqnarray}&&\hat{Z}=\underset{\mu \nu}{\sum}\,Z_{\mu \nu}^{11}\beta _{\mu}^{\dagger}{{\beta}_{\nu}}+\frac{1}{2}\underset{\mu \nu}{\sum}\,Z_{\mu \nu}^{20}\beta _{\mu}^{\dagger}\beta _{\nu}^{\dagger}+\text{h}.\text{c}.\end{eqnarray} \tag{ 35 }$$

The complex matrices Z¹¹ (hermitian) and Z²⁰ (skew-symmetric) are determined by the requirement

$$\begin{eqnarray}{{W}_{1}}={{W}_{0}}\exp \left[\text{i}\left(\begin{array}{c} {{Z}^{11}} & {{Z}^{20}} \\ -{{Z}^{20\,\ast}} & -{{Z}^{11\,\ast}} \end{array}\right)\right]={{W}_{0}}\exp \left[\text{i}Z\right].\end{eqnarray} \tag{ 36 }$$

The relation above can be used to generate a Bogoliubov transformation W₁ given an initial one W₀ and a set of Z¹¹ and Z²⁰ coefficients. Not surprisingly, the number of free parameters in Z¹¹ and Z²⁰ is the same as in a general Bogoliubov transformation W after taking into account the condition ${{W}^{+}}\sigma W=\sigma$. Therefore, these matrices can be used as independent variational parameters. Starting from an initial Bogoliubov transformation W and considering infinitesimal variations $\delta W$ characterized by the matrix Z (that is, Z¹¹ and Z²⁰), we obtain to first order in Z an explicit expression for the variation of the Bogoliubov amplitude W as $\delta W=\text{i}WZ$. It follows that the variation $\delta \mathcal{R}$ of the generalized density $\mathcal{R}=W\mathbb{R}{{W}^{\dagger}}$ under small variations $\delta W$ can be expressed in terms of the Z matrix as $\delta \mathcal{R}=\text{i}\left(\mathcal{Z}\mathcal{R}-\mathcal{R}\mathcal{Z}\right)$ with $\mathcal{Z}=WZ{{W}^{\dagger}}$. From $\delta E=\frac{1}{2}\text{tr}\left(\mathcal{H}\delta \mathcal{R}\right)$ for the variation of the energy, we thus find $\delta E=\frac{\text{i}}{2}\text{tr}\left(\left[\mathcal{R},\mathcal{H}\right]\mathcal{Z}\right)$. The variational principle condition $\delta E=0$ yields the HFB equation $\left[\mathcal{R},\mathcal{H}\right]=0$. Finally, calculating the variation of the energy $\delta E$ to second order in $\mathcal{Z}$ tell us that, if $\mathcal{Z}$ is chosen as $\text{i}\eta {{\left[\mathcal{R},\mathcal{H}\right]}^{\dagger}}$, where η is a sufficiently small step size, then $\delta E$ is negative. This result represents the essence of the gradient method commonly used in numerical implementations of the HFB equation and discussed in more detail in section 4.1.3; see also [85] for a complete presentation.

Constraints.

By construction, the quasiparticle vacuum is not an eigenstate of the particle number operator. In practice, it is thus always necessary to impose a constraint on the average value of the particle number operator, both for protons and neutrons. Solving the HFB equation (21) subject to constraints is also particularly important in fission studies. It is required in the evaluation of the potential energy surface as a function of the collective coordinates used to characterize the fission process. The handling of constraints in the formalism outlined above is straightforward since it only requires replacing the HFB matrix by the auxiliary operator ${{\mathcal{H}}_{ij}}\to {{\mathcal{H}}_{ij}}-{\sum}^{}{{\lambda}_{\alpha}}{{\mathcal{O}}_{\alpha,ij}}$ where ${{\mathcal{O}}_{\alpha,ij}}$ are the matrix elements of the constraint operators ${{\hat{O}}_{\alpha}}$ in the doubled basis and ${{\lambda}_{\alpha}}$ are Lagrange multipliers. One-body operators such as, e.g. the particle number operator or multipole moments, have the same generic structure as in (29). For example, solving the HFB equation with a constraint on particle number implies that the HFB matrix becomes

$$\begin{eqnarray}\mathcal{H}=\left(\begin{array}{c} h-\lambda & \Delta \\ -{{ \Delta }^{\ast}} & -{{h}^{\ast}}+\lambda \end{array}\right),\end{eqnarray} \tag{ 37 }$$

with λ the Fermi energy. In the previous equation, $h-\lambda$ is a shorthand notation for ${{h}_{ij}}-\lambda {{\delta}_{ij}}$. Nearly all operators needed in fission are one-body operators, and their expectation value is thus given by $\langle {{\hat{O}}_{\alpha}}\rangle =\text{tr}\left[{{\hat{O}}_{\alpha}}\rho \right]$. The figure 5 shows an example of constrained HFB solutions in the particular case of a single collective variable. In the figure, the quadrupole deformation ${{\beta}_{2}}$ is obtained from the axial quadrupole moment ${{\hat{Q}}_{20}}$ used as a constraint through ${{\beta}_{2}}=\sqrt{5/(16\pi )}4\pi /\left(3r_{0}^{2}{{A}^{5/3}}\right)\langle {{\hat{Q}}_{20}}\rangle$. Section 2.3 discusses in more detail typical collective variables used in fission calculations.

Zoom In Zoom Out Reset image size

To solve the non-linear HFB equation several approaches are used. A very popular one is the iterative method, where the density matrix of the nth iteration is used to compute the $\mathcal{H}$ matrix for the (n +1)th one. Diagonalization of this matrix produces a new density matrix and the process is repeated until the generalized density matrix remains stationary up to the desired precision. Another popular method is based on the direct minimization of the energy using any of the variants of the gradient method [86]. Finally, the imaginary time method is also popular in the particular case of the HF +BCS equation. The advantages and disadvantages of the three of them will be discussed in more details in section 4.1.3.

As immediately visible from the form (37) of the HFB matrix, solving the HFB equations requires handling matrices twice larger than the simpler HF theory. In addition, the HFB spectrum is unbounded from below and from above, which leads to practical difficulties when trying to solve the HFB equation on a lattice—see section 4.1.2. These are some of the reasons why the BCS approximation to the HFB equation is sometimes preferred in applications.

The idea is to first write the HFB matrix $\mathcal{H}$ in a special basis—the single particle basis where the HF Hamiltonian h is diagonal with eigenvalues ${{e}_{\mu}}$ (the HF single particle energies). In other words, we find the transformation matrix D such that hD  =  De where e is the matrix of eigenvalues ${{e}_{\mu}}$. The BCS approximation then consists in imposing that the skew-symmetric tensor ${{ \Delta }_{\mu {{\mu}^{\prime}}}}$ of (23) be in normal form in that basis. Specifically, we impose that it can only connect states with ${{\mu}^{\prime}}=\bar{\mu}$, where the single particle state $\bar{\mu}$ is the partner of state μ under the time-reversal operator: ${{ \Delta }_{\mu {{\mu}^{\prime}}}}\sim {{ \Delta }_{\mu \bar{\mu}}}{{\delta}_{\bar{\mu}{{\mu}^{\prime}}}}$. If we denote by $\Delta$ the diagonal matrix with elements ${{ \Delta }_{\mu}}\equiv {{ \Delta }_{\mu \bar{\mu}}}$, we find that the HFB matrix (37) becomes

$$\begin{eqnarray}{{\mathcal{H}}_{\text{BCS}}}=\left(\begin{array}{c} e-\lambda & 0 & 0 & \Delta \\ 0 & e-\lambda & - \Delta & 0 \\ 0 & -{{ \Delta }^{\ast}} & -e+\lambda & 0 \\ {{ \Delta }^{\ast}} & 0 & 0 & -e+\lambda \end{array}\right).\end{eqnarray} \tag{ 38 }$$

This matrix can easily be rearranged and block-diagonalized by a BCS transformation, which is equivalent to solving the traditional 'gap equation' of the BCS theory,

$$\begin{eqnarray}&&\left({{e}_{\mu}}-\lambda \right)\left(u_{\mu}^{2}-v_{\mu}^{2}\right)+2{{ \Delta }_{\mu}}{{u}_{\mu}}{{v}_{\mu}}=0.\end{eqnarray} \tag{ 39 }$$

The size of the matrix to be diagonalized, h, is reduced by a factor of two as compared to HFB, which in turns lowers the computational cost by a factor of eight. Several test calculations, for example [87, 88], have shown that the HF +BCS approach is a reasonable approximation to the full HFB treatment for fission calculations. This conclusion does not hold any more for weakly-bound nuclei close to the neutron dripline, where the BCS approximation induces a coupling with a non-physical gas of neutrons as discussed in [89].

There is a relation between the quasiparticle energies ${{E}_{\mu}}$ of the HFB method and the single particle energies ${{e}_{\mu}}$, the Fermi energy λ and the matrix element of the pairing field ${{ \Delta }_{\mu}}$, ${{E}_{\mu}}=\sqrt{{{\left({{e}_{\mu}}-\lambda \right)}^{2}}+ \Delta _{\mu}^{2}}$. It follows that the coefficients of the BCS transformation ${{u}_{\mu}}$ and ${{v}_{\mu}}$ acquire the very simple form

$$\begin{eqnarray}&&\left. \begin{array}{c} u_{\mu}^{2} \\ v_{\mu}^{2} \end{array}\right\}=\frac{1}{2}\left[1\pm \frac{{{e}_{\mu}}-\lambda}{\sqrt{{{\left({{e}_{\mu}}-\lambda \right)}^{2}}+ \Delta _{\mu}^{2}}}\right].\end{eqnarray} \tag{ 40 }$$

In the case of a pure pairing force with constant matrix elements  −G the ${{ \Delta }_{\mu {{\mu}^{\prime}}}}$ matrix is already in normal form with constant ${{ \Delta }_{\mu}}$ matrix elements ${{ \Delta }_{\mu}}= \Delta =G{\sum}_{\nu >0}{{u}_{\nu}}{{v}_{\nu}}$. The ${{v}_{\nu}}$ and ${{u}_{\nu}}$ are again the coefficients of the BCS transformation. In this case, however, the attractive pairing interaction can only be active within a restricted window of single particle levels, being zero elsewhere.

As already emphasized, most nuclear energy functionals used in fission studies are explicitly derived from the expectation value of effective, density-dependent, local two-body nuclear potentials at the HFB approximation. For the particle–hole channel, that is, the component ${{E}^{\text{ph}}}[\rho ]$ of the full functional (16), the two most popular functionals are based on the Skyrme and Gogny two-body effective potentials. Both EDFs are intended to be applicable throughout the nuclear chart using a limited set of parameters adjusted to global nuclear properties.

The Skyrme EDF is the energy functional $E[\rho ]$ derived from the expectation value of the Skyrme potential introduced in [90] on a Slater determinant. The current standard form of the two-body Skyrme interaction is

$$\begin{eqnarray}&&{{\hat{V}}_{12}}\left(\boldsymbol{r}_{{1}},\boldsymbol{r}_{{2}}\right)={{t}_{0}}\left(1+{{x}_{0}}{{\hat{P}}_{\sigma}}\right)\delta \left(\boldsymbol{r}_{{1}}-\boldsymbol{r}_{{2}}\right)\\ &&+\frac{1}{2}{{t}_{1}}\left(1+{{x}_{1}}{{\hat{P}}_{\sigma}}\right)\left[{{\widehat{\boldsymbol{k}^{{\prime}}}}^{2}}\delta \left(\boldsymbol{r}_{{1}}-\boldsymbol{r}_{{2}}\right)+\delta \left(\boldsymbol{r}_{{1}}-\boldsymbol{r}_{{2}}\right){{\widehat{\boldsymbol{k}}}^{2}}\right]\\ &&+{{t}_{2}}\left(1+{{x}_{2}}{{\hat{P}}_{\sigma}}\right)\widehat{\boldsymbol{k}^{{\prime}}}\centerdot \delta \left(\boldsymbol{r}_{{1}}-\boldsymbol{r}_{{2}}\right)\widehat{\boldsymbol{k}}\\ &&+\text{i}{{W}_{0}}\widehat{\boldsymbol{S}}\centerdot \left[\widehat{\boldsymbol{k}^{{\prime}}}\wedge \delta \left(\boldsymbol{r}_{{1}}-\boldsymbol{r}_{{2}}\right)\widehat{\boldsymbol{k}}\right]\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & + & \frac{1}{6}{{t}_{3}}\left(1+{{x}_{3}}{{{\hat{P}}}_{\sigma}}\right){{\rho}^{\alpha}}\left(\boldsymbol{R}_{{12}}\right)\delta \left(\boldsymbol{r}_{{1}}-\boldsymbol{r}_{{2}}\right), \end{array}\end{eqnarray} \tag{ 42 }$$

where

$$\begin{eqnarray}&&\widehat{\boldsymbol{k}}=-\frac{\text{i}}{2}\left(\boldsymbol{}{{\nabla}_{1}}-\boldsymbol{}{{\nabla}_{2}}\right),~~~\widehat{\boldsymbol{k}^{{\prime}}}=\frac{\text{i}}{2}\left({{\overleftarrow{\boldsymbol{}\nabla}}_{1}}-{{\overleftarrow{\boldsymbol{}\nabla}}_{2}}\right)\end{eqnarray} \tag{ 43 }$$

are the relative momentum operator, $\boldsymbol{R}_{{12}}=\left(\boldsymbol{r}_{{1}}+\boldsymbol{r}_{{2}}\right)/2$, ${{\hat{P}}_{\sigma}}$ is the spin exchange operator and $\widehat{\boldsymbol{S}}=\boldsymbol{\sigma }_{{1}}+\boldsymbol{\sigma }_{{2}}$ is the total spin. By convention, $\widehat{\boldsymbol{k}^{{\prime}}}$ acts on the left. The Skyrme potential contains a phenomenological density-dependent term (term proportional to t₃). Originally this term was introduced with $\alpha =1$ and x₃  =  1 in [91] to simulate the effect of three-body forces, in particular with respect to the saturation of nuclear forces. Following [92], both α and x₃ are now usually taken as phenomenological adjustable parameters. The expectation value of the potential (42) on a Slater determinant can be expressed as a functional of the local density $\rho \left(\boldsymbol{r}\right)$ and various other local densities derived from $\rho \left(\boldsymbol{r}\right)$. The potential energy is then given by

$$\begin{eqnarray}&&{{E}_{\text{Skyrme}}}={\int}^{}{{\text{d}}^{3}}\boldsymbol{r}\underset{t=0,1}{\sum}\,\left[\mathcal{H}_{t}^{\text{even}}\left(\boldsymbol{r}\right)+\mathcal{H}_{t}^{\text{odd}}\left(\boldsymbol{r}\right)\right],\end{eqnarray} \tag{ 44 }$$

where t  =  0 refers to isoscalar energy densities, and t  =  1 to isovector ones. The contribution $\mathcal{H}_{t}^{\text{even}}\left(\boldsymbol{r}\right)$ to the total energy density depends only on time-even densities,

$$\begin{eqnarray}&&\mathcal{H}_{t}^{\text{even}}\left(\boldsymbol{r}\right)=C_{t}^{\rho}\rho _{t}^{2}+C_{t}^{ \Delta \rho}{{\rho}_{t}} \Delta {{\rho}_{t}}+C_{t}^{\tau}{{\rho}_{t}}{{\tau}_{t}}+C_{t}^{J}{\mathsf {J}}_{t}^{2}+C_{t}^{\nabla J}{{\rho}_{t}}\boldsymbol{}\nabla \centerdot \boldsymbol{J}_{{t}},\end{eqnarray} \tag{ 45 }$$

while the contribution $\mathcal{H}_{t}^{\text{odd}}\left(\boldsymbol{r}\right)$ depends only on time-odd densities,

$$\begin{eqnarray}&&\mathcal{H}_{t}^{\text{odd}}\left(\boldsymbol{r}\right)=C_{t}^{s}\boldsymbol{s}_{t}^{2}+C_{t}^{ \Delta s}\boldsymbol{s}_{{t}}\centerdot \Delta \boldsymbol{s}_{{t}}+C_{t}^{T}\boldsymbol{s}_{{t}}\centerdot \boldsymbol{T}_{{t}}+C_{t}^{j}\boldsymbol{j}_{t}^{2}+C_{t}^{\nabla j}\boldsymbol{s}_{{t}}\centerdot \left(\boldsymbol{}\nabla \wedge \boldsymbol{j}_{{t}}\right).\end{eqnarray} \tag{ 46 }$$

The various densities involved in these expressions are: the kinetic energy density ${{\tau}_{t}}$; the spin current density ${{{\mathsf {J}}}_{t}}$ (rank 2 tensor); the vector part of the spin current density $\boldsymbol{J}_{{t}}$ (which is obtained by the tensor contraction $\boldsymbol{J}_{{\kappa,t}}={\sum}_{\mu \nu}{{\epsilon}_{\mu \nu \kappa}}{{\text{J}}_{\mu \nu,t}}$); the spin density $\boldsymbol{s}_{{t}}$; the spin kinetic energy density $\boldsymbol{T}_{{t}}$ and the momentum density $\boldsymbol{j}_{{t}}$. Their actual expressions as a function of the one-body density matrix can be found in [79, 93, 94]. The Skyrme energy density is the most general scalar that can be formed by coupling the fields derived from the one-body density matrix up to second order in derivatives as discussed in [95, 96]. There have been attempts to generalize the EDF (45) and (46) by going beyond the second order derivative [97–100] or by adding local, zero-range three-body forces to the potential given by (42) [101].

The Gogny force is the most widely used non-relativistic effective finite-range nuclear potential. It can be viewed as a Skyrme potential where the zero range central potential is replaced by the sum of two Gaussians in the spatial part. In fact, the central part of the Skyrme force can be obtained by expanding a Gaussian central and spin–orbit potential up to second order in momentum space as shown in [91]. The main advantage of the finite range is that the matrix elements are free from ultraviolet divergences in the particle–particle channel, which allows using the same potential in both particle–hole and particle–particle channels without introducing a window of active particles in the particle–particle channel. In its original formulation of [102], the Gogny force reads

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\hat{V}}}_{12}}\left(\boldsymbol{r}_{{1}},\boldsymbol{r}_{{2}}\right) & = & \underset{i=1}{\overset{2}{\sum}}\,\,\left({{W}_{i}}+{{B}_{i}}{{{\hat{P}}}_{\sigma}}-{{H}_{i}}{{{\hat{P}}}_{\tau}}-{{M}_{i}}{{{\hat{P}}}_{\sigma}}{{{\hat{P}}}_{\tau}}\right)\,{{\text{e}}^{-\frac{{{\left(\boldsymbol{r}_{{1}}-\boldsymbol{r}_{{2}}\right)}^{2}}}{\mu _{i}^{2}}}} \\ {} & + & \text{i}{{W}_{\text{LS}}}\widehat{\boldsymbol{S}}\centerdot \left[\widehat{\boldsymbol{k}^{{\prime}}}\wedge \delta \left(\boldsymbol{r}_{{1}}-\boldsymbol{r}_{{2}}\right)\widehat{\boldsymbol{k}}\right] \\ {} & + & {{t}_{3}}\left(1+{{x}_{0}}{{{\hat{P}}}_{\sigma}}\right)\,{{\rho}^{\alpha}}\left(\boldsymbol{R}_{{12}}\right)\delta \left(\boldsymbol{r}_{{1}}-\boldsymbol{r}_{{2}}\right). \end{array}\end{eqnarray} \tag{ 47 }$$

Note that the energy of two-body finite-range potentials such as the Gogny could also be put in the form (44), the difference being that the energy densities $\mathcal{H}\left(\boldsymbol{r}\right)$ would be expressed themselves not directly as functionals of the local density $\rho \left(\boldsymbol{r}\right)$, but as integrals over $\boldsymbol{r}^{{\prime}}$ involving the non-local density $\rho \left(\boldsymbol{r},\boldsymbol{r}^{{\prime}}\right)$.

For both the Skyrme and Gogny forces, the energy density functional is obtained from a central and spin orbit potential plus a phenomenological density-dependent term. Note that we have not discussed the inclusion of explicit tensor potentials in either the Skyrme or Gogny forces, see [94, 103, 104] for detailed discussions on that topic. In the spirit of the Kohn–Sham theory, however, the EDF does not need to be derived from any potential and could be parametrized directly as a functional of the density. In the Barcelona–Catania–Paris–Madrid (BCPM) functional of [105, 106], this idea is only applied to the volume part of the functional, as it can easily be constrained by nuclear and neutron matter properties. The resulting functional of the density is supplemented with a term from a spin–orbit potential and another potential providing the surface term, leading to

$$\begin{eqnarray}&&E[\rho ]=T+{{E}^{\text{SO}}}[\rho ]+E_{\text{int}}^{\infty}[\rho ]+E_{\text{int}}^{\text{FR}}[\rho ],\end{eqnarray} \tag{ 48 }$$

where T is the kinetic energy, ${{E}^{\text{SO}}}$ is the spin–orbit energy density obtained from a zero range spin–orbit potential, $E_{\text{int}}^{\infty}$ is the bulk part given in terms of a fitting polynomial of the density with parameters adjusted to reproduce nuclear matter properties, and $E_{\text{int}}^{\text{FR}}$ is the surface energy obtained from a finite-range Gaussian potential. Similar ideas had also been pursued earlier by Fayans and collaborators in [107, 108]. In their work, the parametrization of $E_{\text{int}}^{\infty}$ was different from the BPCM functional, and the $E_{\text{int}}^{\text{FR}}$ term was derived from a Yukawa potential instead of a Gaussian.

In realistic calculations, the Coulomb potential must also be included for protons. The direct contribution of this potential to the energy does not pose any particular problem. The exchange term is usually treated in the Slater approximation although several studies such as [109] have shown the impact of both Coulomb exchange and the associated Coulomb anti-pairing effect in collective inertias.

As mentioned in section 2.2.1, the nuclear EDF comprises two components, ${{E}^{\text{ph}}}[\rho ]$ in the p.h. channel and ${{E}^{\text{pp}}}\left[\rho,\kappa,{{\kappa}^{\ast}}\right]$ in the p.p. channel. Once again, pairing functionals are most often obtained by taking the expectation value of an effective (two-body) potential on the quasiparticle vacuum. Users of the Gogny force often take the same potential for the pairing channel. In the case of the Skyrme EDF, it is customary to consider simple pairing forces that can be put into functionals of the local pairing density $\tilde{\rho}\left(\boldsymbol{r}\right)$ introduced in [110]. The full one-body pairing density $\tilde{\rho}\left(\boldsymbol{r}\sigma,\boldsymbol{r}^{{\prime}}{{\sigma}^{\prime}}\right)$, depending on spatial ($\boldsymbol{r}$) and spin (σ) coordinates, is related to the pairing tensor through $\kappa \left(\boldsymbol{r}\sigma,\boldsymbol{r}^{{\prime}}{{\sigma}^{\prime}}\right)=2{{\sigma}^{\prime}}\tilde{\rho}\left(\boldsymbol{r}\sigma,\boldsymbol{r}^{{\prime}}-{{\sigma}^{\prime}}\right)$ ($\sigma,{{\sigma}^{\prime}}=\pm 1/2$) and has similar symmetry properties to the one-body density matrix ρ. A commonly used pairing force is the density-dependent, zero range potential

$$\begin{eqnarray}&&{{\hat{V}}_{\text{pair}}}\left(\boldsymbol{r}_{{1}},\boldsymbol{r}_{{2}}\right)=\underset{t=n,p}{\sum}\,V_{t}^{0}\left[1-V_{t}^{1}{{\left(\frac{{{\rho}_{0}}\left(\boldsymbol{R}_{{12}}\right)}{{{\rho}_{c}}}\right)}^{\alpha}}\right]\delta \left(\boldsymbol{r}_{{1}}-\boldsymbol{r}_{{2}}\right),\end{eqnarray} \tag{ 49 }$$

where $V_{t}^{0}$, $V_{t}^{1}$, ${{\rho}_{c}}$, α are adjustable parameters and ${{\rho}_{0}}$ the isoscalar one-body density. While $V_{t}^{0}$ controls the strength of the pairing interaction, ${{\rho}_{c}}$ represents the average density inside the nucleus (it is often set at 0.16 fm⁻³) and, therefore, $V_{t}^{1}$ controls the surface or volume nature of the pairing interaction. Such a schematic pairing force was originally introduced in [111]. Note that its zero-range requires introducing a cut-off in the number of active particles used to define the densities. In both the BCPM and Fayans functionals, a similar zero-range density-dependent interaction was adopted [107, 112].

Irrespective of the mathematical form of the functional, nuclear EDFs contain free parameters that must be adjusted to experimental data. In the case of the Skyrme force, these are the (t, x) parameters plus W₀ and α; for Skyrme EDFs, the ${{C}^{u{{u}^{\prime}}}}$ coupling constants (see [113] for an alternative representation of the Skyrme EDF where volume coupling constants are expressed as a function of nuclear matter characteristics); in the case of the Gogny force, the parameters are the ${{W}_{i}},{{B}_{i}},{{H}_{i}},{{M}_{i}}$ parameters of the central part, plus W_LS and x₀. Fitting these parameters on select nuclear data is an example of an inverse problem in statistics, and several strategies have been explored in the past; see [114, 115] for a discussion.

In the case of the Skyrme EDF, there are hundreds of different parametrizations; see [116] for a review. However, not all of these parametrizations are suitable for fission studies, where the ability of the EDF to reproduce deformation properties is essential. The two most widely used Skyrme EDFs to date are the SkM* parametrization of Bartel et al [117] and the UNEDF1 parametrization of Kortelainen et al [118]. In both cases, the functionals were fitted by considering fission data in actinides, namely the height of the fission barrier in ²⁴⁰Pu for SkM* and the value of a few fission isomer excitation energies for UNEDF1. The case of the Gogny force is similar, although the number of Gogny parametrizations is much more limited since only five parametrizations have been published so far [102, 119–122]. Just as for the Skyrme force, the original D1 and D1' parametrizations of the Gogny force could not reproduce accurately enough the height of the fission barrier in ²⁴⁰Pu. The D1S parametrization proposed by Berger in [123] was a slightly modified version of D1 where the surface energy coefficient in nuclear matter was modified to decrease the fission barrier. Since then, the D1S parametrization has been the most popular force for fission. Note that there are extensions of the Gogny force to include density dependencies in each term as in [83] or to add a tensor force with a finite-range spatial part, see e.g. [104, 124]. In the BCPM case there is a single set of parameters, adjusted to nuclear matter properties (bulk part) and to the binding energies of even–even nuclei (surface term). Fission information has not been included in the fit.

As mentioned in the introduction, applications of the neutron-induced fission of actinides with fast neutrons (kinetic energies of the order of 14 MeV) involve excitation energies of the fissioning nucleus that can be as large as 20 MeV or more. At such excitation energies, the potential energy surface of the nucleus can be markedly different from what is obtained from constrained HF +BCS or HFB calculations. Similar and higher excitation energies can be reached during the formation of the compound nucleus in heavy ion fusion reactions, which are one of the primary mechanisms to synthesise heavy elements as discussed in [125]. As shown in [56], the level density increases exponentially with E^*, and at excitation energies beyond a few MeV, it becomes impossible to individually track excited states using direct methods. Finite-temperature density functional theory thus provides a convenient toolbox to quantify the evolution of nuclear deformation properties as a function of excitation energy. The inclusion of temperature in microscopic studies of fission has been done within the Hartree–Fock approach (FT-HF) in [126, 127], the HF +BCS approach in [17] and the fully-fledged HFB approach (FT-HFB) in [16, 17, 128, 129]. The reader can also refer to [129] for additional references to the use of temperature in semi-microscopic methods. Note that in practical applications, the temperature must be related to the excitation energy of the nucleus, which is not entirely trivial; see [129] for a detailed discussion.

Various elements of the FT-HF and FT-HFB theory can be found in [81, 130–136]. Let us recall that the nucleus is assumed to be in a state of thermal equilibrium at temperature T. It is then characterized by a statistical density operator $\hat{D}$, which contains all information on the system. In particular, once the density operator is known, the expectation value of any operator $\hat{F}$ can be computed by taking the trace $\text{Tr}\,\hat{D}\hat{F}$ in any basis of the Fock space. Here, the notation 'Tr' refers to the statistical trace by contrast to the trace 'tr' used before with respect to a single particle basis. If we further assume the system can be described by a grand canonical ensemble (the average value of the energy and particle numbers are fixed), then the application of the principle of maximum entropy yields the following generic form of the density operator as

$$\begin{eqnarray}&&\hat{D}=\frac{1}{Z}{{\text{e}}^{-\beta \left(\hat{H}-\lambda \hat{N}\right)}},\end{eqnarray} \tag{ 50 }$$

where Z is the grand partition function, $Z=\text{Tr}\,{{\text{e}}^{-\beta \left(\hat{H}-\lambda \hat{N}\right)}}$, $\beta =1/kT$, $\hat{H}$ is the (true) Hamiltonian of the system, λ the Fermi level and $\hat{N}$ the number operator. The demonstration of (50) is given in [81, 137].

In practice, the form (50) is not very useful in nuclear physics, where the true Hamiltonian of the system is not known. The mean-field approximation provides a first simplification. It consists in replacing the true Hamiltonian $\hat{H}-\lambda \hat{N}$ by a simpler, quadratic form $\hat{K}$ of particle operators,

$$\begin{eqnarray}&&{{\hat{D}}_{\text{MF}}}=\frac{1}{{{Z}_{\text{MF}}}}{{\text{e}}^{-\beta \hat{K}}},~~~~{{Z}_{\text{MF}}}=\text{Tr}\,{{\text{e}}^{-\beta \hat{K}}}.\end{eqnarray} \tag{ 51 }$$

Given a generic basis $|i\rangle$ of the single particle space, with c_i and $c_{i}^{\dagger}$ the corresponding single particle operators, the operator $\hat{K}$ can be written formally

$$\begin{eqnarray}&&\hat{K}=\frac{1}{2}\underset{ij}{\sum}\,K_{ij}^{11}c_{i}^{\dagger}{{c}_{j}}+\frac{1}{2}\underset{ij}{\sum}\,K_{ij}^{22}{{c}_{i}}c_{j}^{\dagger}+\frac{1}{2}\underset{ij}{\sum}\,K_{ij}^{20}c_{i}^{\dagger}c_{j}^{\dagger}\\ &&+\frac{1}{2}\underset{ij}{\sum}\,K_{ij}^{02}{{c}_{i}}{{c}_{j}}.\end{eqnarray} \tag{ 52 }$$

At this point, the matrices K¹¹, K²², K²⁰ and K⁰² in (52) are still unknown.

The next step is to take advantage of the Wick theorem for ensemble averages [138]. This theorem establishes that when the density operator has the form (51), statistical traces defining the expectation value of operators can be expressed uniquely in terms of the generalized density $\mathcal{R}$, which is defined by

$$\begin{eqnarray}\mathcal{R}=\left(\begin{array}{c} \langle c_{j}^{\dagger}{{c}_{i}}\rangle & \langle {{c}_{j}}{{c}_{i}}\rangle \\ \langle c_{j}^{\dagger}c_{i}^{\dagger}\rangle & \langle {{c}_{j}}c_{i}^{\dagger}\rangle, \end{array}\right),\end{eqnarray} \tag{ 53 }$$

where $\langle \centerdot \rangle$ refers to the statistical average: for instance, $\langle c_{j}^{\dagger}{{c}_{i}}\rangle =\text{Tr}\left(\hat{D}c_{j}^{\dagger}{{c}_{i}}\right)={{\rho}_{ij}}$ is the one-body density matrix, $\langle {{c}_{j}}{{c}_{i}}\rangle =\text{Tr}\left(\hat{D}{{c}_{j}}{{c}_{i}}\right)={{\kappa}_{ij}}$ is the pairing tensor, etc. The consequence of the Wick theorem is that once the generalized density $\mathcal{R}$ of (53) is known, the expectation value $\langle \hat{F}\rangle$ of an arbitrary operator $\hat{F}$ can be computed by

$$\begin{eqnarray}&&\langle \hat{F}\rangle =\text{Tr}\left({{\hat{D}}_{\text{MF}}}\hat{F}\right)\to \langle \hat{F}\rangle =\frac{1}{2}\text{tr}\left(\mathcal{R}F\right).\end{eqnarray} \tag{ 54 }$$

In other words, the computation of the statistical trace over a many-body basis of the Fock space can be substituted by the much simpler operation of taking the trace within the quasiparticle space. In effect, this implies that the Wick theorem transfers the information content about the system from the density operator ${{\hat{D}}_{\text{MF}}}$ into the generalized density $\mathcal{R}$.

The final step is thus to determine the generalized density. This is where we take advantage of the fact that the matrix elements of the operator $\hat{K}$ in (52) are arbitrary and can be taken as variational parameters of the theory. We thus determine them by requiring that the grand potential be minimum with respect to variations $\delta \hat{K}$. This leads to the identification of K with the usual HFB matrix $\mathcal{H}$, $K=\mathcal{H}$ where $\mathcal{H}$ has the generic form (22), and to the relation

$$\begin{eqnarray}&&\mathcal{R}=\frac{1}{1+\exp \left(\beta \mathcal{H}\right)},\end{eqnarray} \tag{ 55 }$$

where $\beta =1/kT$. Equation (55) is the FT-HFB equation; see [81, 134] for the demonstration. It establishes a self-consistency condition between the HFB matrix $\mathcal{H}$ and the generalized density $\mathcal{R}$. In practice, this self-consistency condition is most easily satisfied in the basis where $\mathcal{H}$ is diagonal. In this basis (the usual quasiparticle basis of the HFB equations), one easily shows the following relation

$$\begin{eqnarray}&&\text{Tr}\left(\hat{D}\beta _{\nu}^{\dagger}{{\beta}_{\mu}}\right)=\frac{1}{1+{{\text{e}}^{\beta {{E}_{\mu}}}}}{{\delta}_{\mu \nu}}={{f}_{\mu \nu}}{{\delta}_{\mu \nu}},\end{eqnarray} \tag{ 56 }$$

with ${{E}_{\mu}}$ the quasiparticle energy, i.e. the eigenvalue μ of $\mathcal{H}$. This result allows one to extract the expression $\mathbb{R}$ of the generalized density in the quasiparticle basis, recall (15). By applying the Bogoliubov transformation, $\mathcal{R}=W\mathbb{R}{{W}^{\dagger}}$, one then finds the following expression for the density matrix and pairing tensor in the single particle basis,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\rho}_{ij}} & = & ={{\left({{V}^{\ast}}(1-f){{V}^{T}}\right)}_{ij}}+{{\left(Uf{{U}^{\dagger}}\right)}_{ij}}, \end{array}\end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\kappa}_{ij}} & = & ={{\left({{V}^{\ast}}(1-f){{U}^{T}}\right)}_{ij}}+{{\left(Uf{{V}^{\dagger}}\right)}_{ij}}, \end{array}\end{eqnarray} \tag{ 58 }$$

where U and V are the matrices of the Bogoliubov transformation.

The finite-temperature extension of the HFB theory poses two difficulties. While in the HFB theory at zero temperature the one-body density matrix and two-body pairing tensors are always localized for systems with negative Fermi energy [89, 110], this is not the case at finite temperature. In particular quasiparticles with ${{E}_{\mu}}>-\lambda$ bring a non-localized contribution to the one-body density matrix through the $\left(Uf{{U}^{\dagger}}\right)$ term. This effect is discussed in [126, 129, 139, 140]. In addition, in the statistical description of the system by a grand canonical ensemble, only the average value of the energy and of the particle number (or any other constrained observable) are fixed. Statistical fluctuations are also present [141]. They increase with temperature and decrease with system size [137]. The FT-HFB theory thus gives only the most probable solution within the grand-canonical ensemble, the one that corresponds to the lowest free energy. Mean values and deviations around the mean values of any observable $\hat{\mathcal{O}}$ should in principle be taken into account. In the classical limit, they are given by

$$\begin{eqnarray}&&\bar{\mathcal{O}}=\frac{{\int}^{}{{\text{d}}^{N}}\boldsymbol{q}\;\mathcal{O}\left(\boldsymbol{q}\right){{\text{e}}^{-\beta F\left(T,\boldsymbol{q}\right)}}}{{\int}^{}{{\text{d}}^{N}}\boldsymbol{q}\;{{\text{e}}^{-\beta F\left(T,\boldsymbol{q}\right)}}}.\end{eqnarray} \tag{ 59 }$$

Such integrals are computed across the whole collective space defined by the variables $\boldsymbol{q}=\left({{q}_{1}},\ldots,{{q}_{N}}\right)$ and require the knowledge of the volume element ${{\text{d}}^{N}}\boldsymbol{q}$. This was done for instance in [136]. Other possibilities involve functional integral methods as in [142].

Potential energy surfaces should in principle be corrected to account for beyond mean field correlations. In particular, all broken symmetries (translational, rotational, parity, particle number) should be restored, which yields an additional correlation energy to the system. Moreover, the variational HFB equation should, in principle, be solved after the projection on good quantum numbers has been performed (variation after projection, VAP).

The general problem of symmetry restoration is most naturally formulated in the framework of the SCMF model: given an effective Hamiltonian, one can associate a projector operator ${{\hat{P}}_{g}}$ to any broken symmetry and compute ${{E}_{g}}=\langle \Phi |\hat{H}{{\hat{P}}_{g}}| \Phi \rangle /\langle \Phi |{{\hat{P}}_{g}}| \Phi \rangle$. When $| \Phi \rangle$ is the symmetry-breaking HFB state minimizing the HFB energy, E_g corresponds to the projection after variation (PAV) result; if $| \Phi \rangle$ is determined from the equations obtained after varying E_g, the energy is the VAP result. The correlation energy is defined as the difference ${{E}_{g}}-{{E}_{0}}$, where ${{E}_{0}}=\langle \Phi |\hat{H}| \Phi \rangle$. Unfortunately, projection techniques, especially the VAP, are computationally very expensive when multiple symmetries are broken (as in fission). In addition, the thorough analysis of [143–146] showed that they are, strictly speaking, ill-defined as soon as $\hat{H}$ contains density-dependent terms. Finally, in spite of the recent work of [147–149], it is not really clear how such correlation energies can be rigorously computed in a strict EDF framework where no Hamiltonian is available.

In practice, the situation depends on the symmetry considered:

• (i)  
  Translational invariance: The kinetic energy of the centre of mass must be subtracted to account for the correlation energy due to restoration of translational invariance [61]. The correlation energy is most often computed as $-\langle \widehat{\boldsymbol{P}}_{\text{cm}}^{2}\rangle /2mA$, which is a first-order approximation of the VAP result as shown in [150]. Here, $\widehat{\boldsymbol{P}}={\sum}_{i}{{\widehat{\boldsymbol{p}}}_{i}}$ is the total linear momentum of the system and the expectation value is taken on the quasiparticle vacuum. In heavy nuclei, the study of [13] showed that its value can vary by about 1 MeV as a function of deformation.
• (ii)  
  Rotational invariance: For the reasons mentioned above, the correlation energy induced by angular momentum projection is also often taken into account by approximate formulas. For large deformations of the intrinsic reference state, the rotational correction energy can be approximated by ${{\epsilon}_{R}}=-{{\langle \boldsymbol{J}\rangle}^{2}}/\left(2\mathcal{J}\right)$, where $\mathcal{J}$ is the nuclear moment of inertia (which depends on the deformation) [61]. In this expression, the analyses of [150–152] showed that one should in principle use the Peierls–Yoccoz moment of inertia of [153] for the denominator although many authors use the Thouless–Valatin one. Typical values of the rotational correction range from zero MeV for spherical nuclei to 7–8 MeV for strongly quadrupole deformed configurations, as illustrated in the bottom panel of figure 6.
• (iii)  
  Reflection symmetry: Asymmetric fission is explained by invoking potential energy surfaces where the nuclear shape breaks reflection symmetry. Parity projection is, therefore, required to restore left–right symmetry. Due to the discrete nature of the symmetry (only two states involved) the correlation energy is negligible for the typically large values of the octupole moment in fission as shown in [154].
• (iv)  
  Particle number: By definition, the quasiparticle vacuum does not conserve particle number and this symmetry should also be restored. There have been very few studies of the impact of this kind of correlation energy on the PES, and the few available results, for instance of [155], are only based on the Lipkin–Nogami approximate particle number restoration scheme. Based on results obtained in other applications, this correlation energy could modify the values of the typical quantities characterizing the PES (barrier height, fission isomer excitation energy, etc) by at most a couple of MeV. The impact on the collective inertia, however, could be significant. In addition to this, the particle number breaking intrinsic states are averages of wave functions with different numbers of protons and neutrons and therefore symmetry restoration is critical to recover nuclear properties strongly dependent on particle number.

Zoom In Zoom Out Reset image size

The issue of how to describe those correlation energies in the transition from the one-fragment regime to the two-fragments one is also of paramount importance in order to determine the kinetic energy distribution of the fragments. The works of [156–159] are among the few attempts to describe such transitions in the HFB framework.

In nuclear DFT, the traditional parametrization of the nuclear shape is based on mass multipole moments ${{q}_{\lambda \mu}}$. These quantities are computed in the intrinsic frame of reference of the nucleus as expectation values of the operators

$$\begin{eqnarray}&&{{\hat{Q}}_{\lambda \mu}}={{C}_{\lambda \mu}}{\int}^{}{{\text{d}}^{3}}\boldsymbol{r}\;{{r}^{\lambda}}{{Y}_{\lambda \mu}}(\theta,\varphi ),\end{eqnarray} \tag{ 60 }$$

where ${{Y}_{\lambda \mu}}(\theta,\varphi )$ are the usual spherical harmonics, and ${{C}_{\lambda \mu}}$ are arbitrary coefficients introduced for convenience; see for example [160] for one particular choice. Since these operators are spin-independent, one-body operators, their expectation value is simply

$$\begin{eqnarray}&&{{q}_{\lambda \mu}}=\text{tr}\left({{\hat{Q}}_{\lambda \mu}}\rho \right)={{C}_{\lambda \mu}}{\int}^{}{{\text{d}}^{3}}\boldsymbol{r}\;\rho \left(\boldsymbol{r}\right){{r}^{\lambda}}{{Y}_{\lambda \mu}}(\theta,\varphi ).\end{eqnarray} \tag{ 61 }$$

The mass multipole moments ${{q}_{\lambda \mu}}$ are the analogues to the deformation parameters ${{\alpha}_{\lambda \mu}}$ of the nuclear surface parametrization of (2) used in the MM approach. In fact, since multipole moments scale like the mass of the nucleus, it is sometimes advantageous to convert them to the A-independent deformation parameters ${{\alpha}_{\lambda \mu}}$. The most commonly used technique to do so, which is recalled in [61], is to insert a constant density ${{\rho}_{0}}$ in (61), compute the volume integral for a surface defined by (2) and expand the result to first order in ${{\alpha}_{\lambda \mu}}$. The result is thus only valid for small deformations. At large deformations, other methods can be used—see, e.g. [161] and references therein for a discussion. Note that there is an important difference between the explicit shape parametrization of the MM models and the choice of multipole moments as collective variables: in the nuclear DFT the multipole moments that are not constrained take a possibly non-zero value so that the energy is minimal, whereas in the MM approach they are zero.

In the context of fission, by far the most important collective variable is the axial quadrupole moment q₂₀, which represents the elongation of the nucleus. The degree of triaxiality of nuclear shapes is captured by either q₂₂ or the ratio $\text{tan}\gamma \propto {{q}_{22}}/{{q}_{20}}$ (the exact relation depends on the convention chosen for the normalization ${{C}_{\lambda \mu}}$ of the multipole moments). Several studies, e.g. in [15, 87, 162, 163], have confirmed that triaxiality is particularly important to lower the first fission barriers of actinides. The mass asymmetry of fission fragments, especially in actinides, can be well reproduced by introducing non-zero values of the axial mass octupole moment q₃₀. The hexadecapole moment q₄₀ has been mostly used as a collective variable in the study of the neutron-induced fission of actinides: starting from the ground-state, the two-dimensional calculations of PES in the $\left({{q}_{20}},{{q}_{40}}\right)$ plane reported in [120, 123, 162, 164] showed the existence of two valleys, the fission (high q₄₀-values) and fusion (lower q₄₀ values) valleys; see also figure 28 p 49. The figure 7 gives a visual representation of the impact of each of these multipole moments on the nuclear shape.

Zoom In Zoom Out Reset image size

Let us insist that in the HFB theory the expectation value of any multipole moment can take non-zero values if symmetries allow it: for example, even when q₂₀ is the only collective variable (=the only constraint on the HFB solution), all other multipole moments will vary along the q₂₀ path in such a way as to minimize locally the total energy. By taking advantage of the non-linear properties of the HFB equations and switching on/off constraints along the iterative process, it is therefore often possible to compute an N-dimensional PES that is guaranteed to be a local minimum in the full variational space. In practice, however, one often encounters situations where this local minimum is the 'wrong' one. This point is illustrated in figure 8 for a toy-model two-dimensional collective space. In panel (a), calculations are initialized with the solution at point A and each value of q₂₀ is obtained by starting the calculation with the solution at a lower q₂₀ value. The resulting path follows the left valley in the (${{q}_{20}},{{Q}_{xx}}$) space and rejoins the right valley only when the barrier between the two vanishes; conversely, in panel (b) calculations are initialized at point B and the step-by-step process follows the right valley. The resulting trajectories are markedly different and show a discontinuity in energy. Panel (c) illustrates a possible continuous trajectory. These points were discussed extensively in [165].

Zoom In Zoom Out Reset image size

In addition to the set $\left({{q}_{20}},{{q}_{22}},{{q}_{30}},{{q}_{40}}\right)$ of multipole moments, a constraint on the size of the neck between the pre-fragments has often been used, in particular at very large elongations; see for instance [129, 162, 166–169]. The standard form of this operator is ${{\hat{Q}}_{N}}=\exp \left(-{{\left(\boldsymbol{r}-\boldsymbol{r}_{{\text{neck}}}\right)}^{2}}/{{a}^{2}}\right)$, with a an arbitrary range and $\boldsymbol{r}_{{\text{neck}}}$ the position of the point with the lowest density between the two fragments. The figure 9 shows the effect of decreasing the value of q_N near the scission point of ²⁴⁰Pu. This constraint is often used as a way to continuously approach the final configuration of two separated fragments.

Zoom In Zoom Out Reset image size

Although multipole moments are widely used to characterize nuclear shapes in low-energy nuclear structure in general, and nuclear fission in particular, they are in fact not particularly well adapted to capture some of the most relevant features of fission. In particular, the charge and mass of the fission fragments computed at scission from PES generated with multipole moments are most of the time non-integers. Although particle number projection techniques were used in [170], they were not applied on a large scale, as e.g. in the determination of fission product yields. In addition, Younes and Gogny noticed in their two-dimensional study of fission mass distributions for ²³⁹Pu(n,f) and ²³⁵U(n,f) presented in [171] that several fragmentations $\left({{Z}_{L}},{{N}_{L}}\right),\left({{Z}_{R}},{{N}_{R}}\right)$ were missing along the scission line. This is a consequence of computing a potential energy surface in a restricted collective space: discontinuities of the surface, especially at scission, become critical, and the choice $\left({{q}_{20}},{{q}_{30}}\right)$ of collective variables does not guarantee that the proper fragmentations will be recovered. For these reasons, the authors suggested adapting the practice of the MM approach by using as collective variables the distance between the two pre-fragments d and the mass asymmetry between the fragments $\xi =\left({{A}_{R}}-{{A}_{L}}\right)/A$. These quantities can be computed by introducing the spatial operators

$$\begin{eqnarray}\begin{array}{*{35}{l}} \hat{d} & = & \frac{1}{{{A}_{R}}}zH\left(z-{{z}_{\text{neck}}}\right)-\frac{1}{{{A}_{L}}}z\left[1-H\left(z-{{z}_{\text{neck}}}\right)\right], \end{array}\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \hat{\xi} & = & 2H\left(z-{{z}_{\text{neck}}}\right)-1, \end{array}\end{eqnarray} \tag{ 63 }$$

where H(x) is the Heaviside step function. As a result of this choice, they obtained a much more continuous potential energy surface. The figure 10 shows the two PES, in the $\left({{q}_{20}},{{q}_{30}}\right)$ and in the $(D,\xi )$ variables, side-by-side. Scission configurations are much better mapped out in the $(D,\xi )$ parametrization.

Zoom In Zoom Out Reset image size

In nuclear DFT, the amount of pairing correlations in the wave function is in principle automatically determined by the minimum energy principle. Altering the amount of pairing correlations in the wave function increases the energy by an amount that depends on the specific properties of the state considered (essentially, the level density) but is typically in the range of 1–2 MeV. Therefore, artificially modifying the amount of pairing correlations can modify the shape of the PES. Calculations in [162, 172] have shown that increasing pairing correlations decreases the fission barrier and leads to scission occurring at lower elongations. Most important is the strong dependence of the collective inertia on the pairing gap already pointed out in [173, 174]. The inertia is inversely proportional to the square of the pairing gap parameter and [175–178] showed that larger pairing gaps imply a smaller inertia and therefore shorter fission half-lives.

There are several possibilities to define collective variables that characterize the amount of pairing correlations in the system:

• (i)  
  The mean value of the number of particles fluctuation operator $\langle \Delta {{\hat{N}}^{2}}\rangle$ has been used in the recent calculations of [179]. The main drawback is that this is a two-body operator and therefore the computation of quantities of interest are more involved. Also its two-body character prevents the use of standard formulas used to compute the inertia.
• (ii)  
  The average value of the pairing gap $\langle \Delta \rangle$ and the gauge angle ϕ associated with the particle number operator have been mostly studied within semi-microscopic approaches based on a phenomenological mean-field complemented by a BCS description of pairing, see, e.g. [180, 181]. The formalism could in principle be extended to a full HFB framework with realistic pairing forces by defining the average value of the pairing gap as the mean value of the Cooper pair creation operator ${{P}^{\dagger}}={\sum}_{k}c_{k}^{\dagger}c_{{\bar{k}}}^{\dagger}$ times a convenient strength parameter G, i.e. $\langle \Delta \rangle =G\langle {{P}^{\dagger}}\rangle$. In this way, the HFB theory with one-body constraints and all the subsequent developments discussed below for one-body operators can be straightforwardly implemented. However, including the gauge angle ϕ may cause problems when using density-dependent forces because of the singularities analysed in [143–146].
• (iii)  
  Another possibility is to use as starting point the Lipkin–Nogami (LN) method that adds to the HFB matrix a term $-{{\lambda}_{2}} \Delta {{\hat{N}}^{2}}$ with ${{\lambda}_{2}}$ in principle determined by the LN equation. Using instead ${{\lambda}_{2}}$ as a free parameter allows to change the mean value of $\Delta {{\hat{N}}^{2}}$ at will and therefore the strength of pairing correlations. This is the choice used in [178].

Historically, the concept of scission takes its origin in the liquid drop (LD) picture of the nucleus and reflects the fact that for very large deformations, the LD potential energy can be a multi-valued function of the deformation parameters, with at least one of the solutions corresponding to two separate fragments as exemplified in [9, 41, 60]. This is illustrated, for example, in the figure 4 p 7 for the parametrization of the LD in terms of Cassini ovals. In the LD approach, these multi-valued regions originate from the finite number of collective variables (=deformations) and/or the lack of bijectivity between a set of parameters and a given geometrical shape.

In DFT, such multi-modal potential energy surfaces are encountered when working in finite collective spaces as recalled in section 2.3.1. For example, studies of hot fission of actinide nuclei published in [120, 123, 162] showed that the least energy fission pathway from ground-state to scission follows the fission valley, which lies higher in energy than the fusion valley. For given values of the axial quadrupole and octupole moments, these two valleys differ by the value of the hexadecapole moment q₄₀, and nuclear shapes in the fusion valley correspond to two well-separated fission fragments. Similar multi-modal potential energy surfaces have been observed in superheavy nuclei, for example in [186, 187]. Discontinuities in the energy (or any relevant collective variable) for large elongations of the fissionning nucleus are usually the first tell-tale signal of the transition to the scission region.

However, these discontinuities are also the manifestation of the truncated nature of the collective space: when additional collective variables such as, e.g. the size of the neck q_N, are used, they can disappear. The transition from a compact shape to very loosely joined fragments can thus be continuous. This point is illustrated in the top panel of figure 12 adapted from [162]. The graph shows the energy of ²⁴⁰Pu as a function of q_N at the point q₂₀  =  345 b. This point is located just before the discontinuity in energy along the most likely fission path of figure 29 p 49. When the q_N collective variable is used as a constraint, the system can go continuously to two separated fragments, as illustrated by the density contours of the bottom panel of the figure.

Zoom In Zoom Out Reset image size

In such cases, one needs a specific criterion to define when exactly one may consider the two fragments as fully separated. To this purpose, one can use the value of the density between the two fragments as in [188]: the fragments are deemed separated if $\parallel \rho {{\parallel}_{\text{max}}}\leqslant {{\rho}_{\text{sc}.}}$. Alternatively, one may use the expectation value of the neck operator discussed in section 2.3.1. This quantity gives a measure of the number of particles in a slice of width a centred on the neck position. The decision of considering the fragments as separated would be based on the condition $\langle {{\hat{q}}_{N}}\rangle \leqslant {{q}_{\text{sc}.}}$. Whatever the criterion retained, however, there remains a part of arbitrariness in the definition of scission: how to choose the value of ${{q}_{\text{sc}.}}$ (or ${{\rho}_{\text{sc}.}}$)?

Recently, there were attempts in [129, 162, 189] to use topological methods to characterize scission in a less arbitrary manner. The main idea is to map the density fields ${{\rho}_{\text{n}}}\left(\boldsymbol{r}\right)$ and ${{\rho}_{\text{p}}}\left(\boldsymbol{r}\right)$ into an abstract contour net, and infer properties of those fields from the connectivity properties of these nets. Mathematically rigorous, the method can identify from the PES the regions where the variations of densities within the pre-fragments are commensurate with those in the compound nucleus. As a result of this work, the authors proposed to redefine scission as a region characterized by a range of values rather than by a fixed value of either $\langle {{\hat{q}}_{N}}\rangle$ or $\parallel \rho {{\parallel}_{\text{max}}}$.

By design, geometrical definitions of scission only reflect static properties, and do not take into account the fact that the split is caused by a time-dependent competition between the repulsive Coulomb and the attractive nuclear force. As a result of this competition, however, scission may occur even when the two pre-fragments are separated by a sizeable neck. To mock up the dynamics of scission in static MM calculations, it was thus proposed already in [190] to use as a criterion for scission the ratio between the Coulomb and the nuclear forces in the neck region. These forces can be computed by taking the derivatives of the potential energy with respect to the relevant collective variable. Such techniques have been later extended in [191] by adding an additional phenomenological neck potential. This 'dynamical' definition of scission was adapted in the DFT framework by the authors [29, 184].

Both the geometrical or dynamical definition of scission are, however, semi-classical, in the sense that they ignore quantum mechanical effects in the neck region. This limitation was highlighted by Younes and Gogny in a couple of papers [184, 192]. In the framework of DFT, the degrees of freedom of the fission fragments are the one-body density matrix and the pairing tensor, which are themselves obtained from the quasiparticle wave functions. Near scission, one can introduce a localization indicator (simply related to the spatial occupation of quasiparticle wave functions) to partition the whole set of quasiparticles into two subsets belonging to either one of the pre-fragments. The total one-body density matrix of the compound nucleus is then decomposed

$$\begin{eqnarray}\begin{array}{*{35}{l}} \rho \left(\boldsymbol{r}\sigma,\boldsymbol{r}^{{\prime}}{{\sigma}^{\prime}}\right) & =\underset{ij}{\sum}\,\underset{\mu \in \{1\}}{\sum}\,V_{i\mu}^{\ast}{{V}_{j\mu}}{{\phi}_{i}}\left(\boldsymbol{r}\sigma \right)\phi _{j}^{\ast}\left(\boldsymbol{r}^{{\prime}}{{\sigma}^{\prime}}\right)\\ +\underset{ij}{\sum}\,\underset{\mu \in \{2\}}{\sum}\,V_{i\mu}^{\ast}{{V}_{j\mu}}{{\phi}_{i}}\left(\boldsymbol{r}\sigma \right)\phi _{j}^{\ast}\left(\boldsymbol{r}^{{\prime}}{{\sigma}^{\prime}}\right), \end{array}\end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & ={{\rho}_{1}}\left(\boldsymbol{r}\sigma,\boldsymbol{r}^{{\prime}}{{\sigma}^{\prime}}\right)+{{\rho}_{2}}\left(\boldsymbol{r}\sigma,\boldsymbol{r}^{{\prime}}{{\sigma}^{\prime}}\right). \end{array}\end{eqnarray} \tag{ 65 }$$

Similar expressions hold for the pairing tensor $\kappa \left(\boldsymbol{r}\sigma,\boldsymbol{r}^{{\prime}}{{\sigma}^{\prime}}\right)$. With these definitions, it becomes possible to analyse fission fragment properties, including their intrinsic energy and their interaction energy within the DFT framework; see also [162] for details. In [184], Younes and Gogny made the crucial observation that the fission fragment density distributions thus extracted have large tails that extend into the other fragment and reflect the quantum entanglement between the two fragments as shown in figure 13. Because of these large tails, there was a substantial nuclear interaction energy between the fragments even when the size of the neck was very small. Similarly, the Coulomb interaction energy was much too high compared to its experimental value.

Zoom In Zoom Out Reset image size

Most importantly, because the HFB solutions are invariant under a unitary transformation of the quasiparticle operators, it is possible to choose a representation in which this degree of entanglement is minimal, as discussed in [162, 184, 192] leading to a 'quantum localization' of the fission fragments. This freedom in choosing the representation of the HFB solutions is the analogue to localization techniques used in quantum chemistry and is discussed in some details in [162]. The implementation of the quantum localization procedure yields much more realistic estimates of fission fragment properties at scission, even for neck size up to ${{q}_{N}}\approx 0.3$. Figure 13 illustrates the impact of choosing a representation that better reproduces the asymptotic conditions of two separated fragments: by reducing the tails by about an order of magnitude, the interaction energy is reduced by a factor 3 and is close to 0 for thin necks.

Recently, the localization method was extended at finite temperature in [129]. Figure 14 illustrates the impact of minimizing the tails at different temperatures for the same case of ²⁴⁰Pu most likely fission, only with the SkM* Skyrme potential instead of the Gogny D1S. Although the practical implementation differs from the zero-temperature case because the generalized density matrix is not diagonal any longer after a rotation of the quasiparticles, it is still possible to reduce the tails of the fragment densities without changing the global properties of the compound nucleus. As the temperature increases, however, this procedure becomes more and more difficult, especially for well-entangled fragments, because of the coupling to the continuum; see section 2.2.4.

Zoom In Zoom Out Reset image size

The generator coordinate method (GCM) is a general quantum many-body technique designed to encapsulate collective correlations in the wave function. It is based on the expansion of the unknown many-body wave function of the system on a basis of known many-body states. The variational principle is used to determine the set of expansion coefficients. The technique is closely related to the configuration interaction (CI) method popular in quantum chemistry and known in nuclear physics as the 'shell model'. Usually, basis states are continuous functions of a finite set of coordinates (such as deformation parameters, Euler rotation angles, etc) whereas in the CI method, they can be unrelated to each other (for instance a set of two quasiparticle excitations). In this section, we focus on the GCM as a tool to extract a collective inertia tensor. The time-dependent extension of the GCM also provides a powerful tool to extract fission fragment distributions, and will be presented separately in section 3.3.3.

In the GCM, the general ansatz for the wave function is

$$\begin{eqnarray}&&| \Psi \rangle ={\int}^{}\text{d}\boldsymbol{q}\;f\left(\boldsymbol{q}\right)| \Phi \left(\boldsymbol{q}\right)\rangle,\end{eqnarray} \tag{ 69 }$$

where $| \Phi \left(\boldsymbol{q}\right)\rangle$ represents a set of known wave functions depending on a general label $\boldsymbol{q}$ that can include not only continuous but also discrete variables. Also the integral has to be taken in a broad sense as representing either sums over discrete values of $\boldsymbol{q}$, genuine integrals or an admixture of the two. In fission studies, the $| \Phi \left(\boldsymbol{q}\right)\rangle$ are usually quasiparticle vacua obtained by solving the HFB equation with constraints on a set of n operators ${{\hat{Q}}_{\alpha}},\alpha =1,\ldots,n$. As before, the boldface symbol $\boldsymbol{q}$ represents the set of collective variables such that $\langle \Phi \left(\boldsymbol{q}\right)|{{\hat{Q}}_{\alpha}}| \Phi \left(\boldsymbol{q}\right)\rangle ={{q}_{\alpha}}$. Applying the variational principle to the energy with the amplitudes $f\left(\boldsymbol{q}\right)$ as variational parameters leads to the Hill–Wheeler equations,

$$\begin{eqnarray}&&{\int}^{}\text{d}\boldsymbol{q}^{{\prime}}\;h\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)n\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)f\left(\boldsymbol{q}^{{\prime}}\right)=E{\int}^{}\text{d}\boldsymbol{q}^{{\prime}}\;n\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)f\left(\boldsymbol{q}^{{\prime}}\right).\end{eqnarray} \tag{ 70 }$$

It is an integral equation with the norm kernels,

$$\begin{eqnarray}&&n\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)=\langle \Phi \left(\boldsymbol{q}\right)| \Phi \left(\boldsymbol{q}^{{\prime}}\right)\rangle,\end{eqnarray} \tag{ 71 }$$

and energy kernels,

$$\begin{eqnarray}&&h\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)=\frac{\langle \Phi \left(\boldsymbol{q}\right)|\hat{H}| \Phi \left(\boldsymbol{q}^{{\prime}}\right)\rangle}{\langle \Phi \left(\boldsymbol{q}\right)| \Phi \left(\boldsymbol{q}^{{\prime}}\right)\rangle}.\end{eqnarray} \tag{ 72 }$$

Since the set of wave functions $| \Phi \left(\boldsymbol{q}\right)\rangle$ is in general not orthogonal, the $f\left(\boldsymbol{q}\right)$ amplitudes cannot be interpreted as probability amplitudes, and $|\,f\left(\boldsymbol{q}\right){{|}^{2}}$ is not a probability density. In order to obtain probability amplitudes, the set $| \Phi \left(\boldsymbol{q}\right)\rangle$ has to be orthogonalized using standard techniques of linear algebra to obtain what are called 'natural states' $|\rm{\tilde \Phi }\left(\boldsymbol{q}\right)\rangle$. These states are defined by folding $| \Phi \left(\boldsymbol{q}\right)\rangle$ with the inverse of the square root of the norm kernel,

$$\begin{eqnarray}&&|\rm{\tilde \Phi }\left(\boldsymbol{q}\right)\rangle ={\int}^{}\text{d}\boldsymbol{q}^{{\prime}}{{\left[n\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)\right]}^{-1/2}}| \Phi \left(\boldsymbol{q}^{{\prime}}\right)\rangle.\end{eqnarray} \tag{ 73 }$$

The square root of the norm kernel is defined as

$$\begin{eqnarray}&&n\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)={\int}^{}\text{d}\boldsymbol{q}^{{\prime \prime}}{{\left[n\left(\boldsymbol{q},\boldsymbol{q}^{{\prime \prime}}\right)\right]}^{1/2}}{{\left[n\left(\boldsymbol{q}^{{\prime \prime}},\boldsymbol{q}^{{\prime}}\right)\right]}^{1/2}}\end{eqnarray} \tag{ 74 }$$

which corresponds to a Cholesky decomposition of the positive-definite norm kernel.

Collective Schrödinger equation.

The connection of the GCM with the theory of fission comes from the reduction of the Hill–Wheeler equation (70) to a collective Schrödinger-like equation (CSE) of the collective coordinates, which naturally leads to the definition of a collective inertia. Following [12, 198–200], the CSE is derived after assuming that the norm overlap (71) is a sharply peaked function of the coordinate difference $\boldsymbol{s}=\boldsymbol{q}-\boldsymbol{q}^{{\prime}}$ and smoothly depends upon the average $\overline{\boldsymbol{q}}=\frac{1}{2}\left(\boldsymbol{q}+\boldsymbol{q}^{{\prime}}\right)$ in such a way that the norm kernel is well approximated by a Gaussian

$$\begin{eqnarray}&&n\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)=\exp \left(-\frac{1}{2}\boldsymbol{s}{\mathsf {\Gamma }}\left(\overline{\boldsymbol{q}}\right)\boldsymbol{s}\right).\end{eqnarray} \tag{ 75 }$$

In this expression the width ${\mathsf {\Gamma }}$ is to be understood as a rank 2 tensor with components ${{ \Gamma }_{\alpha \beta}}$ and the exponent is thus given explicitly by $-\frac{1}{2}{\sum}_{\alpha \beta}{{ \Gamma }_{\alpha \beta}}{{s}_{\alpha}}{{s}_{\beta}}$. If we assume that the components ${{ \Gamma }_{\alpha \beta}}$ of the width are slowly varying functions of the coordinate $\overline{\boldsymbol{q}}$, they can be related to the norm kernel by ${{ \Gamma }_{\alpha \beta}}=\frac{\partial}{\partial {{q}_{\alpha}}}\frac{\partial}{\partial q_{\beta}^{\prime}}n\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)$. This expression can also be computed using standard linear response techniques from the alternative definition

$$\begin{eqnarray}&&{{ \Gamma }_{\alpha \beta}}=\langle \Phi \left(\boldsymbol{q}\right)|\frac{{\overset{\scriptscriptstyle\leftarrow}{\partial}}}{\partial {{q}_{\alpha}}}\frac{{\vec{\partial}}}{\partial {{q}_{\beta}}}| \Phi \left(\boldsymbol{q}\right)\rangle.\end{eqnarray} \tag{ 76 }$$

which involves the 'momentum operator' $\frac{{\vec{\partial}}}{\partial {{q}_{\beta}}}$ see below. Another, simpler way, is to evaluate the overlap of the two HFB wave functions for near $\boldsymbol{q}$ and $\boldsymbol{q}^{{\prime}}$ values using the Onishi formula [61] and make a local fit to a Gaussian. If the Gaussian overlap approximation (GOA) is valid, it does not make sense to accurately compute the Hamiltonian overlap for values of $\boldsymbol{s}$ greater than the inverse of the square root of the width ${{{\mathsf {\Gamma }}}^{\,1/2}}\left(\overline{\boldsymbol{q}}\right)$³. Therefore, a reasonable approximation is to expand $h\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)$ around $\boldsymbol{s}=0$ (or $\boldsymbol{q}=\boldsymbol{q}^{{\prime}}=\overline{\boldsymbol{q}}$) and keep quadratic terms only,

$$\begin{eqnarray}\begin{array}{*{35}{l}} h\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right) & =h\left(\overline{\boldsymbol{q}},\overline{\boldsymbol{q}}\right) \\ {} & +\,{{h}_{\boldsymbol{q}}}\left(\boldsymbol{q}-\overline{\boldsymbol{q}}\right)+{{h}_{\boldsymbol{q}^{{\prime}}}}\left(\boldsymbol{q}^{{\prime}}-\overline{\boldsymbol{q}}\right) \\ {} & +\,\frac{1}{2}\left[\left(\boldsymbol{q}-\overline{\boldsymbol{q}}\right){{{\mathsf {H}}}_{\boldsymbol{q}\boldsymbol{q}}}\left(\boldsymbol{q}-\overline{\boldsymbol{q}}\right)+\left(\boldsymbol{q}^{{\prime}}-\overline{\boldsymbol{q}}\right){{{\mathsf {H}}}_{\boldsymbol{q}^{{\prime}}\boldsymbol{q}^{{\prime}}}}\left(\boldsymbol{q}^{{\prime}}-\overline{\boldsymbol{q}}\right) \right. \\ {} & \left. \;\;\;\;\;\;\;+\,2\left(\boldsymbol{q}-\overline{\boldsymbol{q}}\right){{{\mathsf {H}}}_{\boldsymbol{q}\boldsymbol{q}^{{\prime}}}}\left(\boldsymbol{q}^{{\prime}}-\overline{\boldsymbol{q}}\right)\right]+\cdots \end{array}\end{eqnarray} \tag{ 77 }$$

In this expression ${{h}_{\boldsymbol{q}}}$ denotes the set of partial derivatives with respect to $\boldsymbol{q}$ of the energy kernels,

$$\begin{eqnarray}&&{{h}_{\boldsymbol{q}}}\equiv \left({{h}_{{{q}_{1}}}},\ldots,{{h}_{qN}}\right),~~~{{h}_{{{q}_{\alpha}}}}=\frac{\partial h\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)}{\partial {{q}_{\alpha}}}{{|}_{\boldsymbol{q}=\boldsymbol{q}^{{\prime}}=\overline{\boldsymbol{q}}}}.\end{eqnarray} \tag{ 78 }$$

This quantity is a vector and products like ${{h}_{\boldsymbol{q}}}\left(\boldsymbol{q}-\overline{\boldsymbol{q}}\right)$ have to be understood as scalar products ${\sum}_{\alpha}{{h}_{{{q}_{\alpha}}}}\left({{q}_{\alpha}}-{{\bar{q}}_{\alpha}}\right)$. On the other hand, ${{{\mathsf {H}}}_{\boldsymbol{q}\boldsymbol{q}^{{\prime}}}}$ denotes the set of second partial derivatives with respect to $\boldsymbol{q}$ and $\boldsymbol{q}^{{\prime}}$ and therefore it is a rank 2 tensor with components ${{H}_{{{q}_{\alpha}}q_{\beta}^{\prime}}}$

$$\begin{eqnarray}&&{{{\mathsf {H}}}_{\boldsymbol{q}\boldsymbol{q}^{{\prime}}}}\equiv {{\left({{H}_{{{q}_{\alpha}}q_{\beta}^{\prime}}}\right)}_{\alpha,\beta =1,\ldots,N}},~~~{{H}_{{{q}_{\alpha}}q_{\beta}^{\prime}}}=\frac{{{\partial}^{2}}h\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)}{\partial {{q}_{\alpha}}\partial q_{\beta}^{\prime}}{{|}_{\boldsymbol{q}=\boldsymbol{q}^{{\prime}}=\overline{\boldsymbol{q}}}},\end{eqnarray} \tag{ 79 }$$

and products like ${{{\mathsf {H}}}_{\boldsymbol{q}\boldsymbol{q}^{{\prime}}}}\left(\boldsymbol{q}-\overline{\boldsymbol{q}}\right)\left(\boldsymbol{q}^{{\prime}}-\overline{\boldsymbol{q}}\right)$ are to be understood as contractions of this tensor with the two vectors $\boldsymbol{q}-\overline{\boldsymbol{q}}$ and $\boldsymbol{q}^{{\prime}}-\overline{\boldsymbol{q}}$, namely ${\sum}_{\alpha \beta}{{H}_{{{q}_{\alpha}}q_{\beta}^{\prime}}}\left({{q}_{\alpha}}-{{\bar{q}}_{\alpha}}\right)\left(q_{\beta}^{\prime}-{{\bar{q}}_{\beta}}\right)$. Using the exponential form of the norm kernels, assuming a constant width ${\mathsf {\Gamma }}$ and the quadratic expansion (77) for the energy kernels, it is possible to reduce the Hill–Wheeler equation to the following CSE in the collective space defined by the variables $\boldsymbol{q}$,

$$\begin{eqnarray}&&\left(-\frac{{{\hslash}^{2}}}{2}\frac{\partial}{\partial \boldsymbol{q}}{{{\mathsf {B}}}_{\text{GCM}}}\left(\boldsymbol{q}\right)\frac{\partial}{\partial \boldsymbol{q}}+V\left(\boldsymbol{q}\right)-{{\epsilon}_{\text{zpe}}}\left(\boldsymbol{q}\right)\right){{g}_{\sigma}}\left(\boldsymbol{q}\right)={{\epsilon}_{\sigma}}{{g}_{\sigma}}\left(\boldsymbol{q}\right).\end{eqnarray} \tag{ 80 }$$

The derivation is given in [12, 61, 199]. The collective mass ${{{\mathsf {M}}}_{\text{GCM}}}\left(\boldsymbol{q}\right)$ is a rank 2 tensor that depends on $\boldsymbol{q}$ and is the inverse to the collective inertia ${{{\mathsf {B}}}_{\text{GCM}}}\left(\boldsymbol{q}\right)$ (curvature) of the collective Hamiltonian when expanded to second order in the variable $\boldsymbol{s}$,

$$\begin{eqnarray}&&{\mathsf {M}}_{\text{GCM}}^{-1}\left(\boldsymbol{q}\right)\equiv {{{\mathsf {B}}}_{\text{GCM}}}\left(\boldsymbol{q}\right)=\frac{1}{2}{{{\mathsf {\Gamma }}}^{-1}}\left({{{\mathsf {H}}}_{\boldsymbol{q}\boldsymbol{q}^{{\prime}}}}-{{{\mathsf {H}}}_{\boldsymbol{q}\boldsymbol{q}}}\right){{{\mathsf {\Gamma }}}^{-1}}.\end{eqnarray} \tag{ 81 }$$

The potential energy $V\left(\boldsymbol{q}\right)$ is the HFB energy and ${{\epsilon}_{\text{zpe}}}\left(\boldsymbol{q}\right)$ is a zero point energy correction,

$$\begin{eqnarray}&&{{\epsilon}_{\text{zpe}}}\left(\boldsymbol{q}\right)=\frac{1}{2}{{{\mathsf {H}}}_{\boldsymbol{q}\boldsymbol{q}^{{\prime}}}}{{{\mathsf {\Gamma }}}^{-1}}\end{eqnarray} \tag{ 82 }$$

given in terms of the contraction of the two tensors to form a scalar. The zero point energy (zpe) correction represents a quantum correction to the classical potential for the collective variables and given by the HFB energy. The zpe correction corresponds to the energy of a Gaussian wave packet [61, 201]. The interpretation of this term is similar to the rotational energy correction but with the rotation angles replaced by the collective variables. The solution of (80) provides the energies ${{\epsilon}_{\sigma}}$ and wave functions ${{g}_{\sigma}}\left(\boldsymbol{q}\right)$ of the collective modes. The wave functions ${{g}_{\sigma}}\left(\boldsymbol{q}\right)$ are related to the amplitudes $f\left(\boldsymbol{q}\right)$ of the Hill–Wheeler equation (70) through the relation

$$\begin{eqnarray}&&{{f}_{\sigma}}\left(\boldsymbol{q}\right)={\int}^{}\text{d}\boldsymbol{q}^{{\prime}}n{{\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)}^{-1/2}}{{g}_{\sigma}}\left(\boldsymbol{q}^{{\prime}}\right),\end{eqnarray} \tag{ 83 }$$

which can be used to compute, for instance, mean values of other observables where the Gaussian approximation is not justified. Apart from the simplification that the local reduction brings to the solution of this problem, an interesting physical picture emerges: the collective dynamics is driven by the behaviour of the potential energy surface (PES) given by the HFB energy as a function of $\boldsymbol{q}$ with some coordinate-dependent quantal corrections ${{\epsilon}_{\text{zpe}}}\left(\boldsymbol{q}\right)$ and inertia ${{{\mathsf {M}}}_{\text{GCM}}}\left(\boldsymbol{q}\right)$. If the $\boldsymbol{q}$ are chosen to be the collective variables driving the nucleus to fission (quadrupole moment, octupole moment, neck, etc) then the probability of tunnelling through the fission barrier is given by the integral of the exponential of the action computed with the collective potential and the collective inertia of the CSE approximation to the Hill–Wheeler equation of the GCM.

Local approximation.

The collective inertia can be evaluated from the energy kernels using numerical differentiation to obtain ${{h}_{\boldsymbol{q}}}$ and ${{{\mathsf {H}}}_{\boldsymbol{q}\boldsymbol{q}^{{\prime}}}}$. Various finite difference schemes can be used and the precision is controlled by the value of $\delta \boldsymbol{q}$ and $\delta \boldsymbol{q}^{{\prime}}$. However, it is also common practice to use the explicit expression for the 'momentum' operator associated to the variable $\boldsymbol{q}$. Following [200, 202], we write

$$\begin{eqnarray}&&\frac{\partial}{\partial {{q}_{\alpha}}}| \Phi \left(\boldsymbol{q}\right)\rangle =\frac{\text{i}}{\hslash}{{\hat{P}}_{{{q}_{\alpha}}}}| \Phi \left(\boldsymbol{q}\right)\rangle.\end{eqnarray} \tag{ 84 }$$

The action of the momentum operator ${{\hat{P}}_{q\alpha}}$ on the quasiparticle vacuum $| \Phi \left(\boldsymbol{q}\right)\rangle$ can be obtained from the Ring and Schuck theorem of [203],

$$\begin{eqnarray}&&{{\hat{P}}_{{{q}_{\alpha}}}}| \Phi \left(\boldsymbol{q}\right)\rangle =\underset{\mu <\nu}{\sum}\,\left[{{\left(P_{{{q}_{\alpha}}}^{20}\right)}_{\mu \nu}}\beta _{\mu}^{\dagger}\beta _{\nu}^{\dagger}-{{\left(P_{{{q}_{\alpha}}}^{20\ast}\right)}_{\mu \nu}}{{\beta}_{\mu}}{{\beta}_{\nu}}\right]| \Phi \left(\boldsymbol{q}\right)\rangle.\end{eqnarray} \tag{ 85 }$$

The quasiparticle matrix elements of the momentum operator ${{\left(P_{{{q}_{\alpha}}}^{20}\right)}_{\mu \nu}}$ are obtained by expanding the HFB solution at point $\boldsymbol{q}+\delta \boldsymbol{q}$ to first order in $\delta \boldsymbol{q}$. They are related to the matrix of the derivatives of the generalized density with respect to the collective variables,

$$\begin{eqnarray}\frac{\partial \mathcal{R}}{\partial {{q}_{\alpha}}}=\frac{\text{i}}{\hslash}\left(\begin{array}{c} 0 & P_{{{q}_{\alpha}}}^{20} \\ -P_{{{q}_{\alpha}}}^{20\ast} & 0 \end{array}\right),\end{eqnarray} \tag{ 86 }$$

and take the generic form

$$\begin{eqnarray}&&\left(\begin{array}{c} {{\left(P_{{{q}_{\alpha}}}^{20}\right)}_{\mu \nu}} \\ -{{\left(P_{{{q}_{\alpha}}}^{20\,\ast}\right)}_{\mu \nu}} \end{array}\right)=\underset{\beta}{\sum}\,{{\left[M_{\alpha \beta}^{(-1)}\right]}^{-1}}\underset{{{\mu}^{\prime}}<{{\nu}^{\prime}}}{\sum}\,{{\left({{\mathcal{M}}^{-1}}\right)}_{\mu \nu {{\mu}^{\prime}}{{\nu}^{\prime}}}}\left(\begin{array}{c} {{\left(q_{\beta}^{20}\right)}_{{{\mu}^{\prime}}{{\nu}^{\prime}}}} \\ {{\left(q_{\beta}^{20\,\ast}\right)}_{{{\mu}^{\prime}}{{\nu}^{\prime}}}} \end{array}\right),\end{eqnarray} \tag{ 87 }$$

where $\mathcal{M}$ is the linear response matrix in the quasiparticle basis,

$$\begin{eqnarray}\mathcal{M}=\left(\begin{array}{c} A & B \\ {{B}^{\ast}} & {{A}^{\ast}} \end{array}\right),\end{eqnarray} \tag{ 88 }$$

and A and B are given by (see, e.g. [61, 81])

$$\begin{eqnarray}&&{{A}_{\mu \nu {{\mu}^{\prime}}{{\nu}^{\prime}}}}=\langle {{\beta}_{\nu}}{{\beta}_{\mu}}\hat{H}\beta _{{{\mu}^{\prime}}}^{\dagger}\beta _{{{\nu}^{\prime}}}^{\dagger}\rangle,~~~{{B}_{\mu \nu {{\mu}^{\prime}}{{\nu}^{\prime}}}}=\langle {{\beta}_{\nu}}{{\beta}_{\mu}}{{\beta}_{{{\nu}^{\prime}}}}{{\beta}_{{{\mu}^{\prime}}}}\hat{H}\rangle.\end{eqnarray} \tag{ 89 }$$

In (87), we have introduced the moments

$$\begin{eqnarray}&&M_{\alpha \beta}^{(-n)}=\left(q_{\alpha}^{20\,\dagger},~q_{\alpha}^{20\,T}\right){{\mathcal{M}}^{-n}}\left(\begin{array}{c} q_{\beta}^{20} \\ q_{\beta}^{20\,\ast} \end{array}\right).\end{eqnarray} \tag{ 90 }$$

where ${{\mathcal{M}}^{-n}}={{\mathcal{M}}^{-1}}\times \ldots {{\mathcal{M}}^{-1}}$ n times. Using the explicit form of the momentum operator we can then express the width tensor in terms of the moments (which are also rank 2 tensors)

$$\begin{eqnarray}&&{\mathsf {\Gamma }}\,\,=\frac{1}{2}{{\left[{{{\mathsf {M}}}^{(-1)}}\left(\boldsymbol{q}\right)\right]}^{-1}}{{{\mathsf {M}}}^{(-2)}}\left(\boldsymbol{q}\right){{\left[{{{\mathsf {M}}}^{(-1)}}\left(\boldsymbol{q}\right)\right]}^{-1}}.\end{eqnarray} \tag{ 91 }$$

The linear response matrices of (89) have been defined only in the case of interactions deriving from a Hamiltonian $\hat{H}$. For density-dependent interactions and generic EDFs that cannot be expressed as mean values of a Hamiltonian, the generalization of (89) involves second derivatives of the energy with respect to the variational parameters. Typical expressions are given in [81] in the most general case. The idea is to use the short version of Thouless theorem [61] $| \Phi (Z)\rangle =\exp \left({\sum}_{\mu <\nu}{{Z}_{\mu \nu}}\beta _{\mu}^{\dagger}\beta _{\nu}^{\dagger}\right)| \Phi (0)\rangle$ to define the density and pairing tensor entering the EDF as functions of the independent variables ${{Z}_{\mu \nu}}$ and $Z_{\mu \nu}^{\ast}$ (for instance, ${{\rho}_{ji}}\left({{Z}^{\ast}},Z\right)=\langle \Phi (Z)|c_{i}^{\dagger}{{c}_{j}}| \Phi (Z)\rangle /\langle \Phi (Z)| \Phi (Z)\rangle$) and define

$$\begin{eqnarray}&&{{A}_{\mu \nu {{\mu}^{\prime}}{{\nu}^{\prime}}}}=\frac{\partial E\left({{Z}^{\ast}},Z\right)}{\partial {{Z}_{{{\mu}^{\prime}}{{\nu}^{\prime}}}}\partial Z_{\nu \mu}^{\ast}}~~~{{B}_{\mu \nu {{\mu}^{\prime}}{{\nu}^{\prime}}}}=\frac{\partial E\left({{Z}^{\ast}},Z\right)}{\partial Z_{{{\nu}^{\prime}}{{\mu}^{\prime}}}^{\ast}\partial Z_{\nu \mu}^{\ast}}.\end{eqnarray} \tag{ 92 }$$

Cranking approximation.

The evaluation (90) of the moments ${{{\mathsf {M}}}^{(-n)}}$ requires inverting the full linear response matrix $\mathcal{M}$. Computationally, this is a daunting task that is often alleviated by approximating $\mathcal{M}$ by a diagonal matrix, which simplifies the inversion problem enormously. This 'cranking approximation' corresponds to neglecting the residual quasiparticle interaction. The diagonal matrix elements are simply the two-quasiparticle excitation energies. With this simplification, (90) reduces to the more manageable form

$$\begin{eqnarray}&&M_{\alpha \beta}^{(-n)}=\underset{\mu <\nu}{\sum}\,\frac{\langle \Phi |\hat{Q}_{\alpha}^{\dagger}|\mu \nu \rangle \langle \mu \nu |{{{\hat{Q}}}_{\beta}}| \Phi \rangle}{{{\left({{E}_{\mu}}+{{E}_{\nu}}\right)}^{n}}},\end{eqnarray} \tag{ 93 }$$

where $| \Phi \rangle$ stands for the quasiparticle vacuum at point $\boldsymbol{q}$ and $|\mu \nu \rangle$ represents a two quasiparticle excitation built on top of that vacuum, i.e. $|\mu \nu \rangle =\beta _{\mu}^{\dagger}\beta _{\nu}^{\dagger}| \Phi \rangle$. The combination of introducing the local momentum operator to substitute for exact numerical differentiations with the cranking approximation is also referred to as the 'perturbative cranking' approximation. In this case, the curvature term ${{{\mathsf {H}}}_{\boldsymbol{q}\boldsymbol{q}}}$ vanishes and some algebra leads to

$$\begin{eqnarray}&&{{{\mathsf {H}}}_{\boldsymbol{q}\boldsymbol{q}^{{\prime}}}}=\frac{1}{2}{{\left[{{{\mathsf {M}}}^{(-1)}}\right]}^{-1}}.\end{eqnarray} \tag{ 94 }$$

Introducing this last result in (81) while taking into account the form (91) for the norm overlap, the collective mass tensor reduces to

$$\begin{eqnarray}&&{{{\mathsf {M}}}_{\text{GCM}}}\left(\boldsymbol{q}\right)=4{\mathsf {\Gamma }}{{{\mathsf {M}}}^{(-1)}}{\mathsf {\Gamma }}.\end{eqnarray} \tag{ 95 }$$

The zero point energy correction becomes

$$\begin{eqnarray}&&{{\epsilon}_{\text{zpe}}}\left(\boldsymbol{q}\right)=\frac{1}{2}{\mathsf {\Gamma }}{\mathsf {M}}_{\text{GCM}}^{-1}\left(\boldsymbol{q}\right).\end{eqnarray} \tag{ 96 }$$

This last expression must again be understood as the contraction of the two tensors. Both expressions (95) and (96) are then used in the evaluation of the action $S\left(\boldsymbol{q}\right)$ that enters the WKB expression of the spontaneous fission half-life (114).

Variable width.

The above discussion has been restricted for simplicity to a constant Gaussian width ${\mathsf {\Gamma }}$. In case this assumption is not strictly valid, a change of variables to a new set of coordinates $\boldsymbol{\eta }$ is required. The new variables are defined to make the width locally constant, [175]

$$\begin{eqnarray}&&\underset{\alpha}{\sum}\,\left(\text{d}\eta \right)_{\alpha}^{2}=\underset{\alpha \beta}{\sum}\,\text{d}{{\left(\boldsymbol{q}-\boldsymbol{q}^{{\prime}}\right)}_{\alpha}}{{ \Gamma }_{\alpha \beta}}\left(\boldsymbol{q}\right)\text{d}{{\left(\boldsymbol{q}-\boldsymbol{q}^{{\prime}}\right)}_{\beta}}.\end{eqnarray} \tag{ 97 }$$

This last expression is reminiscent of the field of differential geometry where the new variables $\boldsymbol{\eta }$ have a constant metric ${{\delta}_{\alpha \beta}}$ and distances in the original variables $\boldsymbol{q}$ have to be measured with the new metric ${\mathsf {\Gamma }}$. In this framework the expressions derived above are still valid if standard derivatives with respect to ${{q}_{\alpha}}$ are replaced by covariant derivatives $\mathcal{D}/\mathcal{D}{{q}_{\alpha}}$ that include in their definition the corresponding Christoffel symbols [12, 175]. For an introduction to differential geometry, the reader may refer to chapter 5 of [204]. As the metric in the variables $\boldsymbol{q}$ is now coordinate-dependent, the volume element of integrals must be modified to account for an extra $\sqrt{\det {\mathsf {\Gamma }}}$, which is the determinant of ${\mathsf {\Gamma }}$. Also the kinetic energy term has to incorporate the bells and whistles of differential geometry involving covariant derivatives. The final expression reads

$$\begin{eqnarray}&&\left(-\frac{{{\hslash}^{2}}}{2}\frac{1}{\sqrt{\det {\mathsf {\Gamma }}}}\frac{\partial}{\partial \boldsymbol{q}}\sqrt{\det {\mathsf {\Gamma }}}\;{{{\mathsf {B}}}_{\text{GCM}}}\left(\boldsymbol{q}\right)\frac{\partial}{\partial \boldsymbol{q}}+V\left(\boldsymbol{q}\right)-{{\epsilon}_{\text{zpe}}}\left(\boldsymbol{q}\right)\right){{g}_{\sigma}}\left(\boldsymbol{q}\right)={{\epsilon}_{\sigma}}{{g}_{\sigma}}\left(\boldsymbol{q}\right).\end{eqnarray} \tag{ 98 }$$

The use of covariant derivatives is still required in the evaluation of the inertia. As the calculation of the Christoffel symbols and its use in the covariant derivatives is rather cumbersome, it is often assumed that the width matrix ${\mathsf {\Gamma }}$, which is used also as the metric of the curved space, varies slowly with the $\boldsymbol{q}$ coordinates. In this case, and given that the Christoffel symbols are defined in terms of partial derivatives of the metric, they can be neglected and therefore the covariant derivative reduces to the standard partial derivative. However the term $\sqrt{\det {\mathsf {\Gamma }}}$ is kept in its original form.

For a given choice of collective variables $\left\{\boldsymbol{q}\right\}$, the ansatz (69) for the GCM wave function leads to a notorious underestimation of the collective inertia ${{{\mathsf {M}}}_{\text{GCM}}}\left(\boldsymbol{q}\right)$, even when computed with exact numerical derivatives. This is illustrated in the top panel of figure 6 p 16, which shows the collective inertia as a function of the quadrupole deformation for both the GCM and ATDHFB prescriptions (in the perturbative cranking approximation). Ring and Schuck recall, in their textbook [61], how one-dimensional GCM calculations using a shift of the centre of mass as collective variable $\boldsymbol{q}$ fail to reproduce the exact result in the exactly solvable case of translational large amplitude collective motion (where the collective mass is just the sum of the masses of all nucleons). Back in 1962, Peierls and Yoccoz had showed in [205] that adding another collective variable corresponding to the momentum of the centre of mass was sufficient to reproduce the exact mass. Qualitatively, a naive implementation of the GCM where the only collective degrees of freedom $\boldsymbol{q}$ are time-even functions (such as, e.g. HFB states under various constraints on multipole moments) is bound to fail, since it does not contain enough information on the actual motion in the collective space, which is controlled by time-odd momenta. In fact, Reinhard and Goeke showed in their review [206], that a 'dynamic GCM' (DGCM), where the set of collective variables $\boldsymbol{q}$ is expanded to include the associated momenta ${{\widehat{\boldsymbol{P}}}_{\boldsymbol{q}}}$, was necessary to provide a more realistic description of collective dynamics. The implementation of the DGCM in practical applications is, however, a lot more involved than the usual GCM.

The adiabatic time-dependent HFB (ATDHFB) approximation of the TDHFB equation thus provides an appealing alternative; see, e.g. [207–216] for a presentation of the theory. The ATDHFB is based on a small velocities expansion of the TDHFB equation,

$$\begin{eqnarray}&&\text{i}\hslash \overset{\centerdot}{{\mathcal{R}}}\,=\left[\mathcal{H},\mathcal{R}\right],\end{eqnarray} \tag{ 99 }$$

where the generalized density matrix $\mathcal{R}(t)$ is given by (14) and the HFB matrix by (22)—both are now time-dependent. As we will show below, this expansion introduces a set of collective 'coordinates', which are time-even generalized densities, and the related collective 'momenta'. Coordinates and momenta differ by their properties with respect to time-reversal symmetry. Once these quantities are defined, the energy of the system can be expressed as the sum of a kinetic part, which is a quadratic function of the momenta, and a potential part—both parts being time-dependent functions. It is then possible to assign a collective inertia associated to any point of the collective space. At this point the problem is reduced to a classical one and it is not possible to describe phenomena involving quantum tunnelling through the barrier. However, we can still resort to the semi-classical description of tunnelling based on the WKB method. Within this scheme it is thus possible to compute spontaneous fission lifetimes using quantities provided by the HFB theory and its time-dependent extension ATDHFB.

Note that all applications of ATDHFB to fission so far have been performed in the particular case where the collective path, i.e. the trajectory of the system in the collective space, is not calculated dynamically by solving the ATDHFB equations, but predefined as a set of constrained HFB calculations. This is so even though (i) there are several possible prescriptions to compute the collective path, as can be seen, e.g. in the studies of [217–219] and (ii) comparisons of dynamically-defined collective paths with constrained ones indicate that there can be significant differences in the collective inertia tensor and zero-point energy corrections; see the examples studied in [158, 215]. The main reason for this approximation is that the ATDHFB theory is mostly used as a tool to compute a realistic collective inertia for WKB-types of calculations, and the actual evolution of the system in collective space is not needed. The adiabatic self-consistent collective model formulated in [220, 221] represents an alternative formulation of the adiabatic approximation to the TDHFB equations that solve many of the formal difficulties but has not been applied to the specific case of fission yet.

The starting point of the ATDHFB method is the expansion of the full time-dependent generalized density $\mathcal{R}(t)$ according to

$$\begin{eqnarray}&&\mathcal{R}(t)={{\text{e}}^{-\text{i}\chi (t)}}{{\mathcal{R}}_{0}}(t){{\text{e}}^{\text{i}\chi (t)}}\end{eqnarray}$$

where ${{\mathcal{R}}_{0}}(t)$ is a time-dependent, time-even generalized density matrix satisfying the standard relation $\mathcal{R}_{0}^{2}={{\mathcal{R}}_{0}}$. $\chi (t)$ is also a time-even operator. The time-even density ${{\mathcal{R}}_{0}}(t)$ is considered as some kind of coordinate variable whereas $\chi (t)$ will be related to its conjugate momentum, which is assumed to be a small quantity. Expanding the density matrix in powers of $\chi (t)$

$$\begin{eqnarray}&&\mathcal{R}(t)={{\mathcal{R}}_{0}}(t)+{{\mathcal{R}}_{1}}(t)+{{\mathcal{R}}_{2}}(t)+\cdots \end{eqnarray} \tag{ 100 }$$

we obtain terms which are either time-odd, such as ${{\mathcal{R}}_{1}}=-\text{i}\left[\chi,{{\mathcal{R}}_{0}}\right]$, or time-even, such as ${{\mathcal{R}}_{0}}$ or ${{\mathcal{R}}_{2}}$. Introducing the above expansion in the TDHFB equation, we obtain a set of two equations, one for time-odd and the other for time-even quantities,

$$\begin{eqnarray}\begin{array}{*{35}{l}} \text{i}\hslash {{\overset{\centerdot}{{\mathcal{R}}}\,}_{0}} & = & \left[{{\mathcal{H}}_{0}},{{\mathcal{R}}_{1}}\right]+\left[{{\mathcal{K}}_{1}},{{\mathcal{R}}_{0}}\right], \end{array}\end{eqnarray} \tag{ 101 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \text{i}\hslash {{\overset{\centerdot}{{\mathcal{R}}}\,}_{1}} & =\left[{{\mathcal{H}}_{0}},{{\mathcal{R}}_{0}}\right]+ & \left[{{\mathcal{H}}_{0}},{{\mathcal{R}}_{2}}\right]+\left[{{\mathcal{K}}_{1}},{{\mathcal{R}}_{1}}\right]+\left[{{\mathcal{K}}_{2}},{{\mathcal{R}}_{0}}\right], \end{array}\end{eqnarray} \tag{ 102 }$$

where ${{\mathcal{H}}_{\mu}}$ and ${{\mathcal{K}}_{\mu}}$ are obtained from the expressions (22) and (29) of $\mathcal{H}$ and $\mathcal{K}$ with the densities ${{\mathcal{R}}_{\mu}}$—all being time-dependent quantities. In the next step, the TDHFB energy is expanded according to (100) in order to obtain a zero-order term, which resembles the HFB energy for the ${{\mathcal{R}}_{0}}$ density, and a second order term reminiscent of a kinetic energy and given by

$$\begin{eqnarray}&&\frac{1}{2}\text{tr}\left\{{{\overset{\centerdot}{{\mathcal{R}}}\,}_{0}}\left[{{\mathcal{R}}_{0}},{{\mathcal{R}}_{1}}\right]\right\}.\end{eqnarray} \tag{ 103 }$$

In this expression, the first term ${{\overset{\centerdot}{{\mathcal{R}}}\,}_{0}}$ plays the role of a generalized velocity, whereas the second one is a momentum-like quantity. The identification of this term with a kinetic energy is at the origin of the definition of the ATDHFB mass. In order to obtain a more explicit definition of the mass, we use (101) to express ${{\mathcal{R}}_{1}}$ in terms of ${{\overset{\centerdot}{{\mathcal{R}}}\,}_{0}}$. To this end it is customary to introduce the ATDHFB basis, which diagonalizes ${{\mathcal{R}}_{0}}(t)$ at all times t and has an analogous block structure as the traditional HFB basis. In the ATDHFB basis, the matrix of $\hat{\chi}$ is noted generically

$$\begin{eqnarray}\chi =\left(\begin{array}{c} {{\chi}^{11}} & {{\chi}^{12}} \\ {{\chi}^{21}} & {{\chi}^{22}} \end{array}\right),\end{eqnarray} \tag{ 104 }$$

and one can show that only the blocks ${{\chi}^{12}}$ and ${{\chi}^{21}}$ of the matrix of $\chi (t)$ are relevant for the dynamics (and they are related through ${{\chi}^{21}}={{\chi}^{12\,\dagger}}$). After some manipulations, one finds

$$\begin{eqnarray}&&\left(\begin{array}{c} {{\chi}^{12}} \\ {{\chi}^{12\,\ast}} \end{array}\right)={{\mathcal{M}}^{-1}}\left(\begin{array}{c} \overset{\centerdot}{{\mathcal{R}}}\,_{0}^{12} \\ \overset{\centerdot}{{\mathcal{R}}}\,_{0}^{12\,\ast} \end{array}\right),\end{eqnarray} \tag{ 105 }$$

where $\mathcal{M}$ is the same linear response matrix as in (88). Inserting this relationship in (103) we end up with an expression that is fully reminiscent of the kinetic energy

$$\begin{eqnarray}&&\mathcal{K}\equiv \frac{1}{2}\left({{\chi}^{12\,\dagger}},~{{\chi}^{12\,T}}\right)\mathcal{M}\left(\begin{array}{c} {{\chi}^{12}} \\ {{\chi}^{12\,\ast}} \end{array}\right)=\frac{1}{2}\left(\overset{\centerdot}{{\mathcal{R}}}\,_{0}^{12,\dagger},~\overset{\centerdot}{{\mathcal{R}}}\,_{0}^{12\,T}\right){{\mathcal{M}}^{-1}}\left(\begin{array}{c} \overset{\centerdot}{{\mathcal{R}}}\,_{0}^{12} \\ \overset{\centerdot}{{\mathcal{R}}}\,_{0}^{12\,\ast} \end{array}\right).\end{eqnarray} \tag{ 106 }$$

This expression shows that the linear response matrix is indeed the matrix of inertia.

As for the GCM, the above expressions are used in a framework where it is assumed that just a few collective coordinates are responsible for the time evolution of the system. This assumption allows the expression of the time derivative of the density matrix as

$$\begin{eqnarray}&&{{\overset{\centerdot}{{\mathcal{R}}}\,}_{0}}=\underset{\alpha}{\sum}\,\frac{\partial {{\mathcal{R}}_{0}}}{\partial {{q}_{\alpha}}}{{\overset{\centerdot}{{q}}\,}_{\alpha}},\end{eqnarray} \tag{ 107 }$$

where the ${{q}_{\alpha}}$ are the relevant collective coordinates and ${{\overset{\centerdot}{{q}}\,}_{\alpha}}$ their time derivatives. The kinetic energy (106) becomes

$$\begin{eqnarray}&&\mathcal{K}=\frac{1}{2}\underset{\alpha \beta}{\sum}\,{{M}_{\alpha \beta}}{{\overset{\centerdot}{{q}}\,}_{\alpha}}{{\overset{\centerdot}{{q}}\,}_{\beta}},\end{eqnarray} \tag{ 108 }$$

with

$$\begin{eqnarray}&&{{M}_{\alpha \beta}}=\frac{\partial {{\mathcal{R}}_{0}}}{\partial {{q}_{\alpha}}}{{\mathcal{M}}^{-1}}\frac{\partial {{\mathcal{R}}_{0}}}{\partial {{q}_{\beta}}},\end{eqnarray} \tag{ 109 }$$

which constitutes the expression of the collective inertia for the relevant collective coordinates ${{q}_{\alpha}}$.

As in the GCM case, the final expression of the ATDHFB mass involves the inverse of the linear response matrix. The only (very recent) attempt to invert this matrix explicitly has been reported in [222]. Most often, the problem is simplified by resorting to the same two approximations encountered earlier and discussed in detail in [216]: the 'cranking approximation' uses the diagonal approximation to the linear response matrix with two quasiparticle energies as diagonal elements. The partial derivatives of the density ${{\mathcal{R}}_{0}}$ with respect to the collective variables ${{q}_{\alpha}}$ are evaluated numerically by obtaining the HFB wave function with the constraints $\boldsymbol{q}+\delta \boldsymbol{q}$ and using finite difference approximations to the partial derivatives⁴. However, the partial derivatives of the generalized density can also be obtained by applying linear response theory, which leads to an explicit expression of the partial derivatives involving again the inverse of the linear response matrix. If the same cranking approximation is used as before, we then obtain the 'perturbative cranking' formula for the ATDHFB mass,

$$\begin{eqnarray}&&{{{\mathsf {M}}}_{\text{ATDHFB}}}={{\hslash}^{2}}{{\left[{{{\mathsf {M}}}^{(-1)}}\right]}^{-1}}{{{\mathsf {M}}}^{(-3)}}{{\left[{{{\mathsf {M}}}^{(-1)}}\right]}^{-1}},\end{eqnarray} \tag{ 110 }$$

with the moments defined in (93).

The validity of the perturbative cranking approximation (replacing exact derivatives by a linearization) has recently been tested in the calculation of fission pathways of superheavy elements in [223]. The figure 15 shows the fission pathways in ²⁶⁴Fm in the $\left({{q}_{20}},{{q}_{22}}\right)$ collective space obtained under different approximations: the path obtained by only considering the lowest energy is marked 'static'; the path obtained by minimizing the action while taking a constant inertia is marked 'const'; the paths obtained by minimizing the action with the collective inertia computed at the cranking approximation either directly or perturbatively are denoted by 'C+DPM', 'C^p +DPM' and 'C+RM', 'C^p +RM' respectively. There is a clear, qualitative difference between the perturbative and non-perturbative treatment of collective inertia.

Zoom In Zoom Out Reset image size

The WKB approximation is often used to get an estimate of the transmission coefficient through the fission barrier [197]. To get an idea of how it works, let us consider first the one-dimensional case. The basic idea is to substitute the classical expression for the momentum $p=\sqrt{2m\left(E-V(x)\right)}$ into the quantum mechanical identity

$$\begin{eqnarray}&&\frac{{{\text{d}}^{2}} \Psi (x)}{\text{d}{{x}^{2}}}=-\frac{{{p}^{2}}}{{{\hslash}^{2}}} \Psi (x).\end{eqnarray} \tag{ 111 }$$

This second-order differential equation is then solved under the assumption that the wave function can be expressed in the generic form

$$\begin{eqnarray} &&\Psi (x)=A(x){{\text{e}}^{\text{i}\phi (x)}}.\end{eqnarray}$$

Inserting this ansatz in (111) and assuming that the amplitude A(x) varies slowly with x one obtains

$$\begin{eqnarray}&&A(x)=\frac{C}{\sqrt{p(x)}},\end{eqnarray}$$

and

$$\begin{eqnarray}&&\phi (x)=\pm \frac{1}{\hslash}{\int}^{}\text{d}x\;p(x),\end{eqnarray}$$

which is valid in the 'classical' region where $E\geqslant V(x)$. The extension to classically-forbidden regions is straightforward, requiring the introduction of a complex momentum p(x) that leads to an exponential wave function in that region

$$\begin{eqnarray} &&\Psi (x)=\frac{{{C}_{\pm }}}{\sqrt{|\,p(x)|}}\exp \left(\pm \frac{1}{\hslash}{\int}^{}\text{d}x\sqrt{2m\left(V(x)-E\right)}\right).\end{eqnarray} \tag{ 112 }$$

The positive sign in the exponent yields an exponentially increasing tunnelling probability and therefore the corresponding amplitude C₊ has to be very small. By keeping only the term with the minus sign for the wave function in the classically-forbidden region, we obtain for the transmission coefficient

$$\begin{eqnarray}&&T=\,| \Psi (b){{|}^{2}}/| \Psi (a){{|}^{2}}=\exp \left(-\frac{2}{\hslash}{\int}_{a}^{b}\text{d}x\sqrt{2m\left(V(x)-E\right)}\right),\end{eqnarray} \tag{ 113 }$$

where a and b are the inner and outer turning points at the barrier corresponding to the energy E.

The extension of these ideas to fission is not entirely straightforward because of the infinite degrees of freedom of the nuclear many-body system. In the adiabatic approximation, however, the use of the WKB formula is somewhat justified owing to the reduction in the number of relevant collective variables. Clearly, the role of the mass of the particle in the one-dimensional problem recalled above should be played by the collective inertia tensor, and the potential energy should be replaced by the HFB energy along the collective path (possibly supplemented by quantum corrections such as the rotational energy correction, zero-point energies, etc). The actual path used to compute the transmission coefficient is determined by invoking the last action principle mentioned in the introduction to this chapter: The physical path is the one for which the classical action (67) is minimal. The spontaneous fission lifetime $\tau _{1/2}^{\text{SF}}$ is then usually given by the inverse of the product of the transmission probability T times the number of assaults to the barrier per unit time ν

$$\begin{eqnarray}&&\tau _{1/2}^{\text{SF}}=\frac{1}{\nu}\exp \left(\frac{2}{\hslash}{\int}_{a}^{b}\text{d}s\sqrt{2B(s)\left(V(s)-{{E}_{0}}\right)}\right).\end{eqnarray} \tag{ 114 }$$

Usually $1/\nu$ is estimated assuming that the zero point energy correction of the nucleus in its ground-state is of the order of one MeV, which leads to the value $1/\nu ={{10}^{-21}}$ s. Other authors (e.g. [20]) prefer to compute this parameter as the zero-point energy $\hslash \omega$ of the ground-state of the potential well corresponding to the collective variable q. For the mass, the two expressions for the collective inertia obtained previously in the framework of the GCM and ATDHFB methods are used, usually within the perturbative cranking approximation. Some results have recently been published in [223] without the perturbative treatment of the derivatives but, to our knowledge, there has been no calculation where the exact masses (involving the inverse of the full linear response matrix) have been used.

As expected, the results for $\tau _{1/2}^{\text{SF}}$ are very sensitive to the particular approximation used, especially for nuclei with large $\tau _{1/2}^{\text{SF}}$ values. This is a straightforward consequence of the exponentiation of the action, leading to an exponential dependence on the potential, energy and inertia tensor. In applications of the WKB method, the collective potential V(q) is usually supplemented by the various corrections discussed in detail in section 2.2.5. Strictly speaking, zero-point energy corrections (ZPE) should only be computed when the GCM framework is used to compute the collective inertia. Some authors have argued that even in the ATDHFB case, where no ZPE correction is present by construction, some sort of ZPE can still be used by taking the GCM form and replacing the GCM by the ATDHFB mass, see for instance [19, 20, 223]. In any case, it is commonly assumed as in [202] that the ZPE corrections associated to the quadrupole moment vary slowly with the collective variable and therefore only represent a displacement of the energy origin.

Finally, the E₀ parameter is often taken as the HFB ground-state energy. However it has been argued, e.g. in [168, 224], that the dynamics of the collective degree of freedom has an associated zero-point energy correction on top of the potential minimum (the HFB ground-state energy). This zero-point energy is often taken as a phenomenological parameter varying in the range 0.5–1 MeV. Some authors prefer to compute it using the formulas obtained in the GCM formalism [20].

Although the WKB method is the most popular choice to compute fission half-lives, several authors have considered alternative methods based on functional integrals, see for instance [225–231]. Such techniques offer, at least in principle, the possibility to be extended to arbitrary many-body systems without relying on an explicit choice of collective coordinates (and the underlying adiabaticity hypothesis). Qualitatively, it can be thought of as an extension of time-dependent density functional theory in the classically-forbidden region 'below the barrier'. In this section, we only summarize some of the main features of the theory by following [227] in the simplified case of a one-dimensional problem.

We thus assume that the potential energy of the system has the typical fission-like features shown in figure 16. We then adopt a functional integral representation of the quantum-mechanical evolution operator $\hat{U}$, and following [227] write

$$\begin{eqnarray}\begin{array}{*{35}{l}} \underset{n}{\sum}\,\frac{1}{E-{{E}_{n}}+\text{i}\epsilon} & =\text{Tr}\left[\frac{1}{\hat{H}-E+\text{i}\epsilon}\right] \end{array}\end{eqnarray} \tag{ 115 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & =-\text{i}{\int}_{0}^{\infty}\text{d}T\;{{\text{e}}^{\text{i}ET}}{\int}^{}\text{d}q\;\langle q|{{\text{e}}^{-\text{i}\hat{H}T}}|q\rangle \end{array}\end{eqnarray} \tag{ 116 }$$

with q the one dimensional coordinate, and the summation in the left-hand side extends over all eigenstates of the Hamiltonian. The poles of the resolvent as a function of E give the energy of bound states and resonances. In the classically-forbidden region, no bound states are possible and the resonances appear as complex energies with an imaginary part providing the width $\Gamma$ of the state.

Zoom In Zoom Out Reset image size

In order to estimate the integral over q, the integrand is first converted to a Feynman path integral

$$\begin{eqnarray}&&{\int}^{}\text{d}q\langle q|{{\text{e}}^{-\text{i}\hat{H}T}}|q\rangle ={\int}^{}\text{d}q{\int}^{}\mathcal{D}\left[q(t)\right]{{\text{e}}^{\text{i}S\left[q(t)\right]}},\end{eqnarray} \tag{ 117 }$$

where S[q(t)] is the classical action. This Feynman path integral is then computed in the static path approximation (SPA). The SPA picks out only the path q₀(t) that is periodic both with respect to the coordinate q and the momentum p, q(0)  =  q(T) and p(0)  =  p(T), and corresponds to a minimization of the classical action. The remaining integral over T is computed using again the SPA, which provides the relation $E=-\partial S/\partial T$ connecting E with the classical energy. The contribution to the resolvent is proportional (the proportionality factor depends on second order corrections to the SPA that we do not discuss here) to ${{\text{e}}^{\text{i}W(E)}}$ where $W(E)=E{{T}_{\text{cl}}}+{{S}_{\text{cl}}}\left({{T}_{\text{cl}}}\right)$ and T_cl is the classical period of motion for an energy E. Summing over all integer multiples of the period at an energy E gives the contribution ${{\text{e}}^{\text{i}W(E)}}/\left(1-{{\text{e}}^{\text{i}W(E)}}\right)$ to the resolvent for the orbits in the classically-allowed region around the local minimum. The poles of this quantity are at those energies satisfying $W(E)=2\pi n$, which is the WKB quantization energy formula up to a factor π. As argued in [226], better treatment of the omitted factor in front of the phase provides the missing π. In the classically-forbidden region, we can repeat the same kind of arguments by going into imaginary time, $t\to \text{i}\tau$. Still neglecting the quadratic corrections to the SPA, the total contribution of all paths to the resolvent is then given in [225] by

$$\begin{eqnarray}&&\text{Tr}\left[\frac{1}{\hat{H}-E+\text{i}\epsilon}\right]=\frac{{{\text{e}}^{\text{i}{{W}_{1}}(E)}}-{{\text{e}}^{-{{W}_{2}}(E)}}}{1-{{\text{e}}^{\text{i}{{W}_{1}}(E)}}-{{\text{e}}^{-{{W}_{2}}(E)}}},\end{eqnarray} \tag{ 118 }$$

where W₁(E) is the action in the classical allowed region and W₂(E) is its generalization in the classically-forbidden region,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{W}_{1}}(E)=2{\int}_{{{q}_{1}}}^{{{q}_{2}}}\text{d}q\sqrt{2m\left(E-V(q)\right)}, \\ {{W}_{2}}(E)=2{\int}_{{{q}_{2}}}^{{{q}_{3}}}\text{d}q\sqrt{2m\left(V(q)-E\right)}. \end{array}\end{eqnarray} \tag{ 119 }$$

The expression of the trace is slightly different if corrections are taken into account and is given in [227]. The poles of the resolvent are now the solutions of $1-{{\text{e}}^{\text{i}{{W}_{1}}(E)}}-{{\text{e}}^{-{{W}_{2}}(E)}}=0$ and therefore lie in the complex plane. Assuming that the last term is small, the complex energy solutions are given by ${{E}_{n}}+\frac{\text{i}}{2}{{ \Gamma }_{n}}$ where E_n is the solution of the WKB energy quantization condition. The width of the resonance is given by

$$\begin{eqnarray}&&{{ \Gamma }_{n}}=\frac{{{\omega}_{\text{cl}}}\left({{E}_{n}}\right)}{\pi}{{\text{e}}^{-{{W}_{2}}\left({{E}_{n}}\right)}},\end{eqnarray} \tag{ 120 }$$

which is a factor of 2 smaller than the WKB formula in the large W₂ limit. Again, this factor is recovered when including quadratic corrections, as argued in [227]. It can also be shown that the region to the right of the barrier does not contribute significantly to the resolvent.

These ideas can in principle be extended to a quantum many-body system as discussed in [226, 227, 232]. The peculiarity of this approach is that the SPA to the Feynman path integral corresponds to the mean field solution instead of the classical trajectories. Aside from that, the method requires solving many TDHFB-like equations in imaginary time with periodic boundary conditions to describe the dynamics below the barrier. A non trivial and still unresolved issue is how to connect periodic trajectories in classically-allowed regions to the solutions under the barrier. In one dimension there is only one way to do so but in a multidimensional case the connection is far from trivial. Early studies in [229] point to the important role of symmetry-breaking. The connection to the standard WKB formula is still missing although the results of [232] seem to suggest that the WKB formula might grasp some of the physics involved.

Although not directly connected with functional integral methods, the tunnelling through a multidimensional barrier has also been studied in a two-dimensional model using a semi-classical approximation with complex classical trajectories [233]. Finally, we also mention the recent attempt in [231] to compute tunnelling probabilities by considering a complex absorbing potential: although the method was tested on simple fission model Hamiltonians, it could in principle be extended to more microscopic collective Hamiltonians of the form (131).

The simplest and most drastic of all approximations is to forgo the quantum nature of the nucleus and treat fission dynamics in a semi-classical way. We do not intend to give a comprehensive description of the various stochastic methods used to describe nuclear dynamics, and refer the reader to the review [234] by Abe et al where this topic is discussed in great detail. Here, our goal is simply to recall how some of these techniques have been applied to describe induced fission, especially fission fragment distributions, particle evaporation and fission probabilities.

We introduce the conjugate momenta $\boldsymbol{p}$ of the collective variables $\boldsymbol{q}$ driving fission. At any time t, the nucleus is represented by a point in phase space with coordinates $\left(\boldsymbol{q}(t),\boldsymbol{p}(t)\right)$ where the potential energy is $V\left(\boldsymbol{q}\right)$. If the total energy of the nucleus is E, then the local excitation energy is ${{E}^{\ast}}\left(\boldsymbol{q}\right)=E-V\left(\boldsymbol{q}\right)$ (note that for such a classical treatment of nuclear dynamics, the excitation energy must be positive ${{E}^{\ast}}\left(\boldsymbol{q}\right)\geqslant 0$ for all $\boldsymbol{q}$). The dynamics of the system can be represented in several ways:

• The Langevin equations directly give the position of the system in phase space at any time t. They read
  $$\begin{eqnarray}\begin{array}{*{35}{l}} {{\overset{\centerdot}{{q}}\,}_{\alpha}} & =\underset{\beta}{\sum}\,{{B}_{\alpha \beta}}{{p}_{\beta}}, \end{array}\end{eqnarray} \tag{ 122 }$$
  $$\begin{eqnarray}\begin{array}{*{35}{l}} {{\overset{\centerdot}{{p}}\,}_{\alpha}} & =-\frac{1}{2}\underset{\beta \gamma}{\sum}\,\frac{\partial {{B}_{\beta \gamma}}}{\partial {{q}_{\alpha}}}{{p}_{\beta}}{{p}_{\gamma}}-\frac{\partial V}{\partial {{q}_{\alpha}}}-\underset{\beta \gamma}{\sum}\,{{ \Gamma }_{\alpha \beta}}{{B}_{\beta \gamma}}{{p}_{\gamma}}+\underset{\beta}{\sum}\,{{ \Theta }_{\alpha \beta}}{{\xi}_{\beta}}(t), \end{array}\end{eqnarray} \tag{ 123 }$$
  with ${\mathsf {B}}\left(\boldsymbol{q}\right)\equiv {{B}_{\alpha \beta}}\left(\boldsymbol{q}\right)$ the tensor of inertia, ${\mathsf {\Gamma }}\left(\boldsymbol{q}\right)\equiv {{ \Gamma }_{\alpha \beta}}\left(\boldsymbol{q}\right)$ the coordinate-dependent friction tensor (not to be confused with the width in the GCM) and $V\left(\boldsymbol{q}\right)$ the potential energy in the collective space. The Langevin equations are non-deterministic owing to the presence of the random force $\boldsymbol{\xi }(t)$. The strength of this random force is controlled by the parameter ${\mathsf {\Theta }}\equiv {{ \Theta }_{\alpha \beta}}$. In applications of the Langevin equation to fission such as, e.g. in [63, 235, 236], this parameter is usually related to the friction tensor through the fluctuation-dissipation theorem, ${\sum}_{k}{{ \Theta }_{ik}}{{ \Theta }_{kj}}={{ \Gamma }_{ij}}T$, with $T\equiv T\left(\boldsymbol{q}\right)$ a local nuclear temperature related to the excitation of the system ${{E}^{\ast}}\left(\boldsymbol{q}\right)$ at point q. Furthermore, it is often assumed that the random variable is a Gaussian white noise process characterized by
  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle {{\xi}_{\alpha}}\rangle =0, \end{array}\end{eqnarray} \tag{ 124 }$$
  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle {{\xi}_{\alpha}}(t){{\xi}_{\beta}}\left({{t}^{\prime}}\right)\rangle =2{{\delta}_{\alpha \beta}}\delta \left(t-{{t}^{\prime}}\right), \end{array}\end{eqnarray} \tag{ 125 }$$
  where in this equation, $\langle.\rangle$ refer to statistical averaging. This absence of memory for $\boldsymbol{\xi }$ implies that the Langevin equations represent a Markovian stochastic process, i.e. the value of the random force at time t does not depend on previous values at time ${{t}^{\prime}}<t$; see [137] for an introduction. It can be thought of as a random walk on the PES defined by $V\left(\boldsymbol{q}\right)$. Each such random walk from some initial state to a properly defined scission configuration defines a fission event; repeating the procedure, e.g. by Monte-Carlo sampling, allows one to reconstruct the full distribution of fission fragments.
• The Kramers equation gives the probability distribution function for the nucleus to be at a given point in phase space. It reads
  $$\begin{eqnarray}\begin{array}{*{35}{l}} \overset{\centerdot}{{f}}\,\left(\boldsymbol{q},\boldsymbol{p},t\right) & =\underset{\alpha \beta}{\sum}\,\left[-{{B}_{\alpha \beta}}{{p}_{\beta}}\frac{\partial}{\partial {{q}_{\alpha}}}+\left(\frac{\partial V}{\partial {{q}_{\alpha}}}+\underset{\delta}{\sum}\,\frac{1}{2}\frac{\partial {{B}_{\beta \gamma}}}{\partial {{q}_{\alpha}}}{{p}_{\beta}}{{p}_{\delta}}\right)\frac{\partial}{\partial {{p}_{\alpha}}} \right. \\ {} & \left. +\underset{\gamma}{\sum}\,{{ \Gamma }_{\alpha \beta}}{{B}_{\beta \gamma}}\frac{\partial}{\partial {{p}_{\alpha}}}{{p}_{\delta}}+{{D}_{\alpha \beta}}\frac{{{\partial}^{2}}}{\partial {{p}_{\alpha}}\partial {{p}_{\beta}}}\right]f\left(\boldsymbol{q},\boldsymbol{p},t\right), \end{array}\end{eqnarray} \tag{ 126 }$$
  where ${{D}_{\alpha \beta}}={{ \Gamma }_{\alpha \beta}}T$ and T is the temperature. Since the Kramers equation does not contain an explicit random term, it lends itself to analytic approximations. On the other hand, its generalization to N-dimensional collective spaces is numerically more involved since it is a second-order differential equation in terms of collective variables.

Historically, the Kramers equation was essentially used to extract an analytic expression for the induced fission width ${{ \Gamma }_{f}}$ by considering fission as a diffusion process over a potential barrier approximated by an N-dimensional quadratic surface (of the type $\frac{1}{2}{\sum}_{ij}{{q}_{i}}{{q}_{j}}$ where q_i are the collective variables), for instance in [237–239]. More recently, progress in computing capabilities has enabled solving the full Langevin equations in several dimensions. Results reported by various groups differ mostly in the number and type of collective degrees of freedom, as well as the prescription for the friction tensor ${\mathsf {\Gamma }}$. Most studies have been performed in the high temperature regime, where the potential energy surface is approximated by a liquid-drop-like formula, for example [235, 240–242]. The Langevin equations have also been solved for macroscopic-microscopic PES in the limiting case of strong friction (strongly damped Brownian motion) in [63, 64, 243]. Note that both the Langevin and Kramers equations involve the collective potential energy $V\left(\boldsymbol{q}\right)$ and the inertia tensor ${\mathsf {B}}\left(\boldsymbol{q}\right)$. In the recent work of [244], these quantities were computed using the ATDHFB formula, and Langevin dynamics was solved from the outer turning point to scission in order to extract spontaneous fission fragment distributions.

Time-dependent density functional theory (TDDFT) provides a fully microscopic approach to describe real time fission dynamics. It is a reformulation of the many-body time-dependent Schrödinger equation as discussed in [245, 246]. If one considers an interacting electron gas in a time-dependent external potential, then the Runge–Gross existence theorem of [247] asserts that given an initial state all properties of the system can be expressed as a functional of the (time-dependent) local one-body density, providing the potential satisfies certain regularity conditions. Just as in the static case, the Kohn–Sham scheme can in principle be applied so that the TDDFT equations turn into simple time-dependent Hartree-like equations. However, again as in the static case, the form of the exchange-correlation time-dependent potential ${{\hat{v}}_{xc}}\left(\boldsymbol{r},t\right)$ is not known. What is called the adiabatic approximation in TDDFT (not to be confused with the adiabatic approximation in fission theory) consists in assuming that ${{\hat{v}}_{xc}}\left(\boldsymbol{r},t\right)$ has the same functional dependence on the density as at t  =  0,

$$\begin{eqnarray}&&\hat{v}_{xc}^{\text{adiabatic}}\left(\boldsymbol{r},t\right)={{\hat{v}}_{xc}}[\rho ]\left(\boldsymbol{r},t\right).\end{eqnarray} \tag{ 127 }$$

As for the Hohenberg–Kohn theorem of static DFT, there is no direct analogue of the Runge–Gross theorem for self-bound nuclear systems characterized by symmetry-breaking intrinsic densities; see also [248] for a discussion of self-bound systems with symmetry-conserving internal densities. In spite of this, the popular time-dependent Hartree–Fock (TDHF) and time-dependent Hartree–Fock–Bogoliubov (TDHFB) are de facto adaptations of adiabatic TDDFT in nuclear physics, and we give below a very brief presentation of each of these methods.

The TDHF equation

$$\begin{eqnarray}&&\text{i}\hslash \overset{\centerdot}{{\rho}}\,=[h,\rho ]\end{eqnarray} \tag{ 128 }$$

can be obtained by enforcing that the many-body wave-function remains a Slater determinant at all times [61]. Given an initial density at time t₀, which is typically obtained by solving the static HF equations under a set of constraints, solving the TDHF equation provides the full time-evolution of the system. If the initial condition is such that the system has enough excitation energy—a point discussed in detail in [249, 250]—this time evolution may follow the system past the scission point and lead to two separated, excited fission fragments. In some way, TDHF is the microscopic analogue of the Langevin equation in that it simulates a single fission event in real time. One of the earliest applications of TDHF in [251] was made by Negele and collaborators to study the fission of actinides. At the time, a number of approximations were needed such as axial and reflection symmetry, no spin–orbit potential, and a coarse spatial grid. Progress in computing has enabled more realistic simulations including the full Skyrme potential and 3D geometries as in [183]. In all cases, the initial point must have a deformation larger than the 'dynamical fission threshold' introduced in [250] in order for the system to fission. The existence of such a threshold was explained by Bulgac and collaborators in [182] as the consequence of neglecting pairing correlations. Because TDHF does not rely on the hypothesis of adiabaticity, it is expected to give a much more realistic description of the scission point, in particular of fission fragment properties. The initial results on fragment total kinetic and excitation energies reported in [182, 183] are very promising.

As has already been emphasized several times, pairing correlations play a crucial role in fission. While, in principle, the TDDFT equations could provide the exact time-evolution of the system with only a functional of the density ρ, our ignorance of the form of this functional forces us to introduce explicitly a Kohn–Sham scheme based on symmetry-breaking reference states and a (time-dependent) pairing tensor κ. This problem is identical to the static case. Nuclear dynamics is then described by the TDHFB equation,

$$\begin{eqnarray}&&\text{i}\hslash \overset{\centerdot}{{\mathcal{R}}}\,=\left[\mathcal{H},\mathcal{R}\right].\end{eqnarray} \tag{ 129 }$$

While formally analogous to the TDHF equations, the TDHFB equation is substantially more involved numerically. Indeed, the number of (partially) occupied orbitals at each time t is much larger than the number of nucleons. In spherical symmetry, the TDHFB equation has recently been solved without specific approximations [252]. In deformed nuclei, it has been solved in the canonical basis in [253–255], but no application of this formalism to fission has been performed yet. The first pioneer calculation of neutron-induced fission with full TDHFB, which the authors refer to as time-dependent local density approximation, has also been reported in [182]. We note that, as for the static problem, the TDHFB equation can be approximated by the simpler TDHF +BCS limit as done, for instance in [170, 256]. However, the TDHF +BCS approach does not respect the continuity equation, contrary to TDHFB, which may lead to non-physical results in specific cases such as particle emission. The authors of [256] advocated using a simplified version of TDHF +BCS where occupation numbers are frozen to their initial value and do not change with time.

If the TDHF equations are the analogue of the classical Langevin equations, the time-dependent generator coordinate method (TDGCM) can be viewed as the microscopic translation of the Kramers equation (126). The TDGCM is a straightforward extension of the static GCM introduced in section 3.1.1, where the ansatz for the solution to the time-dependent many-body Schrödinger equation takes the form

$$\begin{eqnarray}&&| \Psi (t)\rangle ={\int}^{}\text{d}\boldsymbol{q}\;f\left(\boldsymbol{q},t\right)| \Psi \left(\boldsymbol{q}\right)\rangle.\end{eqnarray} \tag{ 130 }$$

As with (69), the functions $| \Psi \left(\boldsymbol{q}\right)\rangle$ are known many-body states parametrized by a vector of collective variables $\boldsymbol{q}$, which most often are chosen as the solutions to the static HFB equations under a set of constraints $\boldsymbol{q}$, see also section 2.3 for a discussion of collective variables.

Inserting the ansatz (130) in the variational principle (121) provides the time-dependent analogue of (70), where the only difference is that the functions $f\left(\boldsymbol{q},t\right)$ are now time-dependent. All applications of the TDGCM method have been made by further assuming the Gaussian overlap approximation for the norm kernel. Generalizing the procedure given in section 3.1.1 in the case where the overlap kernel does not depend on the collective variable (constant metric), we find the time-dependent collective Schrödinger equation

$$\begin{eqnarray}&&\text{i}\hslash \frac{\partial}{\partial t}g\left(\boldsymbol{q},t\right)=\left[-\frac{{{\hslash}^{2}}}{2}\underset{\alpha \beta}{\sum}\,\frac{\partial}{\partial {{q}_{\alpha}}}{{B}_{\alpha \beta}}\left(\boldsymbol{q}\right)\frac{\partial}{\partial {{q}_{\beta}}}+V\left(\boldsymbol{q}\right)\right]g\left(\boldsymbol{q},t\right),\end{eqnarray} \tag{ 131 }$$

where the function $g\left(\boldsymbol{q},t\right)$ is related to the weight function $f\left(\boldsymbol{q},t\right)$ of (130) according to

$$\begin{eqnarray}&&g\left(\boldsymbol{q},t\right)={\int}^{}\text{d}\boldsymbol{q}^{{\prime}}\;f\left(\boldsymbol{q}^{{\prime}},t\right){{\left[n\left(\boldsymbol{q},\boldsymbol{q}^{{\prime}}\right)\right]}^{1/2}},\end{eqnarray} \tag{ 132 }$$

(see (73) for the definition of the square root of the norm) and contains all the information about the dynamics of the system. As before, the rank 2 tensor ${\mathsf {B}}\left(\boldsymbol{q}\right)\equiv {{B}_{\alpha \beta}}\left(\boldsymbol{q}\right)$ is the collective inertia of the system in the collective space, and $V\left(\boldsymbol{q}\right)$ is the potential energy.

Equation (131) implies a continuity equation for the quantity $|g\left(\boldsymbol{q},t\right){{|}^{2}}$,

$$\begin{eqnarray}&&\frac{\partial}{\partial t}|g\left(\boldsymbol{q},t\right){{|}^{2}}=-{{\nabla}_{\boldsymbol{q}}}\centerdot \boldsymbol{J}\left(\boldsymbol{q},t\right).\end{eqnarray} \tag{ 133 }$$

This equation can be derived in perfect analogy with standard one-body quantum mechanics, see for example the derivation in [196]. It suggests that $|g\left(\boldsymbol{q},t\right){{|}^{2}}$ can be interpreted as a probability amplitude for the system to be characterized by the collective variables $\boldsymbol{q}$ at time t. Consequently, the vector $\boldsymbol{J}\left(\boldsymbol{q},t\right)$ is a current of probability in perfect analogy with the results of one-body quantum mechanics (see, e.g. equation (IV.9) in [196]),

$$\begin{eqnarray}&&\boldsymbol{J}\left(\boldsymbol{q},t\right)=\frac{\hslash}{2\text{i}}{\mathsf {B}}\left(\boldsymbol{q}\right)\left[{{g}^{\ast}}\left(\boldsymbol{q},t\right){{\nabla}_{\boldsymbol{q}}}g\left(\boldsymbol{q},t\right)-g\left(\boldsymbol{q},t\right){{\nabla}_{\boldsymbol{q}}}{{g}^{\ast}}\left(\boldsymbol{q},t\right)\right].\end{eqnarray} \tag{ 134 }$$

When more than one collective variable are involved, the coordinates of the current of probability are

$$\begin{eqnarray}&&{{J}_{\alpha}}\left(\boldsymbol{q},t\right)=\frac{\hslash}{2\text{i}}\underset{\beta =1}{\overset{N}{\sum}}\,{{B}_{\alpha \beta}}\left(\boldsymbol{q}\right)\left[{{g}^{\ast}}\left(\boldsymbol{q},t\right)\frac{\partial g}{\partial {{q}_{\beta}}}\left(\boldsymbol{q},t\right)-g\left(\boldsymbol{q},t\right)\frac{\partial {{g}^{\ast}}}{\partial {{q}_{\beta}}}\left(\boldsymbol{q},t\right)\right].\end{eqnarray} \tag{ 135 }$$

Therefore, like the classical Kramers equations, the TDGCM equation give the evolution of the flow of probability in the collective space. The probability amplitude $|g\left(\boldsymbol{q},t\right){{|}^{2}}$ and the current (135) are the key quantities to extract fission fragment distributions in the TDGCM + GOA approach to nuclear fission. Based on the adequate identification of scission configurations, see discussion in section 2.4 p 21, it is possible to estimate the probability of a given scission configuration at point $\boldsymbol{q}$ by simply calculating the integrated flux of the probability current through the scission hyper-surface at that same point $\boldsymbol{q}$. If we define the integrated flux $F(\xi,t)$ through an oriented surface element ξ as

$$\begin{eqnarray}&&F(\xi,t)={\int}_{t=0}^{t}\text{d}T{{{\int}^{}}_{\boldsymbol{q}\in \xi}}\boldsymbol{J}\left(\boldsymbol{q},t\right)\centerdot \text{d}\boldsymbol{S}.\end{eqnarray} \tag{ 136 }$$

then, following [171, 257, 258], the fission fragment mass yield for mass A is

$$\begin{eqnarray}&&y(A)\propto \underset{\xi \in \mathcal{A}}{\sum}\,\underset{t\to +\infty}{{\lim}}\,F(\xi,t),\end{eqnarray} \tag{ 137 }$$

where $\mathcal{A}$ is the set of all oriented hyper-surface elements ξ belonging to the scission hyper-surface such that one of the fragments has mass A. In practice, fission fragment mass yields are normalized,

$$\begin{eqnarray}&&Y(A)=\frac{y(A)}{\underset{A}{\sum}\,y(A)}\end{eqnarray} \tag{ 138 }$$

Basis expansion techniques are ubiquitous in practical applications of quantum mechanics and quantum many-body theory. Expanding the solutions to the Schrödinger or Dirac equation on a basis of known functions yields a linear eigenvalue problem that can be solved very efficiently. In particular, the method is oblivious to the local or non-local character of the underlying nuclear potential. Moreover, the formulation of beyond mean-field extensions such as, e.g. the generator coordinate method or projection techniques is straightforward.

In nuclear science, the eigenfunctions of the one-centre harmonic oscillator (HO) have historically played a special role. They are known analytically in spherical, cylindrical, Cartesian coordinates, among others. Talmi, Moshinski and Talman showed long ago in [261–263] that any product of two HO basis functions could be expanded into a sum of single HO functions, which allows for the exact separation of the centre of mass and relative motion for two-body potentials. Using the harmonic oscillator basis also greatly simplifies the calculation of matrix elements of Gaussian potentials, such as, e.g. the Gogny force as highlighted in [102, 264]. Special care must be taken in the evaluation of matrix elements of states with large quantum numbers, e.g. by recurring to known properties of hypergeometric sums as in [167]; see also [265] for a recent account on the calculation of matrix elements of the Gogny force with axially symmetric harmonic oscillator wave functions. High-accuracy expansions of the Coulomb or Yukawa potentials onto a finite sum of Gaussians were also used in [87, 266, 267] for the precise calculation of the Coulomb exchange contribution to the nuclear mean field. Last but not least, the nuclear mean field can be well approximated by a HO, which is one of the reasons for the spectacular success of the phenomenological Nilsson model of nuclear structure [11].

Several DFT solvers based on expanding the HFB solutions on the HO basis have been used in fission studies. Let us mention in particular the codes hfbtho ([268, 269]) and hfodd ([160, 267, 270–274]), which have been released under open source license. Both codes implement generalized Skyrme-like functionals in the particle–hole channel and density-dependent delta-interaction in the particle–particle channel. The latest version of hfodd also implements EDF based on finite-range pseudopotentials such as the Gogny force in both channels. Axial and time-reversal symmetries are built-in in hfbtho, but reflection-asymmetric shapes are possible; by contrast, hfodd breaks all possible symmetries of the nuclear mean-field. The two codes have been carefully benchmarked against one another, and the hfbtho kernel is included as a module of hfodd. These codes were used to study both spontaneous and induced fission; see the following papers [15–17, 20, 129, 162, 186, 216, 223, 275–277]. In applications with the Gogny forces, the codes hfbaxial and hfbtri have been widely used to study spontaneous fission in actinides and superheavy elements [128, 168, 172, 278, 279].

Let us recall that the three-dimensional quantum HO is characterized by its frequency vector $\boldsymbol{\omega }=\left({{\omega}_{x}},{{\omega}_{y}},{{\omega}_{z}}\right)$ (in Cartesian coordinates). The frequency, measured in MeV, is also related to the oscillator length ${{b}_{\mu}}=\sqrt{\hslash /m{{\omega}_{\mu}}}$ (in fm), where m is the mass of a nucleon. While the single particle Hilbert space is of course infinite, practical implementations require truncating the HO basis. This is achieved in different ways. For a spherical basis with ${{\omega}_{x}}={{\omega}_{y}}={{\omega}_{z}}\equiv {{\omega}_{0}}$, one usually imposes a cut-off in the number ${{N}_{\text{shell}}}$ of oscillator shells. Each N-shell contains (N +1)(N +2) degenerate states [11, 56]. Deformed or 'stretched' HO bases are introduced to describe elongated geometries. In these cases, ${{\omega}_{x}}\ne {{\omega}_{y}}\ne {{\omega}_{z}}$, but the condition of volume conservation yields ${{\omega}_{x}}{{\omega}_{y}}{{\omega}_{z}}=\omega _{0}^{3}$ and accordingly $b_{0}^{3}={{b}_{x}}{{b}_{y}}{{b}_{z}}$. Since there is no degeneracy of the HO shells any more, one must introduce additional criteria to truncate the basis. Among the popular choices are the ratios $p={{\omega}_{x}}/{{\omega}_{y}}$ and $q={{\omega}_{x}}/{{\omega}_{z}}$, which define the deformation of the basis. Alternatively, this basis deformation can also be defined by introducing an ellipsoidal liquid drop characterized by the $(\beta,\gamma )$ Bohr deformations; see [61] for the relation between $(\beta,\gamma )$ and the $\left({{\alpha}_{20}},{{\alpha}_{22}}\right)$ parameters of the expansion (2). The quantities p and q can then be expressed as ratios of the radii of each of the principal axes of this ellipsoid, $p={{R}_{y}}/{{R}_{x}}$ and $q={{R}_{z}}/{{R}_{x}}$, each radius being a function of $(\beta,\gamma )$. This method was introduced in [87] and generalized in [270]. An additional number of states ${{n}_{\text{max}}}$ is sufficient to completely determine the basis states.

As important as the truncation schemes are the choices for the oscillator lengths used in the HO basis. Typically, these quantities are used as additional variational parameters that are adjusted to minimize the energy. This search for optimal oscillator lengths is obviously less important the bigger the size of the basis is, as illustrated in figure 17. Therefore, there is always a compromise between using a large basis where the precise value of the oscillator lengths is less relevant but calculations are expensive, or using a smaller basis at the cost of repeating, for the configuration of interest, the HFB calculation several times for different oscillator lengths. In this respect various phenomenological formulas relating the oscillator lengths to the imposed deformation parameters of the nucleus can be used as in [162].

Zoom In Zoom Out Reset image size

Because of the truncation of the basis, the solution of the HFB equation becomes dependent on the characteristics of the basis, namely ${{N}_{\text{shell}}}$, ${{\omega}_{0}}$, and q in the most common case of an axially-deformed basis. This dependence is clearly spurious and disappears in the limit of an infinite basis. In practical calculations, however, its effect must be properly quantified. The figure 17 illustrates the expected size of truncation effects as a function of the oscillator frequencies and maximum spherical shell number. We show the energy of a deformed configuration along the fission path of ²⁴⁰Pu characterized by $\langle {{\hat{Q}}_{20}}\rangle =200b$ and $\langle {{\hat{Q}}_{40}}\rangle =50{{b}^{2}}$. Even at ${{N}_{\text{shell}}}=24$, the energy varies by several hundreds of keV over the range ${{b}_{0}}\in [1.9,2.6]$ fm..

The work reported in [88, 161, 168, 277] showed the impact of basis truncation on fission properties, mostly on the static properties of the potential energy surface such as, e.g. the height of fission barriers. One should bear in mind that truncation errors typically amount to a few hundred keV at the top of the first fission barrier in actinides. Such errors can cause several orders of magnitude changes in spontaneous fission half-lives because of the exponential factor in (114). In addition, the truncation error increases with the mass of the nucleus, especially when pairing correlations are non-zero. Because of the Gaussian asymptotic behaviour of HO wave functions, the convergence of basis expansions in weakly bound nuclei near the dripline is also problematic, see possible alternatives for DFT solvers in, e.g. [280, 281]. This could have a major impact in fission fragment calculations of nuclei involved in the r-process discussed in [282, 283]. Finally, we already mentioned that both spontaneous fission at high excitation energy and neutron-induced fission are often described in the finite-temperature HFB theory, where the density matrix contains a spatially non-localized component. This increases truncation errors accordingly. In figure 18, we show the evolution of the energy as a function of the oscillator length for FT-HFB calculations at T  =  1.0, 1.5, 2.0 MeV for large bases characterized by shell numbers ${{N}_{\text{shell}}}=20$ and ${{N}_{\text{shell}}}=28$. We note that the plateau condition for convergence degrades substantially as nuclear temperature increases: at T  =  1.0 MeV, the energy does not change by more than 50 keV over the range ${{b}_{0}}\in [2.2,2.4]$ fm, while at T  =  2.0 MeV, even a small change of $\delta b=0.05$ fm around the minimum induces variations of energy of 200 keV.

Zoom In Zoom Out Reset image size

In the vast majority of DFT solvers based on the HO basis expansion, the basis functions are the eigenstates of the one-centre, three-dimensional quantum HO. In the nineteen eighties, the French group at CEA Bruyères-le-Châtel developed an axial two-centre DFT solver, in which the basis functions are superpositions of the eigenfunctions of two one-centre HOs shifted by a distance d (adjustable). By construction, the set of all such functions is not orthogonal, which requires a Gram–Schmidt orthogonalization procedure. This represents a disadvantage as the number of linearly independent states depends upon the distance d and therefore wave functions at different elongations are expanded in different sub-spaces of the full Hilbert space. As a consequence, the evolution of observables in the transition from one subspace to the neighbouring ones is not necessarily smooth. On the other hand, this technology is especially advantageous to describe very deformed nuclear shapes and/or two fragments such as the ones encountered near scission. This code has been used to study induced fission in actinide nuclei, see [120, 123, 164, 188, 258, 284] for instance.

The coordinate space formulation of the HFB equation provides one of the simplest methods to remedy the limitations of basis expansion techniques. In this case, the HFB equation (24) becomes

$$\begin{eqnarray}{\int}^{}{{\text{d}}^{3}}\boldsymbol{r}^{{\prime}}\underset{{{\sigma}^{\prime}}}{\sum}\,\left(\begin{array}{c} h\left(\boldsymbol{r}\sigma,\boldsymbol{r}{{\sigma}^{\prime}}\right)-\lambda {{\delta}_{\sigma {{\sigma}^{\prime}}}} & \Delta \left(\boldsymbol{r}\sigma,\boldsymbol{r}{{\sigma}^{\prime}}\right) \\ -{{ \Delta }^{\ast}}\left(\boldsymbol{r}\sigma,\boldsymbol{r}{{\sigma}^{\prime}}\right) & -{{h}^{\ast}}\left(\boldsymbol{r}\sigma,\boldsymbol{r}{{\sigma}^{\prime}}\right)+\lambda {{\delta}_{\sigma {{\sigma}^{\prime}}}} \end{array}\right)\left(\begin{array}{c} {{U}_{\mu}}\left(\boldsymbol{r}{{\sigma}^{\prime}}\right) \\ {{V}_{\mu}}\left(\boldsymbol{r}{{\sigma}^{\prime}}\right) \end{array}\right)\\ &&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,={{E}_{\mu}}\left(\begin{array}{c} {{U}_{\mu}}\left(\boldsymbol{r}\sigma \right) \\ {{V}_{\mu}}\left(\boldsymbol{r}\sigma \right) \end{array}\right),\end{eqnarray} \tag{ 139 }$$

with the mean field $h\left(\boldsymbol{r}\sigma,\boldsymbol{r}^{{\prime}}{{\sigma}^{\prime}}\right)$ and pairing field $\Delta \left(\boldsymbol{r}\sigma,\boldsymbol{r}^{{\prime}}{{\sigma}^{\prime}}\right)$ expressed in coordinate space, see for instance [79] for transformation rules between configuration and coordinate space. While the separation of the Schrödinger equation makes this direct approach very straightforward if spherical symmetry is conserved (see applications in [89, 110]), extensions to deformed nuclei require introducing lattice techniques, as finite differences become either numerically too inaccurate or computationally impractical. There are several different examples of coordinate-space approaches to solving the HFB equation that have been applied to fission studies.

The HFB equation (139) is a particular example of an integro-differential equation. The B-spline collocation method (BSCM) was introduced in nuclear theory more than two decades ago in [285, 286] to solve such equations. In practice, the BSCM has been successfully implemented in cylindrical coordinates for the particular case of axially-deformed nuclei. After discretizing the spatial domain as a set of knots, one introduces a set of interpolating B-spline functions with order M in each knot. Any arbitrary function or operator is then represented at the collocation points defined by the maximum of the spline functions at each knot. With this technique, the HFB equation takes the form of the standard non-linear eigenvalue problem, which can be solved iteratively by successive diagonalizations. This approach was implemented in [287–289]. Vanishing boundary conditions are assumed at the boundaries of the domain (a cylinder of radius R and length 2L).

Such lattice-based techniques achieve a very high numerical precision regardless of the underlying nuclear geometry. Their convergence is essentially characterized by the order M of the B-spline and the maximum value ${{ \Omega }_{\text{max}}}$ of the z-projection of angular momentum (For Hartree–Fock calculations, ${{ \Omega }_{\text{max}}}$ is simply equal to the maximum j-value of occupied single particle orbitals). They are especially suited to dealing with very deformed nuclei, weakly-bound systems or nuclei at high temperature. Comparisons with HO basis expansion techniques given in [114, 115, 161, 289] show that a prohibitively large number of basis states would be needed to reach similar precision. The downside of the BSCM techniques is the larger computational cost: explicit parallelization across a few dozen cores is needed to keep run time within a few hours. In addition, extending the codes to handle either finite-range nuclear potentials or fully triaxial geometries would be costly, and such extensions do not exist yet.

A slightly different lattice technique relies on using a variational method with a set of Lagrange functions introduced in [290]. As in the BSCM, both functions and operators are represented on the resulting mesh, the functions by a vector and the operators by a matrix. This technique is implemented in three-dimensional Cartesian meshes in the code ev of [291, 292]. The use of Lagrange meshes to estimate derivatives delivers high and controllable numerical precision as demonstrated in [293]. In particular, it is possible to limit numerical errors to only a few dozen keV across an arbitrarily large deformation span with even relatively coarse meshes. Note that the implementation of finite-range potentials represents a significant increase in computations, see however [294] for a recent example using the finite range Gogny force.

Another popular technique used in DFT solvers based on the coordinate space representation consists in evaluating derivatives, which are needed to define the Laplacian in the kinetic energy, in momentum space. Consider a Cartesian grid of ${{N}_{x}}\times {{N}_{y}}\times {{N}_{z}}$ points defined in such a way that ${{x}_{i}}=\text{id}x$ for $i=1,\ldots {{N}_{x}}$ (and similarly for the coordinates y and z). Then the momentum is discretized according to $\left. {{p}_{i}}/\hslash =\text{i}2\pi /L\right)$. As is well known, derivatives in momentum space are simple multiplications. The transformation back to coordinate space can be performed by Fast Fourier Transforms. This technique has been used in [13, 14, 33, 34, 79, 182, 295, 296].

We finish this section by mentioning two slightly alternative lattice-based techniques. Finite element analysis was introduced in the series of papers [297–300] to solve the equations of the relativistic mean-field. Although this technique seemed promising, it was not disseminated further. Most recently, the NUCLEI SciDAC collaboration has published in [301] a new DFT solver based on multi-resolution analysis and wavelet expansions. The advantage of multi-resolution is the possibility to impose the numerical precision desired.

The HFB equation (21) is a non-linear problem since the HFB matrix $\mathcal{H}$ depends on the generalized density $\mathcal{R}$—recall (37) and (23). Even in the BCS approximation, the non-linearity of the equation remains since the HF mean field h is still a functional of the one-body density ρ. There are three main strategies to solve such problems: successive diagonalizations, gradient methods or imaginary time evolution. In each case, the algorithm must be initialized. This can be done by solving, e.g. the Schrödinger equation with a Nilsson potential followed by the BCS approximation: this provides an initial one-body density matrix ${{\rho}^{(0)}}$ and pairing tensor ${{\kappa}^{(0)}}$, which define the HF or HFB matrix at the first iteration. Given this initialization, the three algorithms mentioned proceed as follows:

• Successive iterations: In this method the HFB equation is written in the form of a diagonalization problem
  $$\begin{eqnarray}&&\mathcal{H}W=W\mathcal{E}\end{eqnarray} \tag{ 140 }$$
  with W having the structure of (11) and $\mathcal{E}$ of (25). The W⁽ⁱ⁾ matrix at iteration i is used to compute the next iteration of the HFB matrix, ${{\mathcal{H}}^{(i+1)}}$, which is diagonalized to obtain W⁽ⁱ⁺¹⁾. The iterative process is repeated until convergence, i.e. $\parallel {{W}^{(i)}}-{{W}^{(i+1)}}\parallel \leqslant \varepsilon$ within a given accuracy . The convergence of the method is not guaranteed and 'jumping' solutions such that ${{W}^{(i)}}={{W}^{(i+q)}}$ can occur. These deficiencies are usually overcome by annealing the density at iteration i +1 with the density at iteration i
  $$\begin{eqnarray}&&{{\mathcal{R}}^{(i+1)}}\to (1-\alpha ){{\mathcal{R}}^{(i+1)}}+\alpha {{\mathcal{R}}^{(i)}}\end{eqnarray}$$
  with the annealing parameter α. More advanced linear mixing schemes have been introduced in quantum chemistry and ported to nuclear structure in [302]. Constraints are handled very efficiently by using methods based on the augmented Lagrangian method or approximate linear response theory, see [162, 186, 192]. In the special case of the HF + BCS approximation to the HFB equation, the same overall algorithm is applied to diagonalize the HF mean field h instead of the HFB matrix $\mathcal{H}$.
• Gradient methods: The HFB equation is a direct consequence of the variational principle on the HFB energy. Therefore solving the HFB equation is equivalent to finding a minimum of the HFB energy. The gradient method was used as early as the late nineteen seventies in [303, 304] to solve the HFB equations using a phenomenological approach to determine the gradient step and is based on the representation of the HFB equation in terms of the Thouless matrix Z, see p 10. The efficient handling of constraints inherent to the method is especially convenient in fission calculations: during the iterative process one has to ensure that the descendant direction is orthogonal to the gradient of the constraint condition. The orthogonality is imposed by modifying the objective function by subtracting the mean value of the constraint operator multiplied by a Lagrange multiplier that is subsequently used to impose the orthogonality condition at every iteration. The procedure is straightforwardly extended to many constraints. The iteration count of the method is high, although the cost of evaluating the gradient is just slightly higher than evaluating the $\mathcal{H}$ matrix. To reduce the number of iterations, the conjugate gradient method introduced in [305] chooses the steepest descent direction as the 'conjugate' of the previous one. The method is very efficient for a quadratic function but the cost of each iteration increases because of the need to do a line minimization at each step. The Newton or second order method relies on the use of the Hessian matrix to modify the gradient direction in such a way that the minimum is reached in just one iteration for a quadratic function. The iteration count is severely reduced but the cost of evaluating the Hessian is very high (in the HFB case, the Hessian resembles the matrix of linear response theory), as is the evaluation of its inverse. For not so many degrees of freedom the method is competitive as shown in [306] but becomes prohibitively expensive in typical applications in nuclear physics with effective forces. An enormous simplification of the method is achieved if the Hessian matrix is assumed to be diagonal dominant as the most important contribution to the diagonal matrix elements is the sum of two quasiparticle energies. In this way the computation and inversion of the Hessian matrix is enormously simplified. This method was implemented in the hfbaxial and hfbtri DFT solvers for the Gogny force. For an early account of the different gradient method techniques applied to the solution of the HF and HFB problems with an emphasis on the election of the gradient step, see [307]
• Imaginary time method: In the particular case of the HF + BCS approximation, one can also use the imaginary time method. It is based on the general result that the eigenvalues of a Hamiltonian $\hat{H}$ evolve with time through an oscillatory phase factor $\exp \left(-\text{i}/\hslash {{E}_{n}}t\right)$ depending on the eigenvalue E_n. If time is replaced by a complex quantity $t\to -\text{i}\tau$ the complex phase in front of each eigenvalue becomes real and a decaying function with a decay rate that increases with excitation energy. Therefore, applying the time evolution operator with imaginary time to a general linear combination of eigenvalues of the Hamiltonian will converge to the lowest eigenvalue (ground-state energy) after a sufficiently long time. The problem with this method is that the time evolution operator involves the exponential of the true Hamiltonian, typically approximated by the two-body effective Hamiltonian (27), which is very expensive to compute. In practice, the full Hamiltonian $\hat{H}$ is therefore replaced by a simpler approximation: In nuclear physics, this method was introduced in [308] by replacing $\hat{H}$ by the HF Hamiltonian $\hat{h}$. This is the technique implemented in the ev suite of code of [291, 292] to solve the HF + BCS equation. The direct extension of this algorithm to the case of the HFB equation is not trivial because the HFB matrix is unbounded from below: the lowest eigenvalue is infinite resulting in the divergence of the algorithm at $\tau \to +\infty$. In [309], the authors introduced the two-basis method, where the HFB matrix is diagonalized in a small discrete basis composed of the eigenstates of the HF Hamiltonian obtained from the imaginary time evolution.

As emphasized in section 3.3, time-dependent density functional theory (TDDFT) provides, at least on paper, a convenient framework to simulate fission in real time, since it does not require introducing collective variables or scission configurations. Assuming the original Runge–Gross theorem could be extended to the nuclear case, the (unknown) energy functional that gives the exact ground-state energy would only depend on the local, time-dependent, one-body density, and the Kohn–Sham scheme would reduce to solving Skyrme-like TDHF equations.

All existing implementations of TDHF in computer codes are based on the coordinate-space formulation of the HF equations in a box, using both absorption layers on the edges of the box and vanishing boundary conditions. The coordinate space approach is necessary because TDHF solutions at large time can have very extended geometries that cannot be accurately represented by a one-centre basis expansion such as the familiar HO basis. Because of the computational cost of the coordinate space approach, the first implementations of TDHF used several simplifications such as axial symmetry as in [251, 310]. The first application of a fully three-dimensional, Skyrme TDHF calculation was published in [311] in 1997 and was based on the adaptation of the ev8 HF +BCS solver. Very recently, a full 3D implementation of the TDHF equations in coordinate space using the Fourier representation of spatial derivatives has been published [296].

The simplest way to include pairing correlations in TDHF is the BCS approximation. Such a TDHF +BCS solver was developed in [256] and applied to the study of fission of ²⁵⁸Fm in [170]. The solver is also an extension of the ev HF +BCS computer program and relies on the same underlying technology, particularly the use of a 3D Cartesian mesh and a set of Lagrange functions—see the two previous sections 4.1.2 and 4.1.3.

The full TDHFB equation is substantially more involved than the TDHF or TDHF +BCS equation because of the non-zero occupation of high-lying, de-localized quasiparticle states. The full TDHFB equation in nuclei was originally solved in spherical symmetry in [252]. For arbitrary deformations of the nuclear shape, the TDHFB equation has been implemented in the canonical basis, where the density matrix is diagonal and the pairing tensor has the canonical form, in [253–255]. The most advanced implementation of TDHFB today is by Bulgac, Roche and collaborators and was described in [312]. Their code, which is massively parallel and has been ported to GPU architectures, implements a local Kohn–Sham scheme for TDDFT dubbed the time-dependent superfluid local density approximation (equivalent to a TDHFB theory with a functional of the local density $\rho \left(\boldsymbol{r}\right)$ only) and was very recently applied for the first time to the description of neutron-induced fission in [182]. The code uses fast discrete Fourier transforms to evaluate derivatives on a three-dimensional Cartesian lattice, and a multi-step predictor-modifier-corrector method for the time evolution.

In contrast to real-time dynamics described by TDDFT methods, the description of fission dynamics as a large-amplitude collective motion driven by a small set of collective variables is especially suited to the calculation of the distributions of fission fragment properties. The time-dependent generator coordinate method (TDGCM) developed in the 1980s at CEA Bruyères-le-Châtel is currently the only microscopic theory capable of producing realistic fission fragment distributions. Until now, it has only been applied under the Gaussian overlap approximation. The very first applications of this approach to the dynamics of fission were reported in a series of papers by Berger and Gogny in [120, 123, 257]. The first actual calculations of fission fragment mass distributions with this technique were published in [258, 284], with additional results reported in [171]. The implementation of the TDGCM was based on solving the collective Schrödinger equation by using finite differences for derivatives. This choice imposed the use of a regular grid of points that did not offer the possibility to discard regions of the collective space that were irrelevant to the dynamics (for example the points with a very high energy). In practice, applications were restricted to two-dimensional collective spaces, and large computational resources were needed to achieve good numerical accuracy. Very recently, the collaboration between CEA and LLNL developed a new program to solve the TDGCM + GOA equations for an arbitrary number of collective variables using finite element analysis [313].

Computing collective inertia is usually accomplished by using the perturbative cranking formula (110), which only requires moments of the collective operators in quasiparticle space and two-quasiparticle energies. However, as recalled in section 3.1.2 (see also figure 6 p 16), this approximation falls too short and can underestimate the inertia by more than 40%. On the other hand, the exact evaluation of the ATDHFB inertia would require both inverting the linear response matrix $\mathcal{M}$ of (88) and evaluating the partial derivatives of the generalized density matrix with respect to the collective variables, see (109) p 27. There exist analytical formulas giving the partial derivatives (86) as a function of the matrix of the collective variables, but this again requires inverting the linear response matrix as seen from (87). A more convenient and thrifty procedure consists in the numerical evaluation of the derivatives using finite difference formulas such as the popular and accurate centred difference formula. It requires the densities $\mathcal{R}\left(\boldsymbol{q}\pm \delta \boldsymbol{q}\right)$, which can be obtained by any DFT solver by slightly modifying the value of the constraints $\boldsymbol{q}\to \boldsymbol{q}\pm \delta \boldsymbol{q}$. The choice of $\delta \boldsymbol{q}$ obviously depends on the collective variable and usually also on particle number (as $\boldsymbol{q}$ usually has such dependence). This approach has been taken in [216].

Concerning the inversion of the linear response matrix, a possibility used by the authors of [314] is to build the linear response matrix in a reduced subspace of two-quasiparticle excitations and compute the inverse there. This is a computationally intensive task, but the main advantage is that the convergence of the procedure can be easily tested by slightly increasing the size of the subspace and repeating the calculation. Other authors use the fact that the linear response matrix is diagonal dominant, with two-quasiparticle energies on the diagonal, to build an iterative method to compute the action of the inverse on a given vector (that is the only quantity required to compute the inertias). This method has only been used to compute the inertias associated to Goldstone modes in [315]. Yet another approach aiming at the evaluation of the inertias associated to Goldstone modes (Thouless–Valatin moments of inertia) has been recently proposed in [316]. It resorts to the ideas of the finite amplitude method (FAM) to compute the action of the inverse of the linear response matrix on a given vector.