Generalized time-dependent generator coordinate method for induced fission dynamics

Abstract

The generalized time-dependent generator coordinate method (TD-GCM) is extended to include pairing correlations. The correlated GCM nuclear wave function is expressed in terms of time-dependent generator states and weight functions. The particle−hole channel of the effective interaction is determined by a Hamiltonian derived from an energy density functional, while pairing is treated dynamically in the standard BCS approximation with time-dependent pairing tensor and single-particle occupation probabilities. With the inclusion of pairing correlations, various time-dependent phenomena in open-shell nuclei can be described more realistically. The model is applied to the description of saddle-to-scission dynamics of induced fission. The generalized TD-GCM charge yields and total kinetic energy distribution for the fission of ²⁴⁰Pu, are compared to those obtained using the standard time-dependent density functional theory (TD-DFT) approach, and with available data.

Graphical abstract

Keywords

nuclear density functional theory / generator coordinate method / fission dynamics

1 Introduction

The process of spontaneous or induced fission, by which an atomic nucleus breaks into two or more fragments, with the corresponding release of an enormous amount of energy, has become one of the best studied phenomena in nuclear physics. From a modern perspective, the fission process is particularly important as a representative case of large-amplitude collective motion in a self-bound mesoscopic system, that exhibits both classical and quantal characteristics. Fission is also relevant for the stability of superheavy elements, production of short-lived exotic nuclides far from stability, nuclear astrophysics and the mechanism of nucleosynthesis. A wealth of experimental results on nuclear fission have been accumulated over more than eight decades, and a basic understanding of the mechanism gained. Several very successful phenomenological methods have been developed that reproduce, to various degrees of accuracy, low- and medium-energy fission observables. A unified microscopic framework for the description of the entire fission process, from the quasi-stationary initial state to the outer fission barrier (saddle point), and the complicated dissipative dynamics that lead to the emergence of excited fragments at scission, remains a formidable challenge for nuclear theory.

In the last decade significant advances have been made in the development of microscopic models of induced fission dynamics. The two basic approaches are based on the time-dependent generator coordinate method (TD-GCM) and the time-dependent density functional theory (TD-DFT), respectively. In the former, the nuclear wave function is represented by a superposition of generator states that are functions of collective coordinates, and can be applied to an adiabatic description of the entire fission process [1–7]. This approach is fully quantum mechanical but only takes into account collective degrees of freedom and, therefore, cannot be used to describe the highly dissipative dynamics that occurs beyond the saddle point. While the TD-DFT framework consistently incorporates the one-body dissipation mechanism, it can only model a single fission event by propagating nucleons independently toward scission and beyond [8–16]. In the fission case, TD-DFT models the classical evolution of independent nucleons in mean-field potentials, and cannot be applied in the classically forbidden region of the collective space nor does it take into account quantum fluctuations.

In Ref. [17], an implementation of the generalized TD-GCM has been developed and applied to the dynamics of small and large amplitude collective motion of atomic nuclei. Both the generator states and weight functions of the GCM correlated wave function depend on time. The initial generator states are obtained as solution of deformation-constrained self-consistent mean-field equations, and are evolved in time by the standard mean-field equations of nuclear density functional theory (TD-DFT). The TD-DFT trajectories are used as a basis in which the TD-GCM wave function is expanded. As a first application, in Ref. [17] the generalized TD-GCM has been employed in a study of excitation energies and spreading width of giant resonances, and in an initial exploration of the dynamics of induced fission. An equivalent realization of the model has very recently been applied to a study of collective multiphonon states in nuclei [18].

In this work, the generalized TD-GCM is extended with the inclusion of dynamical pairing correlations and employed in a systematic investigation of induced fission dynamics. Pairing is included in the TD-DFT framework in a standard BCS approximation, with time-dependent pairing tensor and occupation probabilities. A basis of TD-DFT fission trajectories, generally non-orthogonal and overcomplete, is used to build the correlated TD-GCM wave function. With the coherent superposition of TD-DFT trajectories in the generalized TD-GCM, fission dynamics is described fully quantum mechanically in an approach that extends beyond the adiabatic approximation of the standard GCM and, simultaneously, includes quantum fluctuations.

In Section 2, we outline the theoretical framework of generalized TD-GCM that also takes into account pairing correlations. Section 3 collects the numerical details of the fission calculation. As an application of the generalized TD-GCM, the dynamics of induced fission of ²⁴⁰Pu is discussed in Section 4. Section 5 summarizes the results and presents a brief outlook for future studies.

2 Theoretical framework: Generalized time-dependent GCM with pairing interactions

The Griffin−Hill−Wheeler (GHW) ansatz for the TD-GCM correlated nuclear wave function reads [5, 19, 20]

(1)$|\mathrm{\Psi }(t)⟩={\displaystyle {\int }_{\mathit{q}}\mathrm{d}\mathit{q}{f}_{\mathit{q}}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩,}$

where the vector $\mathit{q}$ denotes the continuous real generator coordinates that parametrize the collective degrees of freedom. This wave function is a linear superposition of, generally non-orthogonal, many-body generator states $|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩$, and $f_{\mathit{q}}(t)$ are the corresponding complex-valued weight functions. The generalized TD-GCM without the inclusion of pairing correlations has been implemented in the first part of this work [17]. In this study, pairing is also taken into account, and the discretized generator coordinates are the mass multipole moments (axial quadrupole and octupole) of the nucleon density distribution. Thus, the nuclear wave function

(2)$|\mathrm{\Psi }(t)⟩={\displaystyle \sum _{\mathit{q}}{f}_{\mathit{q}}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩,}$

is the solution of the time-dependent equation

(3)$\mathrm{i}\hslash {\mathrm{\partial }}_t|\mathrm{\Psi }(t)⟩=\hat{H}|\mathrm{\Psi }(t)⟩,$

where $\hat{H}$ is the Hamiltonian of the nuclear system. From a time-dependent variational principle [20], one obtains the equation of motion for the weight functions

(4)$\mathrm{i}\hslash \mathcal{N}\dot{f}=(\mathcal{H}-{\mathcal{H}}^{MF})f,$

which, in the discretized collective space, reads

(5)${\displaystyle \sum _{\mathit{q}}\mathrm{i}\hslash {\mathcal{N}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}(t){\mathrm{\partial }}_t{f}_{\mathit{q}}(t)+{\displaystyle \sum _q{\mathcal{H}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}^{MF}(t)f_{\mathit{q}}(t)={\displaystyle \sum _q{\mathcal{H}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}(t)f_{\mathit{q}}(t).}}}$

The time-dependent kernels

(6a)${\mathcal{N}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}(t)=⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩,$

(6b)${\mathcal{H}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}(t)=⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|\hat{H}|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩,$

(6c)${\mathcal{H}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}^{MF}(t)=⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|\mathrm{i}\hslash {\mathrm{\partial }}_t|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩$

include the overlap, the Hamiltonian, and the time derivative of the generator states, respectively.

2.1 Time-dependent quasiparticle vacuum $|{\mathrm{\Phi }}_{\mathbf{q}}(t)⟩$

The evolution in time of the quasiparticle vacuum characterized by a vector of generator coordinates $\mathit{q}$

(7)$|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩={\displaystyle \prod _{k>0}[{\mu }_{\mathit{q},k}(t)+{\nu }_{\mathit{q},k}(t)c_{\mathit{q},k}^{†}(t)c_{\mathit{q},\overline{k}}^{†}(t)]|-⟩}$

is modeled by the time-dependent covariant density functional theory [16, 21], using the time-dependent BCS approximation [22, 23]. In Eq. (7), ${\mu }_{\mathit{q},k}(t)$ and ${\nu }_{\mathit{q},k}(t)$ are the parameters of the transformation between the canonical and quasiparticle bases, and $c_{\mathit{q},k}^{†}(t)$ denotes the creation operator associated with the canonical state ${\varphi }_k^{\mathit{q}}(\mathit{r},t)$. The time evolution of ${\varphi }_k^{\mathit{q}}(\mathit{r},t)$ is determined by the time-dependent Dirac equation

