gbsn

# Microscopic self-consistent description of induced fission: dynamical
pairing degree of freedom

Jie Zhao (赵杰)  Center for Circuits and Systems, Peng Cheng Laboratory,
Shenzhen 518055, China  Tamara Nikšić  Physics Department, Faculty of Science,
University of Zagreb, Bijenička Cesta 32, Zagreb 10000, Croatia  Dario
Vretenar  Physics Department, Faculty of Science, University of Zagreb,
Bijenička Cesta 32, Zagreb 10000, Croatia  State Key Laboratory of Nuclear
Physics and Technology, School of Physics, Peking University, Beijing 100871,
China

###### Abstract

The role of dynamical pairing in induced fission dynamics is investigated
using the time-dependent generator coordinate method in the Gaussian overlap
approximation, based on the microscopic framework of nuclear energy density
functionals. A calculation of fragment charge yields for induced fission of
228Th is performed in a three-dimensional space of collective coordinates
that, in addition to the axial quadrupole and octupole intrinsic deformations
of the nuclear density, also includes an isoscalar pairing degree of freedom.
It is shown that the inclusion of dynamical pairing has a pronounced effect on
the collective inertia, the collective flux through the scission hyper-
surface, and the resulting fission yields, reducing the asymmetric peaks and
enhancing the contribution of symmetric fission, in better agreement with the
empirical trend.

##  I Introduction

Nuclear density functional theory (DFT) is the only microscopic framework that
can be used over the entire table of nuclides in a self-consistent description
of phenomena ranging from ground-state properties and collective excitations,
to large-amplitude nucleonic motion, fission and low-energy collisions. In the
case of nuclear fission Bender  _et al._ (2020), in particular, a fully
quantum mechanical many-body model can be constructed starting from the time-
dependent generator coordinate method (TDGCM) Berger  _et al._ (1991). In this
approach the nuclear wave function is represented by a linear superposition of
many-body generator states that are functions of collective coordinates. In
most cases these coordinates parameterize the shape of the nuclear density.
The Hill-Wheeler equation of motion determines the time evolution of the wave
function in the restricted space of generator states Verriere and Regnier
(2020). By employing the Gaussian overlap approximation (GOA), the GCM Hill-
Wheeler equation reduces to a local, time-dependent Schrödinger equation in
the space of collective coordinates. The microscopic input for the collective
Schrödinger equation, that is, the nuclear potential and collective inertia,
are determined by self-consistent mean-field calculations for a choice of the
energy density functional (EDF) or effective interaction. The TDGCM+GOA method
can be applied to the dynamics of induced fission, starting from the ground
state and following the time evolution of collective degrees of freedom all
the way to scission and the emergence of fission fragments. This framework has
been very successfully implemented in a number of fission studies based on
nonrelativistic Skyrme and Gogny functionals Berger  _et al._ (1991); Verriere
and Regnier (2020); Goutte  _et al._ (2005); Regnier  _et al._ (2016a, 2018,
b, 2017); Verriere  _et al._ (2017); Regnier  _et al._ (2019); Zdeb  _et al._
(2017). These studies have investigated the dependence of the predicted
fission dynamics on the choice of the EDF, initial conditions, form of the
collective inertia, and the definition of scission configurations.

Relativistic energy density functionals Vretenar  _et al._ (2005); Meng  _et
al._ (2006); Meng (2016) have also been employed in the description of both
spontaneous Zhao  _et al._ (2015, 2016) and induced nuclear fission Tao  _et
al._ (2017); Zhao  _et al._ (2019a, b, 2020). The microscopic input for these
studies is generated using either the multidimensionally constrained
relativistic mean-field (MDC-RMF) Lu  _et al._ (2014) or the relativistic
Hartree-Bogoliubov model (MDC-RHB) Zhao  _et al._ (2017). By employing the
TDGCM+GOA collective model, several interesting topics have been explored in
this framework, such as the influence of static pairing correlations on
fission yields, different approximations for the collective inertia tensor,
and finite temperature effects.

Most applications of the TDGCM to fission dynamics have considered a two-
dimensional space of collective coordinates such as, for instance, quadrupole
and octupole shape degrees of freedom. The recently developed computer code
FELIX Regnier  _et al._ (2016a, 2018) offers the possibility of solving the
TDGCM+GOA equation for an arbitrary number of collective variables. In Ref.
Regnier  _et al._ (2017) a preliminary calculation of induced fission dynamics
of 240Pu isotope in the three-dimensional space of shape variables
(quadrupole, octupole and hexadecupole intrinsic deformations) has been
reported. In particular, this model can also be used for a quantitative
analysis of the critical role of dynamical pairing correlations in the process
of induced fission.

The importance of pairing correlations for both spontaneous and induced
nuclear fission has been emphasized in a number of studies Negele  _et al._
(1978); Łojewski and Staszczak (1999); Sadhukhan  _et al._ (2014); Staszczak
_et al._ (1985); Zhao  _et al._ (2016); Bernard  _et al._ (2019); Lazarev
(1987); Sadhukhan  _et al._ (2016); Bulgac  _et al._ (2016); Qiang  _et al._
(2021). For spontaneous fission it has been shown that the coupling between
shape and pairing degrees of freedom has a pronounced effect on the calculated
fission lifetimes Sadhukhan  _et al._ (2014); Zhao  _et al._ (2016). In
particular, when the gap parameter is considered as a dynamical variable,
pairing correlations are generally enhanced thus reducing the effective
inertia and the action integral along the fission path. This effect can
significantly reduce the estimated spontaneous fission lifetimes, and it has
also been noted that pairing fluctuations can restore axial symmetry in the
fissioning system. A study of induced fission of 240Pu, using the microscopic
time-dependent superfluid local density approximation (TD SLDA) Bulgac  _et
al._ (2016), has shown that both shape and pairing modes determine the
dynamics of the final stage of the fission process, from configurations close
to the outer fission barrier to full scission.

The influence of ground-state (static) pairing correlations on charge yields
and total kinetic energy of fission fragments for the case of induced fission
of 226Th isotope was analyzed in Ref. Tao  _et al._ (2017) using the TDGCM+GOA
framework. It has been shown that an increase of the strength of the pairing
interaction, beyond the range determined by empirical pairing gaps obtained
from the experimental masses of neighboring nuclei, reduces the asymmetric
peaks and enhances the symmetric peak in charge yields distribution. This is a
very interesting result, and thus it is important to explore dynamical pairing
correlations in induced fission. In this work we explicitly include the
isoscalar pairing degree of freedom in the space of TDGCM+GOA collective
coordinates, and perform the first realistic three-dimensional calculation of
induced fission of 228Th. The theoretical framework and methods are reviewed
in Sec. II. The details of the calculation and principal results are discussed
in Sec. III. Section IV contains a short summary and outlook for future
studies.

##  II The TDGCM+GOA method

In the TDGCM+GOA framework induced fission is described as a slow adiabatic
process determined by a small number of collective degrees of freedom. The
initial step in modeling the fission of a heavy nucleus is a self-consistent
mean-field (SCMF) calculation of the corresponding deformation energy surface
as a function of few selected collective coordinates. Such a calculation
provides the microscopic input, that is, the single-quasiparticle states,
energies, and occupation factors, that determine the parameters of a local
equation of motion for the collective wave function.

As in our previous studies, here we use the point-coupling relativistic energy
density functional DD-PC1 Nikšić  _et al._ (2008) in the particle-hole
channel, while pairing correlations are taken into account in the Bardeen-
Cooper-Schrieffer (BCS) approximation by a separable pairing force of finite
range Tian  _et al._ (2009):

$V(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime})=G_{0}~{}\delta(\mathbf{R}-\mathbf{R}^{\prime})P(\mathbf{r})P(\mathbf{r}^{\prime})\frac{1}{2}\left(1-P^{\sigma}\right),$
(1)

where $\mathbf{R}=(\mathbf{r}_{1}+\mathbf{r}_{2})/2$ and
$\mathbf{r}=\mathbf{r}_{1}-\mathbf{r}_{2}$ denote the center-of-mass and the
relative coordinates, respectively. $P(\mathbf{r})$ reads

$P(\mathbf{r})=\frac{1}{\left(4\pi
a^{2}\right)^{3/2}}e^{-\mathbf{r}^{2}/4a^{2}}.$ (2)

The parameters of the interaction were originally adjusted to reproduce the
density dependence of the pairing gap in nuclear matter at the Fermi surface
computed with the D1S parameterization of the Gogny force Berger  _et al._
(1991). To reproduce the empirical pairing gaps in the mass region considered
in the present study, the strength parameters of the pairing force have been
increased with respect to the original values by the following factors:
$G_{n}/G_{0}=1.12$ and $G_{p}/G_{0}=1.08$ for neutrons and protons,
respectively.

In the present analysis the self-consistent deformation energy surfaces are
calculated using the multidimensionally constrained relativistic mean-field
(MDC-RMF) model Lu  _et al._ (2012, 2014); Zhou (2016); Zhao  _et al._ (2016)
with constraints on mass multipole moments
$Q_{\lambda\mu}=r^{\lambda}Y_{\lambda\mu}$, and the particle-number dispersion
operator $\Delta\hat{N}^{2}=\hat{N}^{2}-\langle\hat{N}\rangle^{2}$. The
Routhian is therefore defined as

$E^{\prime}=E_{\rm{RMF}}+\sum_{\lambda\mu}\frac{1}{2}C_{\lambda\mu}Q_{\lambda\mu}+\lambda_{2}\Delta\hat{N}^{2},$
(3)

where $E_{\rm RMF}$ denotes the total RMF energy that includes static BCS
pairing correlations. The amount of dynamic pairing correlations can be
controlled by the Lagrange multipliers $\lambda_{2\tau}$ ($\tau=n,p$), Vaquero
_et al._ (2011, 2013); Sadhukhan  _et al._ (2014); Zhao  _et al._ (2016). To
reduce the number of collective degrees of freedom and, therefore, the
considerable computational task, here we only consider isoscalar dynamical
pairing; $\lambda_{2n}=\lambda_{2p}\equiv\lambda_{2}$ is employed as the
collective coordinate. The nuclear shape is parameterized by the deformation
parameters

$\beta_{\lambda\mu}={4\pi\over 3AR^{\lambda}}\langle Q_{\lambda\mu}\rangle,$
(4)

with $R=1.2~{}A^{1/3}$ fm. The shape is assumed to be invariant under the
exchange of the $x$ and $y$ axes, and all deformations $\beta_{\lambda\mu}$
with even $\mu$ can be included simultaneously. The constrained RMF+BCS
equations are solved by an expansion in the axially deformed harmonic
oscillator (ADHO) basis Gambhir  _et al._ (1990). In the present study
calculations have been performed in an ADHO basis truncated to $N_{f}=20$
oscillator shells.

The dynamics of the fission process is thus governed by a local, time-
dependent Schrödinger-like equation in the space of collective coordinates
$\bm{q}$:

$i\hbar\frac{\partial}{\partial t}g(\bm{q},t)=\hat{H}_{\rm
coll}(\bm{q})g(\bm{q},t),$ (5)

where $g(\bm{q},t)$ is the complex wave function of the collective variables
$\bm{q}$ and time $t$. The collective Hamiltonian $\hat{H}_{\rm coll}(\bm{q})$

$\displaystyle\hat{H}_{\rm coll}(\bm{q})=-{\hbar^{2}\over
2}\sum_{ij}{\partial\over\partial q_{i}}B_{ij}(\bm{q}){\partial\over\partial
q_{j}}+V(\bm{q}),$ (6)

governs the time-evolution of the nuclear wave function from an initial state
at equilibrium deformation, up to scission and the formation of fission
fragments. $B_{ij}(\bm{q})$ and $V(\bm{q})$ denote the inertia tensor and
collective potential, respectively, that are computed using the self-
consistent solutions (total energy, single-quasiparticle states and occupation
factors) for the RMF+BCS deformation energy surface. Here we assume axial
symmetry with respect to the axis along which the two fragments eventually
separate, and consider the three-dimensional (3D) collective space of
quadrupole $\beta_{2}$ and octupole $\beta_{3}$ deformation parameters, and
the dynamical pairing coordinate $\lambda_{2}$. The inertia tensor is the
inverse of the mass tensor, that is,
$B_{ij}(\bm{q})=(\mathcal{M}^{-1})_{ij}(\bm{q})$. The mass tensor is
calculated using the adiabatic time-dependent Hartree-Fock-Bogoliubov (ATDHFB)
method in the cranking approximation Baran  _et al._ (2011):

$\mathcal{M}_{ij}^{C}={\hbar^{2}\over
2\dot{q}_{i}\dot{q}_{j}}\sum_{\mu\nu}{F^{i*}_{\mu\nu}F^{j}_{\mu\nu}+F^{i}_{\mu\nu}F^{j*}_{\mu\nu}\over
E_{\mu}+E_{\nu}},$ (7)

where

${F^{i}\over\dot{q}_{i}}=U^{\dagger}{\partial\rho\over\partial
q_{i}}V^{*}+U^{\dagger}{\partial\kappa\over\partial
q_{i}}U^{*}-V^{\dagger}{\partial\rho^{*}\over\partial
q_{i}}U^{*}-V^{\dagger}{\partial\kappa^{*}\over\partial q_{i}}V^{*}\;.$ (8)

$U$ and $V$ are the self-consistent Bogoliubov matrices, and $\rho$ and
$\kappa$ are the corresponding particle and pairing density matrices,
respectively. The derivatives of the densities are calculated using the
Lagrange three-point formula for unequally spaced points Yuldashbaeva  _et
al._ (1999); Baran  _et al._ (2011). The cranking expression Eq. (7) can be
further simplified in the perturbative approach Brack  _et al._ (1972);
Nilsson  _et al._ (1969); Girod and Grammaticos (1979); Bés and Szymański
(1961); Sobiczewski  _et al._ (1969), and this leads to the perturbative
cranking mass tensor:

$\mathcal{M}^{Cp}=\hbar^{2}{\it M}_{(1)}^{-1}{\it M}_{(3)}{\it M}_{(1)}^{-1},$
(9)

where

$\left[{\it M}_{(k)}\right]_{ij}=\sum_{\mu\nu}{\left\langle
0\left|\hat{Q}_{i}\right|\mu\nu\right\rangle\left\langle\mu\nu\left|\hat{Q}_{j}\right|0\right\rangle\over(E_{\mu}+E_{\nu})^{k}}.$
(10)

$|\mu\nu\rangle$ are the two-quasiparticle states with the corresponding
quasiparticle energies $E_{\mu}$ and $E_{\nu}$. The details of the derivation
of cranking formulas for the mass tensor can be found in Ref. Baran  _et al._
(2011).

At this point we introduce an approximation that is not entirely consistent,
but is necessary to reduce the computational task and stabilize the time-
evolution of the collective state. Because we consider the particle-number
dispersion operator $\Delta\hat{N}^{2}=\hat{N}^{2}-\langle\hat{N}\rangle^{2}$,
the pairing part of the mass tensor must be calculated using the non-
perturbative cranking expression Eq. (7). In the recent study of the
differences between the perturbative and non-perturbative ATDHFB collective
masses in the TDGCM+GOA description of induced fission Zhao  _et al._ (2020),
using the axial quadrupole and octupole intrinsic deformation as dynamical
variables, we have shown that the structure of non-perturbative collective
masses is much more complex due to changes in the intrinsic shell structure
across the deformation energy surface, and is characterized by pronounced
isolated peaks located at single-particle level crossings near the Fermi
surface. It has been shown that the choice of non-perturbative cranking
collective mass leads to a reduction of symmetric charge yields and,
generally, to a better agreement with data. Even though in the analysis of
Ref. Zhao  _et al._ (2020) both non-perturbative and perturbative mass tensors
were used in modeling induced fission dynamics, in the present 3D study the
number of mesh points required to accurately calculate all the non-
perturbative collective masses becomes prohibitively large for the available
computational resources. The reason is the occurrence of pronounced peaks in
the collective masses related to single-particle level crossings near the
Fermi surface, and the corresponding abrupt changes of the occupation factors
of single-particle configurations. This would require a major refinement of
the grid, leading to very large number of mesh points and possible
instabilities. Therefore, to be able to quantitatively analyze the effect of
dynamical pairing on induced fission, here we use the perturbative cranking
expression Eq. (9) to calculate the elements of the mass tensor that
correspond to the quadrupole and octupole deformations while, as noted above,
the pairing element of the mass tensor must be computed using the non-
perturbative cranking formula. With the indices $1$, $2$ and $3$ corresponding
to the $\beta_{20}$, $\beta_{30}$, and $\lambda_{2}$ collective coordinates,
the following elements of the mass tensor are used: $\mathcal{M}^{Cp}_{11}$,
$\mathcal{M}^{Cp}_{12}$, $\mathcal{M}^{Cp}_{22}$, $\mathcal{M}^{C}_{33}$, and
we neglect the coupling terms $\mathcal{M}_{13}$ and $\mathcal{M}_{23}$.