(8)$\mathrm{i}{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}}{\varphi }_k^{\mathit{q}}(\mathit{r},t)=[{\hat{h}}^{\mathit{q}}(\mathit{r},t)-{\epsilon }_k^{\mathit{q}}(t)]{\varphi }_k^{\mathit{q}}(\mathit{r},t),$

where ${\epsilon }_k^{\mathit{q}}(t)=⟨{\psi }_k^{\mathit{q}}|{\hat{h}}^{\mathit{q}}|{\psi }_k^{\mathit{q}}⟩$ is the single-particle energy, and the single-particle Hamiltonian ${\hat{h}}^{\mathit{q}}(\mathit{r},t)$ reads

(9)${\hat{h}}^{\mathit{q}}(\mathit{r},t)=\mathit{\alpha }\cdot (\hat{\mathit{p}}-{\mathit{V}}_{\mathit{q}})+V_{\mathit{q}}^0+\beta (m_N+S_{\mathit{q}}).$

Here, $m_N$ is the nucleon mass, and the scalar $S_{\mathit{q}}(\mathit{r},t)$ and four-vector $V_{\mathit{q}}^{\mu }(\mathit{r},t)$ potentials at every instant are determined by the time-dependent densities and currents in the isoscalar−scalar, isoscalar−vector, and isovector−vector channels. In the present study, we employ the point-coupling relativistic energy density functional PC-PK1 [24], and the explicit expressions for the potentials read

(10a)$S_{\mathit{q}}(\mathit{r})={\alpha }_S{\rho }_S^{\mathit{q}}+{\beta }_S({\rho }_S^{\mathit{q}}{)}^2+{\gamma }_S({\rho }_S^{\mathit{q}}{)}^3+{\delta }_S\mathrm{\Delta }{\rho }_S^{\mathit{q}},$

(10b)$\begin{array}{rl}{V}_{\mathit{q}}^{\mu }(\mathit{r})=& {\alpha }_V{j}^{\mathit{q},\mu }+{\gamma }_V(j^{\mathit{q},\mu }{j}_{\mu }^{\mathit{q}})j^{\mathit{q},\mu }+{\delta }_V\mathrm{\Delta }{j}^{\mathit{q},\mu }+{\tau }_3{\alpha }_{TV}{j}_{TV}^{\mathit{q},\mu }\\ & +{\tau }_3{\delta }_{TV}\mathrm{\Delta }{j}_{TV}^{\mathit{q},\mu }+e\frac{1-{\tau }_3}{2}{A}^{\mathit{q},\mu },\end{array}$

where ${\tau }_3$ is the isospin Pauli matrix, and $A^{\mathit{q},\mu }$ is the electromagnetic vector potential. The densities and currents are defined in terms of occupied single-particle wave functions ${\varphi }_k^{\mathit{q}}(\mathit{r},t)$ and the occupation probabilities $n_k^{\mathit{q}}(t)$:

(11a)${\rho }_S^{\mathit{q}}(\mathit{r},t)={\displaystyle \sum _k^{l_q}{n}_{\mathit{q},k}(t){\overline{\varphi }}_k^{\mathit{q}}(\mathit{r},t){\varphi }_k^{\mathit{q}}(\mathit{r},t),}$

(11b)$j^{\mathit{q},\mu }(\mathit{r},t)={\displaystyle \sum _k^{l_q}{n}_{\mathit{q},k}(t){\overline{\varphi }}_k^{\mathit{q}}(\mathit{r},t){\gamma }^{\mu }{\varphi }_k^{\mathit{q}}(\mathit{r},t),}$

(11c)$j_{TV}^{\mathit{q},\mu }(\mathit{r},t)={\displaystyle \sum _k^{l_q}{n}_{\mathit{q},k}(t){\overline{\varphi }}_k^{\mathit{q}}(\mathit{r},t){\gamma }^{\mu }{\tau }_3{\varphi }_k^{\mathit{q}}(\mathit{r},t),}$

where $l_q$ is the number of the canonical basis states. The time evolution of the occupation probability $n_{\mathit{q},k}(t)=|{\nu }_{\mathit{q},k}(t){|}^2$ and pairing tensor ${\kappa }_{\mathit{q},k}(t)={\mu }_{\mathit{q},k}^{\ast }(t){\nu }_{\mathit{q},k}(t)$ is governed by the following equations:

(12a)$\mathrm{i}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{n}_{\mathit{q},k}(t)={\kappa }_{\mathit{q},k}(t){\mathrm{\Delta }}_{\mathit{q},k}^{\ast }(t)-{\kappa }_{\mathit{q},k}^{\ast }(t){\mathrm{\Delta }}_{\mathit{q},k}(t),$

(12b)$\mathrm{i}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\kappa }_{\mathit{q},k}(t)=[{\epsilon }_k^{\mathit{q}}(t)+{\epsilon }_{\overline{k}}^{\mathit{q}}(t)]{\kappa }_{\mathit{q},k}(t)+{\mathrm{\Delta }}_{\mathit{q},k}(t)[2n_{\mathit{q},k}(t)-1]$

(for details, see Refs. [22, 23]). In time-dependent calculations, a monopole pairing interaction is employed, and the gap parameter ${\mathrm{\Delta }}_{\mathit{q},k}(t)$ is defined in terms of single-particle energies and the pairing tensor,

(13)${\mathrm{\Delta }}_{\mathit{q},k}(t)=\left[G{\displaystyle \sum _{k^{\prime }>0}f({\epsilon }_{k^{\prime }}^{\mathit{q}}){\kappa }_{\mathit{q},k^{\prime }}}\right]f({\epsilon }_k^{\mathit{q}}),$

where $f({\epsilon }_k^{\mathit{q}})$ is the cut-off function for the pairing window [23], and $G$ is the pairing strength.

2.2 Overlap kernel ${\mathcal{N}}_{{\mathbf{q}}^{\prime }\mathbf{q}}(t)$

According to Eq. (7), the expression for the overlap kernel Eq. (6a) can be written in the following form:

(14)$\begin{array}{rl}{\mathcal{N}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}(t)=& ⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩\\ =& ⟨-|\prod _{k^{\prime }>0}[{\mu }_{{\mathit{q}}^{\prime },k^{\prime }}^{\ast }(t)+{\nu }_{{\mathit{q}}^{\prime },k^{\prime }}^{\ast }(t)c_{{\mathit{q}}^{\prime },{\overline{k}}^{\prime }}(t)c_{{\mathit{q}}^{\prime },k^{\prime }}(t)]\\ & \cdot \prod _{k>0}[{\mu }_{\mathit{q},k}(t)+{\nu }_{\mathit{q},k}(t)c_{\mathit{q},k}^{†}(t)c_{\mathit{q},\overline{k}}^{†}(t)]|-⟩\\ =& \frac{(-1{)}^{(l_{{\mathit{q}}^{\prime }}-1)l_{{\mathit{q}}^{\prime }}/2}}{{{\displaystyle \prod }}_{k^{\prime }}^{l_{{\mathit{q}}^{\prime }}/2}{{\displaystyle \prod }}_k^{l_{\mathit{q}}/2}{\nu }_{{\mathit{q}}^{\prime },k^{\prime }}^{\ast }{\nu }_{\mathit{q},k}}⟨-|{\beta }_{{\mathit{q}}^{\prime },1}^{†}...{\beta }_{{\mathit{q}}^{\prime },l_{{\mathit{q}}^{\prime }}}^{†}{\beta }_{\mathit{q},1}...{\beta }_{\mathit{q},l_{\mathit{q}}}|-⟩,\end{array}$

where ${\beta }_{\mathit{q},k}^{†}$ is the quasi-particle creation operator associated with the quasiparticle vacuum $|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩$. The overlap between two quasi-particle vacuums can be calculated using the Pfaffian algorithms developed in Refs. [25, 26].