To model the fission dynamics we follow the time-evolution of an initial wave
packet $g(\bm{q},t=0)$ ($\bm{q}\equiv\\{\beta_{2},\beta_{3},\lambda_{2}\\}$),
built as a Gaussian superposition of the quasi-bound states $g_{k}$,

$g(\bm{q},t=0)=\sum_{k}\exp\left({(E_{k}-\bar{E})^{2}\over
2\sigma^{2}}\right)g_{k}(\bm{q}),$ (11)

where the value of the parameter $\sigma$ is set to 0.5 MeV. The collective
states $\\{g_{k}(\bm{q})\\}$ are solutions of the stationary eigenvalue
equation in which the original collective potential $V(\bm{q})$ is replaced by
a new potential $V^{\prime}(\bm{q})$ that is obtained by extrapolating the
inner potential barrier with a quadratic form. The mean energy $\bar{E}$ in
Eq. (11) is then adjusted iteratively in such a way that $\langle
g(t=0)|\hat{H}_{\rm coll}|g(t=0)\rangle=E_{\rm coll}^{*}$, and this average
energy $E_{\rm coll}^{*}$ is set to be 1 MeV above the fission barrier. The
TDGCM+GOA Hamiltonian of Eq. (6), with the original collective potential
$V(\bm{q})$, propagates the initial wave packet in time.

The time propagation is modeled using the TDGCM+GOA computer code FELIX
(version 2.0) Regnier  _et al._ (2018). The time step is $\delta t=5\times
10^{-4}$ zs (1 zs $=10^{-21}$ s), and the charge and mass distributions are
calculated after $4\times 10^{4}$ time steps, which correspond to 20 zs. As in
our recent calculations of Refs. Tao  _et al._ (2017); Zhao  _et al._ (2019a,
b, 2020), the parameters of the additional imaginary absorption potential that
takes into account the escape of the collective wave packet in the domain
outside the region of calculation Regnier  _et al._ (2018) are: the absorption
rate $r=20\times 10^{22}$ s-1 and the width of the absorption band $w=1.0$.

The collective space is divided into an inner region with a single nuclear
density distribution, and an external region that contains two separated
fission fragments. The scission hyper-surface that divides the inner and
external regions is determined by the expectation values of the Gaussian neck
operator $\displaystyle\hat{Q}_{N}=\exp[-(z-z_{N})^{2}/a_{N}^{2}]$, where
$a_{N}=1$ fm and $z_{N}$ is the position of the neck Younes and Gogny (2009).
We define the pre-scission domain by $\langle\hat{Q}_{N}\rangle>3$, and
consider the frontier of this domain as the scission surface. The flux of the
probability current through this hyper-surface provides a measure of the
probability of observing a given pair of fragments at time $t$. Each
infinitesimal surface element is associated with a given pair of fragments
$(A_{L},A_{H})$, where $A_{L}$ and $A_{H}$ denote the lighter and heavier
fragments, respectively. The integrated flux $F(\xi,t)$ for a given surface
element $\xi$ is defined as Regnier  _et al._ (2018)

$F(\xi,t)=\int_{t_{0}}^{t}dt^{\prime}\int_{\bm{q}\in\xi}\bm{J}(\bm{q},t^{\prime})\cdot
d\bm{S},$ (12)

where $\bm{J}(\bm{q},t)$ is the current

$\displaystyle
J_{k}(\bm{q},t)=\hbar\sum_{l}B_{kl}(\bm{q}){\mathrm{Im}}\left(g^{*}\frac{\partial
g}{\partial q_{l}}\right).$ (13)

The yield for the fission fragment with mass $A$ is defined by

$Y(A)\propto\sum_{\xi\in\mathcal{A}}\lim_{t\rightarrow\infty}F(\xi,t).$ (14)

The set $\mathcal{A}(\xi)$ contains all elements belonging to the scission
hyper-surface such that one of the fragments has mass number $A$. The charge
yields are obtained by convoluting the raw flux with a Gaussian function of
the number of particles Regnier  _et al._ (2016b); Zhao  _et al._ (2019a),
with a width of 1.6 units.

##  III Results and discussion

As an illustrative example, the fission of 228Th is considered. For this
nucleus the charge distribution of fission fragments exhibits a coexistence of
symmetric and asymmetric peaks Schmidt  _et al._ (2000). In the first step a
large scale MDC-RMF calculation is performed to generate the potential energy
surface, single-nucleon wave functions and occupation factors in the
$(\beta_{2},\beta_{3},\lambda_{2})$ collective space. The intervals for the
values of the collective variables are: $-1\leq\beta_{2}\leq 7$ with a step
$\Delta\beta_{2}=0.04$; $0\leq\beta_{3}\leq 3.5$ with a step
$\Delta\beta_{3}=0.05$; and $-0.1\leq\lambda_{2}\leq 2.0$ with a step
$\Delta\lambda_{2}=0.1$. The relativistic energy density functional DD-PC1
Nikšić  _et al._ (2008) is used in the particle-hole channel, and particle-
particle correlations are described by the separable pairing force (1) in the
BCS approximation.