2.3 Energy kernel ${\mathcal{H}}_{{\mathbf{q}}^{\prime }\mathbf{q}}(t)$

For the point-coupling relativistic energy density functional PC-PK1 [24], one obtains the expression for the energy kernel ${\mathcal{H}}^{\mathrm{D}\mathrm{F}}(t)$, under the assumption [9] that it only depends on the transition densities at time $t$:

(15)$\begin{array}{rl}{\mathcal{H}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}^{\mathrm{D}\mathrm{F}}(t)=& ⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|{\hat{H}}^{\mathrm{D}\mathrm{F}}|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩=⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩\cdot \int {\mathrm{d}}^{3}r\{{\rho }_{\mathrm{k}\mathrm{i}\mathrm{n}}(\mathit{r},t)\\ & +\frac{{\alpha }_S}{2}{\rho }_S(\mathit{r},t{)}^2+\frac{{\beta }_S}{3}{\rho }_S(\mathit{r},t{)}^3\\ & +\frac{{\gamma }_S}{4}{\rho }_S(\mathit{r},t{)}^4+\frac{{\delta }_S}{2}{\rho }_S(\mathit{r},t)\mathrm{\Delta }{\rho }_S(\mathit{r},t)\\ & +\frac{{\alpha }_V}{2}{j}^{\mu }(\mathit{r},t)j_{\mu }(\mathit{r},t)+\frac{{\gamma }_V}{4}[j^{\mu }(\mathit{r},t)j_{\mu }(\mathit{r},t){]}^2\\ & +\frac{{\delta }_V}{2}{j}^{\mu }(\mathit{r},t)\mathrm{\Delta }{j}_{\mu }(\mathit{r},t)+\frac{{\alpha }_{TV}}{2}{j}_{TV}^{\mu }(\mathit{r},t)\cdot [j_{TV}(\mathit{r},t){]}_{\mu }\\ & +\frac{{\delta }_{TV}}{2}{j}_{TV}^{\mu }(\mathit{r},t)\cdot \mathrm{\Delta }[j_{TV}(\mathit{r},t){]}_{\mu }+\frac{e^2}{2}{j}_p^{\mu }(\mathit{r},t)A_{\mu }(\mathit{r},t)\},\end{array}$

where the densities and currents ${\rho }_{\mathrm{k}\mathrm{i}\mathrm{n}}$, ${\rho }_S$, $j^{\mu }$, $j_{TV}^{\mu }$, and $j_p^{\mu }$ read

(16a)${\rho }_{\mathrm{k}\mathrm{i}\mathrm{n}}(\mathit{r},t)={\displaystyle \sum _{k^{\prime }}^{l_{{\mathit{q}}^{\prime }}}{\displaystyle \sum _k^{l_{\mathit{q}}}{\overline{\varphi }}_{k^{\prime }}^{{\mathit{q}}^{\mathbf{\prime }}}(\mathit{r},t)(-\mathrm{i}\mathit{\gamma }\cdot \mathbf{\nabla }+m_N){\varphi }_k^{\mathit{q}}(\mathit{r},t){\rho }_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q},k^{\prime }k}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}(t),}}$

(16b)${\rho }_S(\mathit{r},t)={\displaystyle \sum _{k^{\prime }}^{l_{{\mathit{q}}^{\prime }}}{\displaystyle \sum _k^{l_{\mathit{q}}}{\overline{\varphi }}_{k^{\prime }}^{{\mathit{q}}^{\mathbf{\prime }}}(\mathit{r},t){\varphi }_k^{\mathit{q}}(\mathit{r},t){\rho }_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q},k^{\prime }k}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}(t),}}$

(16c)$j^{\mu }(\mathit{r},t)={\displaystyle \sum _{k^{\prime }}^{l_{{\mathit{q}}^{\prime }}}{\displaystyle \sum _k^{l_{\mathit{q}}}{\overline{\varphi }}_{k^{\prime }}^{{\mathit{q}}^{\mathbf{\prime }}}(\mathit{r},t){\gamma }^{\mu }{\varphi }_k^{\mathit{q}}(\mathit{r},t){\rho }_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q},k^{\prime }k}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}(t),}}$

(16d)$j_{TV}^{\mu }(\mathit{r},t)={\displaystyle \sum _{k^{\prime }}^{l_{{\mathit{q}}^{\prime }}}{\displaystyle \sum _k^{l_{\mathit{q}}}{\overline{\varphi }}_{k^{\prime }}^{{\mathit{q}}^{\prime }}(\mathit{r},t){\tau }_3{\gamma }^{\mu }{\varphi }_k^{\mathit{q}}(\mathit{r},t){\rho }_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q},k^{\prime }k}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}(t),}}$

(16e)$j_p^{\mu }(\mathit{r},t)={\displaystyle \frac{1-{\tau }_3}{2}}{\displaystyle \sum _{k^{\prime }}^{l_{{\mathit{q}}^{\prime }}}{\displaystyle \sum _k^{l_{\mathit{q}}}{\overline{\varphi }}_{k^{\prime }}^{{\mathit{q}}^{\mathbf{\prime }}}(\mathit{r},t){\gamma }^{\mu }{\varphi }_k^{\mathit{q}}(\mathit{r},t){\rho }_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q},k^{\prime }k}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}(t).}}$

The transition density matrix ${\rho }^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}(t)$ is defined by the following relation

(17)$\begin{array}{rl}{\rho }_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q},k^{\prime }k}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}(t)=& \frac{⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|c_{{\mathit{q}}^{\mathbf{\prime }},k^{\prime }}^{†}(t)c_{\mathit{q},k}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩}{⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩}\\ =& {\nu }_{{\mathit{q}}^{\prime },k^{\prime }}^{\ast }(t){\nu }_{\mathit{q},k}(t)\frac{⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|{\beta }_{{\mathit{q}}^{\mathbf{\prime }},{\overline{k}}^{\prime }}(t){\beta }_{\mathit{q},\overline{k}}^{†}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩}{⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩}.\end{array}$

The numerator of the transition density matrix ${\rho }_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q},k^{\prime }k}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}(t)$ is the overlap between a quasi-particle vacuum with $\left(l_{{\mathit{q}}^{\prime }}-1\right)$ quasi-particle levels and a quasi-particle vacuum with $\left(l_{\mathit{q}}-1\right)$ quasi-particle levels. It can be calculated using the Pfaffian algorithms [25, 26].

For the monopole pairing interaction, the pairing Hamiltonian operator ${\hat{H}}^{\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r}}$ in 3D-lattice space is defined as

(18)${\hat{H}}_{\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r}}=-{\displaystyle \sum _{r_1,s_1>0,r_2,s_2>0}G(c_{r_1,s_1}^{†}{c}_{r_1,{\overline{s}}_1}^{†})(c_{r_2,{\overline{s}}_2}{c}_{r_2,s_2}),}$

where $c_{r_1,s_1}^{†}$ is the creation operator for the lattice coordinate wave function $|r_1,s_1⟩$, and $r_1$ is the index of the lattice point and $s_1$ is the index of the spin. One obtains the expression for the pairing part of energy kernel ${\mathcal{H}}^{\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r}}(t)$ in the 3D-lattice space

(19)$\begin{array}{rl}{\mathcal{H}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}^{\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r}}(t)=& ⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|{\hat{H}}^{\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r}}|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩\\ =& -G⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩\sum _{k_1,k_2,k_3,k_4>0}[f({\epsilon }_{k_1}^{{\mathit{q}}^{\prime }})f({\epsilon }_{k_2}^{\mathit{q}})f({\epsilon }_{k_3}^{{\mathit{q}}^{\prime }})f({\epsilon }_{k_4}^{\mathit{q}}){]}^{1/2}\\ & \times ⟨{\varphi }_{k_1}^{{\mathit{q}}^{\prime }}|{\varphi }_{{\overline{k}}_2}^{\mathit{q}}⟩⟨{\varphi }_{{\overline{k}}_3}^{{\mathit{q}}^{\prime }}|{\varphi }_{k_4}^{\mathit{q}}⟩[{\kappa }_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q},k_2{k}_1}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}(t){]}^{\ast }{\kappa }_{\mathit{q}{\mathit{q}}^{\mathbf{\prime }},k_3{k}_4}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}(t),\end{array}$

where the transition pairing tensor matrix ${\kappa }^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}(t)$ is defined by the following relation:

(20)$\begin{array}{rl}{\kappa }_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q},k^{\prime }k}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}(t)=& \frac{⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|c_{{\mathit{q}}^{\mathbf{\prime }},k^{\prime }}(t)c_{\mathit{q},k}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩}{⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩}\\ =& {\mu }_{{\mathit{q}}^{\prime },k^{\prime }}^{\ast }(t){\nu }_{\mathit{q},k}(t)\frac{⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|{\beta }_{{\mathit{q}}^{\mathbf{\prime }},k^{\prime }}(t){\beta }_{\mathit{q},\overline{k}}^{†}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩}{⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩}.\end{array}$

In the BCS model, $|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩$ is not an eigenstate of the neutron (proton) number operator $\hat{N}$ ($\hat{Z}$), and its expectation value in the collective wave function generally deviates from the desired neutron number $N_0$ (proton number $Z_0$). The method developed in Ref. [27] is used to correct for variations of the nucleon number. The energy kernel finally reads

(21)$\begin{array}{rl}{\mathcal{H}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}(t)=& {\mathcal{H}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}^{\mathrm{D}\mathrm{F}}(t)+{\mathcal{H}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}^{\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r}}(t)-{\lambda }_N^{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}(t)[⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|\hat{N}|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩-N_0]\\ & -{\lambda }_Z^{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}(t)[⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\mathbf{\prime }}}(t)|\hat{Z}|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩-Z_0],\end{array}$

where ${\lambda }_i^{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}(t)$ is defined as the average of the chemical potentials ${\lambda }_i^{{\mathit{q}}^{\mathbf{\prime }}}(t)$ and ${\lambda }_i^{\mathit{q}}(t)$,

(22)${\lambda }_i^{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}(t)={\displaystyle \frac{{\lambda }_i^{{\mathit{q}}^{\mathbf{\prime }}}(t)+{\lambda }_i^{\mathit{q}}(t)}{2}},i=N,Z.$

2.4 Mean-field kernel ${\mathcal{H}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}^{MF}(t)$

From the expression for the time evolution of $|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩$ [16, 21],

(23)$\begin{array}{rl}& \mathrm{i}\hslash {\mathrm{\partial }}_t|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩\\ & =\mathrm{i}\hslash \sum _{k>0}\{{\mathrm{\partial }}_t[c_{\mathit{q},k}^{†}(t)c_{\mathit{q},\overline{k}}^{†}(t)]+{\dot{\mu }}_{\mathit{q},k}(t)+{\dot{\nu }}_{\mathit{q},k}(t)c_{\mathit{q},k}^{†}(t)c_{\mathit{q},\overline{k}}^{†}(t)\}\\ & \cdot \prod _{j\ne k,j>0}[{\mu }_{\mathit{q},j}(t)+{\nu }_{\mathit{q},j}(t)c_{\mathit{q},j}^{†}(t)c_{\mathit{q},\overline{j}}^{†}(t)]|-⟩\\ & =\sum _k^{l_q}[{\hat{h}}^{\mathit{q}}(\mathit{r},t)-{\epsilon }_k^{\mathit{q}}(t)]c_{\mathit{q},k}^{†}(t)c_{\mathit{q},k}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩\\ & +\mathrm{i}\hslash \sum _{k>0}\sqrt{|{\dot{\mu }}_{\mathit{q},k}(t){|}^2+|{\dot{\nu }}_{\mathit{q},k}(t){|}^2}|{\tilde{\mathrm{\Phi }}}_{\mathit{q},k}(t)⟩,\end{array}$

where the Slater determinant $|{\tilde{\mathrm{\Phi }}}_{\mathit{q},k}(t)⟩$ is defined as

(24)$\begin{array}{rl}|{\tilde{\mathrm{\Phi }}}_{\mathit{q},k}(t)⟩=& [\frac{{\dot{\mu }}_{\mathit{q},k}(t)}{\sqrt{|{\dot{\mu }}_{\mathit{q},k}(t){|}^2+|{\dot{\nu }}_{\mathit{q},k}(t){|}^2}}\\ & +\frac{{\dot{\nu }}_{\mathit{q},k}(t)}{\sqrt{|{\dot{\mu }}_{\mathit{q},k}(t){|}^2+|{\dot{\nu }}_{\mathit{q},k}(t){|}^2}}{c}_{\mathit{q},k}^{†}(t)c_{\mathit{q},\overline{k}}^{†}(t)]\\ & \cdot \prod _{j\ne k,j>0}[{\mu }_{\mathit{q},j}(t)+{\nu }_{\mathit{q},j}(t)c_{\mathit{q},j}^{†}(t)c_{\mathit{q},\overline{j}}^{†}(t)]|-⟩.\end{array}$

Eq. (6c) can be written in the form

(25)$\begin{array}{rl}{\mathcal{H}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}^{MF}(t)=& ⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\prime }}(t)|\mathrm{i}\hslash {\mathrm{\partial }}_t|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩\\ =& ⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\prime }}(t)|\sum _k^{l_q}[{\hat{h}}^{\mathit{q}}(\mathit{r},t)-{\epsilon }_k^{\mathit{q}}(t)]c_{\mathit{q},k}^{†}(t)c_{\mathit{q},k}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩\\ & +\mathrm{i}\hslash \sum _{k>0}\sqrt{|{\dot{\mu }}_{\mathit{q},k}(t){|}^2+|{\dot{\nu }}_{\mathit{q},k}(t){|}^2}⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\prime }}(t)|{\tilde{\mathrm{\Phi }}}_{\mathit{q},k}(t)⟩.\end{array}$

By expanding $[{\hat{h}}^{\mathit{q}}(\mathit{r},t)-{\epsilon }_k^{\mathit{q}}(t)]c_{\mathit{q},k}^{†}(t)$ in a complete basis $c_{{\mathit{q}}^{\prime },k^{\prime }}^{†}(t)$,

(26)$\begin{array}{rl}[{\hat{h}}^{\mathit{q}}(\mathit{r},t)-{\epsilon }_k^{\mathit{q}}(t)]c_{\mathit{q},k}^{†}(t)=& \sum _{k^{\prime }}⟨{\varphi }_{k^{\prime }}^{{\mathit{q}}^{\prime }}(\mathit{r},t)|[{\hat{h}}^{\mathit{q}}(\mathit{r},t)\\ & -{\epsilon }_k^{\mathit{q}}(t)]|{\varphi }_k^{\mathit{q}}(\mathit{r},t)⟩c_{{\mathit{q}}^{\prime },k^{\prime }}^{†}(t),\end{array}$

one obtains for ${\mathcal{H}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}^{MF}(t)$ the expression

(27)$\begin{array}{rl}{\mathcal{H}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}^{MF}(t)=& ⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\prime }}(t)|{\mathrm{\Phi }}_{\mathit{q}}(t)⟩\cdot \sum _{k^{\prime }}^{l_{{\mathit{q}}^{\prime }}}\sum _k^{l_{\mathit{q}}}⟨{\varphi }_{k^{\prime }}^{{\mathit{q}}^{\prime }}(\mathit{r},t)|[{\hat{h}}^{\mathit{q}}(\mathit{r},t)\\ & -{\epsilon }_k^{\mathit{q}}(t)]|{\varphi }_k^{\mathit{q}}(\mathit{r},t)⟩{\rho }_{k^{\prime }k}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}(t)\\ & +\mathrm{i}\hslash \sum _{k>0}\sqrt{|{\dot{\mu }}_{\mathit{q},k}(t){|}^2+|{\dot{\nu }}_{\mathit{q},k}(t){|}^2}⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\prime }}(t)|{\tilde{\mathrm{\Phi }}}_{\mathit{q},k}(t)⟩,\end{array}$