Figure 1 displays the 2D projections of the 3D deformation energy surface of
228Th on the collective plane $(\beta_{2},\beta_{3})$, for several values of
the collective coordinate $\lambda_{2}$. Note that the value $\lambda_{2}=0$
corresponds to static BCS pairing, while positive $\lambda_{2}$ equates to
enhanced pairing correlations. Only configurations with $Q_{N}\geq 3$ are
plotted, and the frontier of this domain determines the scission contour. The
red curves correspond to static fission paths of minimum total energy. The
deformation surfaces for $\lambda_{2}=0$ and $\lambda_{2}=0.3$ are very
similar, with a pronounced ridge separating the asymmetric and symmetric
fission valleys. This ridge decreases with increasing values of the pairing
coordinate $\lambda_{2}$. The scission contour for $\lambda_{2}=0$ starts from
an elongated symmetric point at $\beta_{2}\approx 6$, and evolves to a minimal
elongation $\beta_{2}\approx 3$ as reflection asymmetry increases. For larger
values of $\lambda_{2}$ the scission contour is not modified significantly,
although the starting point of the scission contour on the quadrupole axis
shifts to smaller elongations at $\beta_{2}\approx 5$.

Figure 1: (Color online) 2D projections of the deformation-energy manifold of
228Th on the quadrupole-octupole axially symmetric plane, calculated with the
RMF+BCS model based on the functional DD-PC1, for selected values of the
pairing coordinate $\lambda_{2}$. Contours join points on the surface with the
same energy, and the separation between neighboring contours is 1 MeV. The red
curves denote static fission paths of minimum energy for each value of
$\lambda_{2}$.

To illustrate the dynamical effect on the pairing correlations, in Fig. 2 we
display the values of the neutron (upper panel) and proton (lower panel)
pairing gap along the static fission paths, as functions of the quadrupole
coordinate $\beta_{2}$, for several values of the isoscalar pairing collective
coordinate $\lambda_{2}$ (cf. Fig. 1). In the interval of values of
$\lambda_{2}$ considered here, the values of the pairing gaps increase by a
factor $\approx 2-3$ and, characteristically, for $\lambda_{2}\geq 0.6$ all
traces of the underlying shell effects along the static fission path vanish.

Figure 2: (Color online) Pairing gaps for neutrons $\Delta_{n}$ (upper panel)
and protons $\Delta_{p}$ (lower panel) along the static fission path, as
functions of the axial quadrupole deformation, for selected values of the
isoscalar pairing collective coordinate $\lambda_{2}$.

In Fig. 3 we plot the values of the perturbative cranking collective masses
$M^{Cp}_{11}$, $M^{Cp}_{22}$, and the non-perturbative cranking mass
$M^{C}_{33}$, along the static fission path as functions of the quadrupole
deformation, for different pairing collective coordinates $\lambda_{2}$. As
noted above, the indices $1$, $2$ and $3$ refer to the $\beta_{2}$,
$\beta_{3}$ and $\lambda_{2}$ coordinates, respectively. The collective masses
$M^{Cp}_{11}$ and $M^{Cp}_{22}$ exhibit a more complex structure for smaller
values of deformations $\beta_{2}<1$ and decrease for larger deformations. One
notices that $M^{Cp}_{22}$ exhibits a sharper decrease for $\beta_{2}>1$,
compared to $M^{Cp}_{11}$. The behaviour of the $M^{C}_{33}$ collective mass
is much more complex, due to the fact that it has been calculated using the
non-perturbative cranking formula Eq. (7). This leads to the appearance of
prominent peaks at the locations of single-particle level crossing near the
Fermi surface. However, in general also $M^{C}_{33}$ decreases for large
deformations and the peaks are less pronounced. As the value of the dynamical
pairing coordinate $\lambda_{2}$ increases, all three collective masses
decrease and the underlying shell effects gradually vanish. This effect is
consistent with the results shown in Fig. 6 of Ref. Tao  _et al._ (2017).

Figure 3: (Color online) Perturbative cranking masses $\mathcal{M}^{Cp}_{11}$,
$\mathcal{M}^{Cp}_{22}$, and the non-perturbative cranking mass
$\mathcal{M}^{C}_{33}$ (in $\hbar^{2}$ MeV-1) along the static fission path
for several values of $\lambda_{2}$.

The scission contours in the ($\beta_{2},\beta_{3}$) plane are shown in Fig. 4
for several collective pairing coordinates $\lambda_{2}$. The contours are
generally not very different, especially for asymmetric fission. In
particular, scission points that are close to the static fission path are not
sensitive to dynamical pairing. For larger values of $\lambda_{2}$, however,
the scission contour is shifted towards smaller quadrupole deformations
$\beta_{2}$ values for nearly symmetric fission.

Figure 4: (Color online) The scission controur of 228Th in the
($\beta_{2},\beta_{3}$) deformation plane for several values of the collective
pairing coordinate$\lambda_{2}$ .