where ${\dot{\mu }}_{\mathit{q},k}(t)$ and ${\dot{\nu }}_{\mathit{q},k}(t)$ can be derived from Eq. (12), and $⟨{\mathrm{\Phi }}_{{\mathit{q}}^{\prime }}(t)|{\tilde{\mathrm{\Phi }}}_{\mathit{q},k}(t)⟩$ can be obtained by the Pfaffian algorithms [25, 26].

2.5 Collective wave function $g(t)$

Equation (4) is not a collective Schrödinger equation, and the weight function $f_{\mathit{q}}(t)$ is not a probability amplitude of finding the system at the collective coordinate $\mathit{q}$. The corresponding collective wave function $g_{\mathit{q}}(t)$ is defined by the transformation [28]

(28)$g={\mathcal{N}}^{1/2}f,$

where ${\mathcal{N}}^{1/2}$ is the square root of the overlap kernel matrix. Inserting Eq. (28) into Eq. (4), the time evolution of the collective wave function is governed by the equation [20]

(29)$\mathrm{i}\hslash \dot{g}={\mathcal{N}}^{-1/2}(H-H^{MF}){\mathcal{N}}^{-1/2}g+\mathrm{i}\hslash {\dot{\mathcal{N}}}^{1/2}{\mathcal{N}}^{-1/2}g.$

2.6 Observables $\hat{O}$

The kernel of any observable $\hat{O}$

(30)${\mathcal{O}}_{{\mathit{q}}^{\mathbf{\prime }}\mathit{q}}=⟨{\mathrm{\Phi }}_{q^{\prime }}(t)|\hat{O}|{\mathrm{\Phi }}_q(t)⟩$

can be mapped to the corresponding collective operator ${\mathcal{O}}^c$:

(31)${\mathcal{O}}^c={\mathcal{N}}^{-1/2}\mathcal{O}{\mathcal{N}}^{-1/2}.$

The expectation value of the observable $\hat{O}$ in the correlated GHW state reads

(32)$⟨\mathrm{\Psi }(t)|\hat{O}|\mathrm{\Psi }(t)⟩=f^{†}\mathcal{O}f=g^{†}{\mathcal{O}}^{c}g.$

3 Numerical details

In the present calculation that employs the generalized TD-GCM, the mesh spacing of the lattice is 1.0 fm for all directions, and the box size is $L_x\times L_y\times L_z=20\mathrm{f}\mathrm{m}\times 20\mathrm{f}\mathrm{m}\times 60\mathrm{f}\mathrm{m}$. The time-dependent Dirac equation (8) is solved using the predictor−corrector method, and the time-dependent equations (12) and (29) with the Euler algorithm. The step for the time evolution is $6.67\times {10}^{-4}\mathrm{z}\mathrm{s}$. The point-coupling relativistic energy density functional PC-PK1 [24] is used in the particle−hole channel, together with a monopole pairing interaction. The pairing strength parameters: $-0.135$ MeV for neutrons, and $-0.230$ MeV for protons, are determined by the empirical pairing gaps of ²⁴⁰Pu, using the three-point odd-even mass formula [29]. The initial states for the TD-DFT are obtained by self-consistent deformation-constrained relativistic DFT calculations in a three-dimensional lattice space, using the inverse Hamiltonian and Fourier spectral methods [30–32], with the box size: $L_x\times L_y\times L_z=20\mathrm{f}\mathrm{m}\times 20\mathrm{f}\mathrm{m}\times 50\mathrm{f}\mathrm{m}$.

4 Induced fission of ²⁴⁰Pu

In the analysis of induced fission in Ref. [21], we used a consistent microscopic framework, based on TD-GCM and TD-DFT, to model the entire process of induced fission. The TDGCM was used to evolve adiabatically a set of quadrupole and octupole deformation degrees of freedom of ²⁴⁰Pu from the quasi-stationary initial state to the outer barrier and beyond. Starting from an iso-energy contour beyond the outer barrier, for which the probabilities that the collective wave function reaches these points at any given time were calculated using TDGCM, the TDDFT was subsequently used to model the dissipative fission dynamics in the saddle-to-scission phase. It has been shown, however, that the time-dependent GCM+DFT approach cannot reproduce both the distributions of fission yields and kinetic energy at a quantitative level. The principal reason is that quantum fluctuations, that are essential for a quantitative estimate of fission yields, are included in GCM but not in DFT trajectories that model the saddle-to-scission evolution of the fissioning system.

In the present study, we aim to improve the description of the saddle-to-scission phase of the fission process, by performing generalized time-dependent GCM calculations that use a basis of, generally non-orthogonal and overcomplete, TD-DFT fission trajectories to build the correlated TD-GCM wave function. As shown in the first part of this work [17] and in Ref. [21], we consider the illustrative case of fission of ²⁴⁰Pu. In Fig.1, we plot the two-dimensional microscopic self-consistent mean-field deformation energy surface of ²⁴⁰Pu as a function of the axial quadrupole (${\beta }_{20}$) and octupole (${\beta }_{30}$) deformation parameters. The open dots on the energy surface correspond to the iso-energy contour at 1 MeV below the energy of the equilibrium minimum, located beyond the outer fission barrier. The panel on the right displays the normalized probability distribution that the initial TD-GCM wave packet reaches a particular point after a time of 30 zs, as a function of ${\beta }_{30}$. These points are used as initial locations for the TD-DFT calculations that produce a set of self-consistent fission trajectories to be used as generator states of the generalized TD-GCM. We have plotted 25 trajectories that start from points on the iso-energy contour, reached by the initial collective wave function with non-negligible probability.

Fig.1 Self-consistent deformation energy surface of 240Pu in the plane of quadrupole-octupole axially-symmetric deformation parameters, calculated with the relativistic density functional PC-PK1 and a monopole pairing interaction. Contours join points on the surface with the same energy, and the open dots correspond to points on the iso-energy curve at 1 MeV below the energy of the equilibrium minimum. In the panel on the right, the normalized probability that the initial TD-GCM wave packet reaches a particular point after 30 zs, is plotted as function of the octupole deformation parameter. The curves correspond to self-consistent TD-DFT fission trajectories from the initial points, to be used as generator states of the generalized TD-GCM.

Full size|PPT slide

In principle, one proceeds by building the correlated TD-GCM wave function in the basis of all DFT fission trajectories. At present, however, in most cases this is computationally prohibitive. Therefore, in this illustrative calculation, we have grouped the DFT fission trajectories in sets of 5 neighbouring trajectories, denoted by different colours in Fig.1. Of course, this is not an ideal solution as trajectories from different sets might cross, but will suffice to illustrate the effect of including quantum fluctuations in the description of self-consistent evolution of the fissioning system.

We start by considering the time-evolution of the correlated GCM wave function that starts at the point ${\beta }_{20}=2.91$ and ${\beta }_{30}=2.08$ on the deformation energy surface, and is a superposition of 5 TD-DFT trajectories from region 1 (green trajectories in Fig.1). When these trajectories are used as generator states of the generalized TD-GCM, their initial overlaps are relatively large, as shown by the eigenvalues of the overlap kernel in the top panel of Fig.2. One notices that the eigenvalues approach 1 asymptotically with time, which means that the trajectories become orthogonal. This is because after scission they correspond to distinct pairs of fragments with different particle numbers. The second panel of Fig.2 displays the evolution of the five components of the collective TD-GCM wave function. One notices considerable mixing and, therefore, pronounced quantum fluctuations in the period before scission, while the contributions of the square moduli of components $g$ are constant after scission, reflecting the orthogonality of the TD-DFT basis trajectories. In the lower two panels, we plot the time dependence of the quadrupole and octupole deformations on the way to scission and beyond, in comparison to the single TD-DFT trajectory that also starts at ${\beta }_{20}=2.91$ and ${\beta }_{30}=2.08$. The interesting result here is that, while the single TD-DFT trajectory eventually gets trapped in local minima and does not lead to scission (${\beta }_{20}$ and ${\beta }_{30}$ decrease with time), quantum fluctuations induced by admixing configurations that correspond to the other four TD-DFT trajectories are sufficient to produce a scission event.

Fig.2 Fission of 240Pu with the five TD-DFT trajectories that start in region 1. The initial point for the generalized GCM evolution is at ${\beta }_{20}=2.91$ and ${\beta }_{30}=2.08$ on the deformation energy surface. The top panel displays the eigenvalues of the overlap kernel, and the square moduli of components of the TD-GCM collective wave function are shown in the second panel. The time evolution of the quadrupole and octupole deformations on the way to scission and beyond, is compared to the single TD-DFT trajectory starting from the same point, in the two lower panels.

Full size|PPT slide

In the next example, we examine the superposition of 5 fission trajectories that start in region 3 on Fig.1. This is, of course, not an optimal choice as obviously one expects pronounced admixtures of TD-DFT trajectories from regions 2 and 4 as well, but at present we cannot evolve in time a collective wave function expanded in a much larger TD-DFT basis. The evolution of the correlated GCM wave function starts at the point $({\beta }_{20},{\beta }_{30})=(2.30,1.13)$, and the time dependence of the eigenvalues of the overlap kernel, the square moduli of components of the wave function, and the quadrupole and octupole deformations, are shown in Fig.3. In this particular case, starting from $({\beta }_{20},{\beta }_{30})=(2.30,1.13)$, both the TD-DFT and generalized TD-GCM trajectories lead to scission, though the increase of deformation parameters is obviously steeper for the latter. Scission occurs after approximately 1200 fm/c, after which the TD-DFT basis trajectories are orthogonal and the eigenvalues of the overlap matrix are all equal to 1. At scission, each TD-DFT basis trajectory corresponds to a specific number of protons and neutrons in the fission fragments. This means that, at scission and beyond, the square moduli of components of the TD-GCM collective wave function, shown in the second panel of Fig.3, denote the probabilities to observe the light and heavy fragments of the corresponding TD-DFT trajectory.

Fig.3 Same as in the caption to Fig.2 but for the five TD-DFT trajectories that start in region 3 on Fig.1. The initial point for the generalized GCM evolution is at ${\beta }_{20}=2.30$ and ${\beta }_{30}=1.13$ on the deformation energy surface. The vertical dashed line denotes the instant of scission.

Full size|PPT slide

This is illustrated in Fig.4, where we display the probability distributions of charge yields for the TD-DFT trajectories that start in region 3 on the deformation energy surface in Fig.1. In each panel, from top to bottom, the blue bars denote the charge of the light and heavy fragments calculated for the TD-DFT trajectory that starts from the initial point: $({\beta }_{20},{\beta }_{30})=(2.41,1.30),(2.33,1.20),(2.30,1.13),(2.36,1.00)$, and $(2.50,0.96)$, respectively. Obviously, for each TD-DFT trajectory one obtain just one charge number for the light, and one for the heavy fragment. Since these are average yields, that is, no number projection is performed, they are not necessarily integer numbers. In addition, the values are located in a small interval because all five TD-DFT trajectories start from neighbouring points on the deformation energy surface. The red bars, normalized to 1 for the light and heavy fragments, correspond to the charge yields obtained for the generalized TD-GCM trajectory that starts from the same initial point as the TD-DFT trajectory, but is a superposition of all five basis trajectories. In this case, a distribution of charge yields is obtained, rather than a single light and heavy fragment. This is a direct manifestation of quantum fluctuations that are intrinsic to the generalized TD-GCM.

Fig.4 Probability distributions of charge yields for the TD-DFT trajectories that start in region 3 on the deformation energy surface in Fig.1. In each panel, from top to bottom, the blue bars denote the charge of the light and heavy fragments calculated for the TD-DFT trajectory that starts from the initial points: $({\beta }_{20},{\beta }_{30})=(2.41,1.30),(2.33,1.20),(2.30,1.13),$$(2.36,1.00)$, and $(2.50,0.96)$, respectively. The red bars, normalized to 1 for the light and heavy fragments, correspond to the charge yields obtained for the generalized TD-GCM trajectory that starts from the same initial point as the TD-DFT one, but is a superposition of all five basis trajectories.

Full size|PPT slide

To calculate the total charge yields, we multiply each scission event by the probability that the initial GCM collective wave packet has reached the corresponding starting point of the generalized TD-GCM trajectory after 30 zs (panel on the right of Fig.1). When the calculated fractional charge numbers are assigned to the nearest integer charge, and the yields are normalized to $\sum _{Z}Y(Z)=200$, the resulting total charge yields are plotted in the upper panel of Fig.5, in comparison to the experimental fragment charge distribution [33]. The latter correspond to an average excitation energy of 10.7 MeV. Compared to the data, it is evident that the calculated yields reproduce the position of the peaks but not the tails of the experimental distribution. The reason is the very small probabilities that the initial GCM collective wave packet reaches the points on the iso-energy contour that would lead to the corresponding scission events in the tails of the charge distribution. This could change, however, if the GCM probabilities at the iso-energy curve were calculated at a finite temperature that correspond to the experimental excitation energy of 10.7 MeV. As shown in our studies of Refs. [7, 34], by employing the TD-GCM framework based on a finite-temperature extension of nuclear density functional theory, theoretical yields are predicted in considerably better agreement with available data. In this work, however, we have included pairing correlations but not finite temperature effects.

Fig.5 Charge yields for induced fission of 240Pu. The yields calculated with the generalized TD-GCM (a), and TD-DFT (b) methods, are shown in comparison with the experimental charge distribution. The data are from Ref. [33], and correspond to an average excitation energy of 10.7 MeV.

Full size|PPT slide

Since each correlated TD-GCM trajectory is expanded in a basis of only five neighbouring TD-DFT trajectories, the total charge yields are not very different from those obtained using the standard TD-DFT method [21], shown in the lower panel of Fig.5. Taking into account the results displayed in Fig.4, it is apparent that a much larger basis of TD-DFT trajectories is required for the generalized TD-GCM to produce quantitatively different charge yields.

Similar results are also obtained when the TD-DFT and generalized TD-GCM are compared in a calculation of the total kinetic energy distribution, as shown in Fig.6. For one TD-DFT trajectory, the total kinetic energy (TKE) at a finite distance between the fission fragments (≈ 25 fm, at which shape relaxation brings the fragments to their equilibrium shapes) is calculated using the expression

Fig.6 The calculated total kinetic energies of the nascent fragments for induced fission of 240Pu, as functions of the fragment charge. The generalized TD-GCM and TD-DFT results are shown in comparison to the data [35].

Full size|PPT slide

(33)$E_{\mathrm{T}\mathrm{K}\mathrm{E}}={\displaystyle \frac{1}{2}}m{A}_{\mathrm{H}}^{\mathit{q}}{\mathit{v}}_{\mathrm{H},\mathit{q}}^2+{\displaystyle \frac{1}{2}}m{A}_{\mathrm{L}}^{\mathit{q}}{\mathit{v}}_{\mathrm{L},\mathit{q}}^2+E_{\mathrm{C}\mathrm{o}\mathrm{u}\mathrm{l}}^{\mathit{q}},$

where the velocity of the fragment $f=H,L$ reads

(34)${\mathit{v}}_{f,\mathit{q}}={\displaystyle \frac{1}{m{A}_f^{\mathit{q}}}}{\displaystyle {\int }_{V_f^{\mathit{q}}}\mathrm{d}\mathit{r}{\mathit{j}}_{\mathit{q}}(\mathit{r}),}$