In Fig. 5 we compare the theoretical predictions for the charge yields with
the data for photo-induced fission of 228Th. The 3D calculation employs the
collective space built from the deformation $\beta_{2}$, $\beta_{3}$ and
pairing $\lambda_{2}$ coordinates, while the 2D calculation includes only the
shape degrees of freedom $\beta_{2}$ and $\beta_{3}$, and static pairing
correlations. As in our study of Ref. Tao  _et al._ (2017), the calculation in
the 2D collective space corresponds to normal and enhanced static pairing,
that is, to 100% and 110% of the normal pairing strength, determined by the
empirical pairing gaps. The theoretical predictions follow the general trend
of the data, except that our model obviously cannot reproduce the odd-even
staggering of the experimental charge yields. The calculation that includes
only the 2D collective space with a static pairing strength adjusted to
empirical ground-state pairing gaps in this mass region (100%), predicts
yields that are entirely dominated by asymmetric fission with peaks at $Z=35$
and $Z=55$. By increasing static pairing (110%), the asymmetric peaks are
reduced and a contribution of symmetric fission develops, but not strong
enough to reproduce the data. It is interesting to notice that a very similar
distribution of charge yields is predicted by the 3D model calculation that
includes dynamical pairing. On a quantitative level, even the 3D calculation
does not completely reproduce the experimental yields. The model predicts
tails of the asymmetric peaks that are not seen in experiment, and thus fails
to quantitatively match the symmetric contribution. It has to be noted,
however, that in the present study the collective potential and inertia tensor
have been calculated at zero temperature. In our recent study of finite
temperature effects in TDGCM+GOA Zhao  _et al._ (2019a), a calculation of
induced fission of 226Th has shown that, although the model can qualitatively
reproduce the empirical triple-humped structure of the fission charge and mass
distributions already at zero temperature, the position of the asymmetric
peaks and the symmetric-fission yield can be described much better when the
potential and collective mass are determined at a temperature that
approximately corresponds to the internal excitation energy of the fissioning
nucleus.

Figure 5: (Color online) Charge yields for induced fission of 228Th,
calculated in the 3D collective space built from the deformation $\beta_{2}$,
$\beta_{3}$ and dynamical pairing $\lambda_{2}$ coordinates (solid red curve).
The yields are shown in comparison to the results obtained in the 2D space of
shape degrees of freedom $\beta_{2}$ and $\beta_{3}$, with static pairing
correlations adjusted to empirical ground-state pairing gaps (100% pairing
strength), and enhanced by ten percent (110% pairing strength). The data for
photo-induced fission correspond to photon energies in the interval 8-14 MeV,
and peak value of $E_{\gamma}=11$ MeV Schmidt  _et al._ (2000).

Finally, to illustrate the effect of dynamical pairing on the flux of the
probability current through the scission hyper-surface, in Fig. 6 we plot the
time-integrated flux through the scission contour in the
$(\beta_{2},\beta_{3})$ plane, for a given value of the pairing collective
coordinate $\lambda_{2}$

$B(\lambda_{2})\propto\sum_{\xi\in\mathcal{B}}{\lim_{t\to\infty}F(\xi,\lambda_{2},t)}.$
(15)

The set $\mathcal{B}(\xi\equiv\beta_{2},\beta_{3})$ contains all elements of
the scission contour with a given value $\lambda_{2}$. Even though it appears
that dynamical pairing does not significantly modify the scission contour (cf.
Fig. 4), nevertheless its effect on the collective flux and, therefore, on the
occurrence of fission, is remarkable. For negative values of $\lambda_{2}$,
that is, for correlations weaker than static pairing at $\lambda_{2}=0$, the
flux rapidly decreases to zero. For positive values of $\lambda_{2}$ the flux
exhibits a steep increase and a prominent peak at $\lambda_{2}\approx 0.3$.
Note that this value corresponds to an increase of $\approx 20\%$ with respect
to the static proton and neutron pairing gaps (cf. Fig. 2). The collective
flux through the scission contour weakens with a further increase of pairing,
and eventually vanishes for $\lambda_{2}>1$.

Figure 6: Time-integrated collective flux $B(\lambda_{2})$ Eq. (15) through
the scission contour, as a function of the pairing collective coordinate
$\lambda_{2}$.

##  IV Summary

The influence of dynamical pairing degrees of freedom on induced fission has
been investigated in a unified theoretical framework based on the generator
coordinate method (GCM) with the Gaussian overlap approximation (GOA). In an
illustrative calculation of fragment charge yields for induced fission of
228Th, the collective potential and inertia tensor have been computed using
the self-consistent multidimensionally constrained relativistic mean field
model, based on the energy density functional DD-PC1, and with pairing
correlations treated in the BCS approximation with a separable pairing force
of finite range. The fission fragment charge distributions are obtained by
propagating the initial collective state in time with the time-dependent
GCM+GOA. The flux of the probability current through the scission hyper-
surface determines the probability of observing a given pair of fragments.

In this work, for the first time, the dynamics of induced fission has been
consistently described in a three-dimensional space of collective coordinates
that, in addition to the axial quadrupole and octupole intrinsic deformations
of the nuclear density, also includes an isoscalar pairing degree of freedom.
A number of studies has already demonstrated the importance of dynamical
pairing for the calculation of spontaneous fission lifetimes, and static
pairing correlations for modeling induced fission. As this work has also
shown, a much more difficult problem is the inclusion of dynamical pairing
degrees of freedom in a time-dependent description of induced fission. One
expects, of course, that the model becomes more realistic as the Hilbert space
of collective coordinates is expanded. However, as the present analysis has
illustrated, it can be difficult to numerically stabilize the time-evolution
of the fissioning system when shape and pairing collective coordinates are
considered simultaneously in a three-dimensional calculation. Here it was
necessary to compute the collective inertia using different cranking
approximations for the shape and pairing degrees of freedom. Nevertheless,
this approach enables a qualitative study of the effect of dynamical pairing
on induced fission.