and $\mathit{j}(\mathit{r})$ is the total current density. The integration is over the half-volume corresponding to the fragment $f$, and $E_{\mathrm{C}\mathrm{o}\mathrm{u}\mathrm{l}}$ is the Coulomb energy. For the correlated TD-GCM state, the charge number of the fragments and the total kinetic energy of a trajectory is calculated using Eq. (32). In comparison with the experimental results [35], both methods qualitatively reproduce the experimental TKEs for the fragments close to the peaks of the charge yields distribution, but predict TKEs considerably below the experimental values for the tails of the distribution. As already noted in Ref. [21], this is due to the fact that the calculated TKEs do not include the contribution of pre-scission energy because the initial points for the fission trajectories are on the deformation energy surface, 1 MeV below the energy of the equilibrium minimum and, therefore, more than 10 MeV below the ‘physical’ value. Also in this case, the inclusion of finite temperature effects could improve the agreement with experimental results.

To illustrate the possible effect of including finite temperature in the description of induced fission dynamics (cf. our recent study in Ref. [36]), in Fig.7 we plot the generalized TD-GCM and TD-DFT charge yields, calculated in the same way as those shown in Fig.5, but with equal probabilities for all 25 initial points. This is, of course, exaggerated as even at the temperature that corresponds to the experimental excitation energy, the initial points will not be reached with equal probability by the GCM collective wave packet after $30$ zs. However, we notice a rather pronounced effect of the initial probabilities and a much better agreement with experiment, both for the TD-DFT and generalized TD-GCM.

Fig.7 Same as Fig.5, but with equal initial probabilities for each fission trajectory.

Full size|PPT slide

As already noted in the discussion above, the choice to expand each generalized TD-GCM fission trajectory in a basis of only 5 TD-DFT trajectories is dictated by the available computational capabilities. While we hope that in the future it will become possible to perform more realistic calculations, here the effect is illustrated by merging the initial points in regions 1 and 2 into the basis 1−2 with 10 TD-DFT trajectories, as shown in Fig.8. The trajectories in region 1 produce the fragments in the tail of the charge distribution (charge number $Z=36,37,38$, and $56,57,58$), while the trajectories in region 2 mainly contribute to the peak fragments with charge number $Z=40,41$ and $53,54$.

Fig.8 Same as the left panel of Fig.1, but the initial points in regions 1 and 2 are here merged into a larger basis 1−2, which includes 10 TD-DFT trajectories.

Full size|PPT slide

Starting from one of the ten initial points in the region 1−2, and expanding in the basis of 10 TD-DFT trajectories, the corresponding generalized TD-GCM charge yields are obtained. In each panel of Fig.9, from (a) to (j), the blue bars denote the charge of the light and heavy fragments calculated for the TD-DFT trajectory that starts from an initial point in decreasing order of ${\beta }_{30}$, respectively. In panels (a) and (b) no blue bars are shown, and this means that the corresponding trajectory does not lead to scission and the formation of fission fragments. The red bars, normalized to 1 for the light and heavy fragments, correspond to the charge yields obtained for the generalized TD-GCM trajectory that starts from the same initial point as the TD-DFT trajectory, but is a superposition of all ten basis trajectories. We note that, even in the cases in which the TD-DFT trajectory does not end in scission [panels (a) and (b)], the correlated TD-GCM state starting from the same initial point, but expanded in a basis that includes neighbouring TD-DFT trajectories, produces a broad distribution of fission fragments. When compared with the yields shown in Fig.4, where we have shown results obtained with the 5 TD-DFT basis trajectories in region 3, the advantage of expanding the generalized TD-GCM state in a larger basis is obvious, as it results in a much broader distribution.

Fig.9 Probability distributions of charge yields for the TD-DFT trajectories that start in regions 1−2 on the deformation energy surface in Fig.8. In each panel, from (a) to (j), the blue bars denote the charge of the light and heavy fragments calculated for the TD-DFT trajectory that starts from an initial point in decreasing order of ${\beta }_{30}$, respectively. The red bars, normalized to 1 for the light and heavy fragments, correspond to the charge yields obtained for the generalized TD-GCM trajectory that starts from the same initial point as the TD-DFT trajectroy, but is a superposition of all ten basis trajectories.

Full size|PPT slide

5 Summary

In this work, we have extended the generalized TD-GCM framework, developed in Ref. [17], to include also pairing correlations. The correlated GCM wave function of the nuclear system is determined by time-dependent generator states and weight functions. In the particular implementation used in Ref. [17] and in this work, the initial generator states are calculated using the deformation-constrained self-consistent mean-field method, and are evolved in time by the standard equations of density functional theory (DFT). The DFT generator states form a time-dependent basis in which the GCM wave function is expanded.

The particle−hole channel of the interaction is determined by a Hamiltonian derived from a relativistic (covariant) energy density functional. Here we have employed the point-coupling functional PC-PK1 [24]. Pairing correlations are included in a standard BCS approximation, with time-dependent pairing tensor and occupation probabilities, and a monopole pairing force is used, with strength parameters for protons and neutrons determined by empirical pairing gaps. The inclusion of pairing allows for a more complete modelling of various time-dependent phenomena, but is technically much more complicated and, therefore, computationally demanding.

One of the main motivation for expanding the GCM framework with time-dependent generator states, and also for including pairing correlations, is the possibility for a more realistic description of the final, saddle-to-scission, phase of the process of induced fission. To take into account the dissipation of energy of nuclear collective motion into excitation of intrinsic degrees of freedom, that is, the heating of the system as it evolves toward scission, the dynamics of this phase of fission is usually modelled by methods based on time-dependent density-functional theory. However, since TD-DFT based models only describe a single fission trajectory by propagating nucleons independently from a given initial point on the deformation energy surface, they do not take into account quantum fluctuations that arise from a coherent superposition of different trajectories. In the generalized TD-GCM, with TD-DFT trajectories that start from different points on the deformation energy surface as time-dependent basis states in which the correlated nuclear function is expanded, the description of fission dynamics goes beyond the adiabatic approximation of the standard GCM, and also includes quantum fluctuations.

In the present illustrative study, we have applied the new model to the description of saddle-to-scission dynamics of induced fission of ²⁴⁰Pu. The charge yields and total kinetic energy distribution, obtained by the generalized TD-GCM, have been compared to those calculated using the standard TD-DFT approach, and with available data. Although very encouraging results have been obtained for the principal charge yields and the corresponding kinetic energies, a detailed comparison with data indicates the necessity to expand the GCM correlated function in a relatively large basis of TD-DFT fission trajectories. This is, at present, still computationally prohibitive and will require the utilisation of large computing facilities and/or more efficient algorithms. In addition, to take into account the initial excitation energy of the fissioning system, the theoretical formalism will have to be extended to finite temperatures. Work along these lines is in progress and will be reported in a forthcoming publication.

References

[1] H.J. KrappeK. Pomorski, Theory of Nuclear Fission, Berlin, Heidelberg: Springer, 2012

[2] N. Schunck, L. M. Robledo. Microscopic theory of nuclear fission: A review. Rep. Prog. Phys., 2016, 79(11): 116301

[3] W.YounesD. M. GognyJ.F. Berger, A Microscopic Theory of Fission Dynamics Based on the Generator Coordinate Method, Springer Cham, 2019

[4] D. Regnier, N. Dubray, N. Schunck, M. Verrière. Fission fragment charge and mass distributions in ²³⁹Pu(n, f) in the adiabatic nuclear energy density functional theory. Phys. Rev. C, 2016, 93(5): 054611

[5] M. Verriere, D. Regnier. The time-dependent generator coordinate method in nuclear physics. Front. Phys. (Lausanne), 2020, 8: 233

[6] H. Tao, J. Zhao, Z. P. Li, T. Nikšić, D. Vretenar. Microscopic study of induced fission dynamics of ²²⁶Th with covariant energy density functionals. Phys. Rev. C, 2017, 96(2): 024319

[7] J. Zhao, J. Xiang, Z. P. Li, T. Nikšić, D. Vretenar, S. G. Zhou. Time-dependent generator-coordinate-method study of mass-asymmetric fission of actinides. Phys. Rev. C, 2019, 99(5): 054613

[8] C. Simenel, A. Umar. Heavy-ion collisions and fission dynamics with the time-dependent Hartree−Fock theory and its extensions. Prog. Part. Nucl. Phys., 2018, 103: 19

[9] T. Nakatsukasa, K. Matsuyanagi, M. Matsuo, K. Yabana. Time-dependent density-functional description of nuclear dynamics. Rev. Mod. Phys., 2016, 88(4): 045004

[10] P. Stevenson, M. Barton. Low-energy heavy-ion reactions and the Skyrme effective interaction. Prog. Part. Nucl. Phys., 2019, 104: 142

[11] A. Bulgac, P. Magierski, K. J. Roche, I. Stetcu. Induced fission of ²⁴⁰Pu within a real-time microscopic framework. Phys. Rev. Lett., 2016, 116(12): 122504

[12] P. Magierski, K. Sekizawa, G. Wlazlowski. Novel role of superfluidity in low-energy nuclear reactions. Phys. Rev. Lett., 2017, 119(4): 042501

[13] G. Scamps, C. Simenel. Impact of pear-shaped fission fragments on mass-asymmetric fission in actinides. Nature, 2018, 564(7736): 382

[14] A. Bulgac, S. Jin, K. J. Roche, N. Schunck, I. Stetcu. Fission dynamics of ²⁴⁰Pu from saddle to scission and beyond. Phys. Rev. C, 2019, 100(3): 034615

[15] A. Bulgac, S. Jin, I. Stetcu. Nuclear fission dynamics: Past, present, needs, and future. Front. Phys. (Lausanne), 2020, 8: 63

[16] Z. X. Ren, D. Vretenar, T. Nikšić, P. W. Zhao, J. Zhao, J. Meng. Dynamical synthesis of ⁴He in the scission phase of nuclear fission. Phys. Rev. Lett., 2022, 128(17): 172501

[17] B. Li, D. Vretenar, T. Nikšić, P. W. Zhao, J. Meng. Generalized time-dependent generator coordinate method for small- and large-amplitude collective motion. Phys. Rev. C, 2023, 108(1): 014321

[18] P. Marević, D. Regnier, D. Lacroix. Quantum fluctuations induce collective multiphonons in finite Fermi liquids. Phys. Rev. C, 2023, 108(1): 014620

[19] P. G. Reinhard, R. Cusson, K. Goeke. Time evolution of coherent ground-state correlations and the TDHF approach. Nucl. Phys. A, 1983, 398(1): 141

[20] D. Regnier, D. Lacroix. Microscopic description of pair transfer between two superfluid Fermi systems (ii): quantum mixing of time-dependent Hartree−Fock−Bogolyubov trajectories. Phys. Rev. C, 2019, 99(6): 064615

[21] Z. X. Ren, J. Zhao, D. Vretenar, T. Nikšić, P. W. Zhao, J. Meng. Microscopic analysis of induced nuclear fission dynamics. Phys. Rev. C, 2022, 105(4): 044313

[22] S. Ebata, T. Nakatsukasa, T. Inakura, K. Yoshida, Y. Hashimoto, K. Yabana. Canonical-basis time-dependent Hartree−Fock−Bogoliubov theory and linear-response calculations. Phys. Rev. C Nucl. Phys., 2010, 82(3): 034306

[23] G. Scamps, D. Lacroix. Effect of pairing on one- and two-nucleon transfer below the coulomb barrier: A time-dependent microscopic description. Phys. Rev. C, 2013, 87(1): 014605

[24] P. W. Zhao, Z. P. Li, J. M. Yao, J. Meng. New parametrization for the nuclear covariant energy density functional with a point-coupling interaction. Phys. Rev. C, 2010, 82(5): 054319

[25] L. M. Robledo. Sign of the overlap of Hartree−Fock−Bogoliubov wave functions. Phys. Rev. C, 2009, 79(2): 021302

[26] Q. L. Hu, Z. C. Gao, Y. Chen. Matrix elements of one-body and two-body operators between arbitrary HFB multi-quasiparticle states. Phys. Lett. B, 2014, 734: 162

[27] P. Bonche, J. Dobaczewski, H. Flocard, P. H. Heenen, J. Meyer. Analysis of the generator coordinate method in a study of shape isomerism in ¹⁹⁴Hg. Nucl. Phys. A, 1990, 510(3): 466

[28] P. G. Reinhard, K. Goeke. The generator coordinate method and quantised collective motion in nuclear systems. Rep. Prog. Phys., 1987, 50(1): 1

[29] M. Bender, K. Rutz, P. G. Reinhard, J. A. Maruhn. Pairing gaps from nuclear mean field models. Eur. Phys. J. A, 2000, 8(1): 59

[30] Z. X. Ren, S. Q. Zhang, J. Meng. Solving Dirac equations on a 3D lattice with inverse hamiltonian and spectral methods. Phys. Rev. C, 2017, 95(2): 024313

[31] Z. X. Ren, S. Q. Zhang, P. W. Zhao, N. Itagaki, J. A. Maruhn, J. Meng. Stability of the linear chain structure for ¹²C in covariant density functional theory on a 3D lattice. Sci. China Phys. Mech. Astron., 2019, 62(11): 112062

[32] Z. X. Ren, P. W. Zhao, S. Q. Zhang, J. Meng. Toroidal states in ²⁸Si with covariant density functional theory in 3D lattice space. Nucl. Phys. A, 2020, 996: 121696

[33] D. Ramos, M. Caamaño, F. Farget, C. Rodríguez-Tajes, L. Audouin, J. Benlliure, E. Casarejos, E. Clement, D. Cortina, O. Delaune, X. Derkx, A. Dijon, D. Doré, B. Fernańdez-Domínguez, G. de France, A. Heinz, B. Jacquot, A. Navin, C. Paradela, M. Rejmund, T. Roger, M. D. Salsac, C. Schmitt. Isotopic fission-fragment distributions of ²³⁸U, ²³⁹Np, ²⁴⁰Pu, ²⁴⁴Cm, and ²⁵⁰Cf produced through inelastic scattering, transfer, and fusion reactions in inverse kinematics. Phys. Rev. C, 2018, 97(5): 054612

[34] J. Zhao, T. Nikšić, D. Vretenar, S. G. Zhou. Microscopic self-consistent description of induced fission dynamics: Finite-temperature effects. Phys. Rev. C, 2019, 99(1): 014618

[35] M. Caamaño, F. Farget, O. Delaune, K. H. Schmidt, C. Schmitt, L. Audouin, C. O. Bacri, J. Benlliure, E. Casarejos, X. Derkx, B. Fernańdez-Domínguez, L. Gaudefroy, C. Golabek, B. Jurado, A. Lemasson, D. Ramos, C. Rodríguez-Tajes, T. Roger, A. Shrivastava. Characterization of the scission point from fission-fragment velocities. Phys. Rev. C, 2015, 92(3): 034606

[36] B. Li, D. Vretenar, Z. X. Ren, T. Nikšić, J. Zhao, P. W. Zhao, J. Meng. Fission dynamics, dissipation, and clustering at finite temperature. Phys. Rev. C, 2023, 107(1): 014303

Declarations

The authors declare that they have no competing interests and there are no conflicts.

Acknowledgements

This work was supported in part by the High-end Foreign Experts Plan of China, the National Key R&D Program of China (Contract No. 2018YFA0404400), the National Natural Science Foundation of China (Grant Nos. 12070131001, 11875075, 11935003, 11975031, and 12141501), the High-performance Computing Platform of Peking University, the QuantiXLie Centre of Excellence, a project co-financed by the Croatian Government and European Union through the European Regional Development Fund − the Competitiveness and Cohesion Operational Programme (KK.01.1.1.01.0004), and the Croatian Science Foundation under the project Uncertainty quantification within the nuclear energy density framework (IP-2018-01-5987).

RIGHTS & PERMISSIONS

2024 Higher Education Press