This initial study has clearly demonstrated the important effect that
dynamical pairing correlations have on the induced-fission fragment
distribution. In particular, the charge distribution of fission fragments of
228Th is characterized by symmetric and asymmetric peaks, but this structure
cannot be reproduced in a two-dimensional calculation that only includes shape
collective variables. In that case, and with a static pairing strength
adjusted to empirical ground-state pairing gaps, the calculated yields
correspond to a completely asymmetric fission. Only by artificially increasing
the static pairing correlations or, more naturally, by including the dynamical
pairing degree of freedom in the three-dimensional calculation, the asymmetric
peaks get reduced and a contribution of symmetric fission develops in
agreement with the empirical trend. It is also interesting to note that the
time-integrated collective flux through the scission contour in the
$(\beta_{2},\beta_{3})$ plane, exhibits a characteristic functional dependence
on the pairing collective coordinate, with a prominent peak at a value that
correspond to an increase of $\approx 20\%$ with respect to the static pairing
gaps.

Future advances in computational capabilities will open the possibility of
more quantitative applications of multi-dimensional TDGCM+GOA to fission
dynamics. An immediate task will be to consider shape and pairing degrees of
freedom on an equal footing, and consistently compute the corresponding
collective inertia tensor and metric. A more challenging problem is to include
dynamical pairing degrees of freedom in recently developed fission models that
attempt to incorporate restoration of symmetries broken by the intrinsic
densities in constrained mean-field calculations (rotational, reflection, and
particle number symmetry) Marević and Schunck (2020); Verriere  _et al._
(2019, 2021).

###### Acknowledgements.

This work has been supported in part by 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). It has also been supported by the National
Natural Science Foundation of China under Grant No. 12005107 and No. 11790325.

## References

  * Bender  _et al._ (2020) M. Bender, R. Bernard, G. Bertsch, S. Chiba, J. Dobaczewski, N. Dubray, S. A. Giuliani, K. Hagino, D. Lacroix, Z. Li, P. Magierski, J. Maruhn, W. Nazarewicz, J. Pei, S. Péru, N. Pillet, J. Randrup, D. Regnier, P.-G. Reinhard, L. M. Robledo, W. Ryssens, J. Sadhukhan, G. Scamps, N. Schunck, C. Simenel, J. Skalski, I. Stetcu, P. Stevenson, S. Umar, M. Verriere, D. Vretenar, M. Warda, and S. Åberg, Journal of Physics G: Nuclear and Particle Physics 47, 113002 (2020).
  * Berger  _et al._ (1991) J. Berger, M. Girod, and D. Gogny, Comput. Phys. Commun. 63, 365 (1991).
  * Verriere and Regnier (2020) M. Verriere and D. Regnier, Front. Phys. 8, 233 (2020).
  * Goutte  _et al._ (2005) H. Goutte, J. F. Berger, P. Casoli, and D. Gogny, Phys. Rev. C 71, 024316 (2005).
  * Regnier  _et al._ (2016a) D. Regnier, M. Verriere, N. Dubray, and N. Schunck, Comput. Phys. Commun. 200, 350 (2016a).
  * Regnier  _et al._ (2018) D. Regnier, N. Dubray, M. Verriere, and N. Schunck, Comput. Phys. Commun. 225, 180 (2018).
  * Regnier  _et al._ (2016b) D. Regnier, N. Dubray, N. Schunck, and M. Verriere, Phys. Rev. C 93, 054611 (2016b).
  * Regnier  _et al._ (2017) D. Regnier, N. Dubray, N. Schunck, and M. Verriere, EPJ Web Conf. 146, 04043 (2017).
  * Verriere  _et al._ (2017) M. Verriere, N. Dubray, N. Schunck, D. Regnier, and P. Dossantos-Uzarralde, EPJ Web Conf. 146, 04034 (2017).
  * Regnier  _et al._ (2019) D. Regnier, N. Dubray, and N. Schunck, Phys. Rev. C 99, 024611 (2019).
  * Zdeb  _et al._ (2017) A. Zdeb, A. Dobrowolski, and M. Warda, Phys. Rev. C 95, 054608 (2017).
  * Vretenar  _et al._ (2005) D. Vretenar, A. Afanasjev, G. Lalazissis, and P. Ring, Phys. Rep. 409, 101 (2005).
  * Meng  _et al._ (2006) J. Meng, H. Toki, S. Zhou, S. Zhang, W. Long, and L. Geng, Prog. Part. Nucl. Phys. 57, 470 (2006).
  * Meng (2016) J. Meng, ed.,  _Relativistic Density Functional for Nuclear Structure_ , International Review of Nuclear Physics, Vol. 10 (World Scientific, 2016).
  * Zhao  _et al._ (2015) J. Zhao, B.-N. Lu, T. Niksic, and D. Vretenar, Phys. Rev. C 92, 064315 (2015).
  * Zhao  _et al._ (2016) J. Zhao, B.-N. Lu, T. Niksic, D. Vretenar, and S.-G. Zhou, Phys. Rev. C 93, 044315 (2016).
  * Tao  _et al._ (2017) H. Tao, J. Zhao, Z. P. Li, T. Niksic, and D. Vretenar, Phys. Rev. C 96, 024319 (2017).
  * Zhao  _et al._ (2019a) J. Zhao, T. Niksic, D. Vretenar, and S.-G. Zhou, Phys. Rev. C 99, 014618 (2019a).
  * Zhao  _et al._ (2019b) J. Zhao, J. Xiang, Z.-P. Li, T. Niksic, D. Vretenar, and S.-G. Zhou, Phys. Rev. C 99, 054613 (2019b).
  * Zhao  _et al._ (2020) J. Zhao, T. Nikšić, D. Vretenar, and S.-G. Zhou, Phys. Rev. C 101, 064605 (2020).
  * Lu  _et al._ (2014) B.-N. Lu, J. Zhao, E.-G. Zhao, and S.-G. Zhou, Phys. Rev. C 89, 014323 (2014).
  * Zhao  _et al._ (2017) J. Zhao, B.-N. Lu, E.-G. Zhao, and S.-G. Zhou, Phys. Rev. C 95, 014320 (2017).
  * Negele  _et al._ (1978) J. W. Negele, S. E. Koonin, P. Möller, J. R. Nix, and A. J. Sierk, Phys. Rev. C 17, 1098 (1978).
  * Łojewski and Staszczak (1999) Z. Łojewski and A. Staszczak, Nucl. Phys. A 657, 134 (1999).
  * Sadhukhan  _et al._ (2014) J. Sadhukhan, J. Dobaczewski, W. Nazarewicz, J. A. Sheikh, and A. Baran, Phys. Rev. C 90, 061304 (2014).
  * Staszczak  _et al._ (1985) A. Staszczak, A. Baran, K. Pomorski, and K. Böning, Phys. Lett. B 161, 227 (1985).
  * Bernard  _et al._ (2019) R. Bernard, S. A. Giuliani, and L. M. Robledo, Phys. Rev. C 99, 064301 (2019).
  * Lazarev (1987) Y. A. Lazarev, Phys. Scr. 35, 255 (1987).
  * Sadhukhan  _et al._ (2016) J. Sadhukhan, W. Nazarewicz, and N. Schunck, Phys. Rev. C 93, 011304 (2016).
  * Bulgac  _et al._ (2016) A. Bulgac, P. Magierski, K. J. Roche, and I. Stetcu, Phys. Rev. Lett. 116, 122504 (2016).
  * Qiang  _et al._ (2021) Y. Qiang, J. C. Pei, and P. D. Stevenson, Phys. Rev. C 103, l031304 (2021).
  * Nikšić  _et al._ (2008) T. Nikšić, D. Vretenar, and P. Ring, Phys. Rev. C 78, 034318 (2008).
  * Tian  _et al._ (2009) Y. Tian, Z. Y. Ma, and P. Ring, Phys. Lett. B 676, 44 (2009).
  * Lu  _et al._ (2012) B.-N. Lu, E.-G. Zhao, and S.-G. Zhou, Phys. Rev. C 85, 011301(R) (2012).
  * Zhou (2016) S.-G. Zhou, Phys. Scr. 91, 063008 (2016).
  * Vaquero  _et al._ (2011) N. L. Vaquero, T. R. Rodriguez, and J. L. Egido, Phys. Lett. B 704, 520 (2011).
  * Vaquero  _et al._ (2013) N. L. Vaquero, J. L. Egido, and T. R. Rodrǵuez, Phys. Rev. C 88, 064311 (2013).
  * Gambhir  _et al._ (1990) Y. Gambhir, P. Ring, and A. Thimet, Ann. Phys. 198, 132 (1990).
  * Baran  _et al._ (2011) A. Baran, J. A. Sheikh, J. Dobaczewski, W. Nazarewicz, and A. Staszczak, Phys. Rev. C 84, 054321 (2011).
  * Yuldashbaeva  _et al._ (1999) E. Yuldashbaeva, J. Libert, P. Quentin, and M. Girod, Phys. Lett. B 461, 1 (1999).
  * Brack  _et al._ (1972) M. Brack, J. Damgaard, A. S. Jensen, H. C. Pauli, V. M. Strutinsky, and C. Y. Wong, Rev. Mod. Phys. 44, 320 (1972).
  * Nilsson  _et al._ (1969) S. G. Nilsson, C. F. Tsang, A. Sobiczewski, Z. Szymanski, S. Wycech, C. Gustafson, I.-L. Lamm, P. Moller, and B. Nilsson, Nucl. Phys. A 131, 1 (1969).
  * Girod and Grammaticos (1979) M. Girod and B. Grammaticos, Nucl. Phys. A 330, 40 (1979).
  * Bés and Szymański (1961) D. Bés and Z. Szymański, Nucl. Phys. 28, 42 (1961).
  * Sobiczewski  _et al._ (1969) A. Sobiczewski, Z. Szymański, S. Wycech, S. G. Nilsson, J. R. Nix, C. F. Tsang, C. Gustafson, P. Möller, and B. Nilsson, Nucl. Phys. A 131, 67 (1969).
  * Younes and Gogny (2009) W. Younes and D. Gogny, Phys. Rev. C 80, 054313 (2009).
  * Schmidt  _et al._ (2000) K.-H. Schmidt, S. Steinhauser, C. Bockstiegel, A. Grewe, A. Heinz, A. Junghans, J. Benlliure, H.-G. Clerc, M. de Jong, J. Muller, M. Pfutzner, and B. Voss, Nucl. Phys. A 665, 221 (2000).
  * Marević and Schunck (2020) P. Marević and N. Schunck, Phys. Rev. Lett. 125, 102504 (2020).
  * Verriere  _et al._ (2019) M. Verriere, N. Schunck, and T. Kawano, Phys. Rev. C 100, 024612 (2019).
  * Verriere  _et al._ (2021) M. Verriere, N. Schunck, and D. Regnier, Phys. Rev. C 103, 054602 (2021).