Enhanced Generator Coordinate Method: eGCM

The early attempts to describe molecular spectra by W. Heitler, F. London, J. C. Slater, L. Pauling, and others used hybridized localized atomic orbitals, centered on atoms. Later F. Hund, R. S. Mulliken, and J. Lennard-Jones introduced delocalized molecular orbitals, which proved more flexible in practice. The atomic orbitals centered at atom positions describe electrons in various states, e.g. the $\sigma-$ and $\pi-$orbitals used to describe the valence electrons in a carbon atom. At each atomic site this set of orbitals is not complete and higher energy excitations of the system are thus not included. In the case of carbon atoms the excitation of the deep lying electrons with principal quantum numbers $n=1$ and $l=0$ or to levels with $n\geq 3$ are not allowed. Whether this is a good approximation naturally depends on the specific problem considered, which typically encompasses relatively low energy local excitations of the quantum system under consideration. These ideas propagated further in condensed matter theory in order to describe the band structure of solids Aschroft and Mermin (1976) and lately in the high-temperature superconductivity in the Hubbard model Dong et al. (2022). Localized or delocalized orbitals describe electrons in various quantum states, including the spin degrees of freedom.

These ideas also inspired J. A. Wheeler and his students D. A. Hill and J. J. Griffin Hill and Wheeler (1953); Griffin and Wheeler (1957) to introduce a description of nuclear LACM, nuclear fission in particular, for which they coined the term GCM. These authors considered however a simplified version of the molecular and condensed matter frameworks, by using at “each site,” in that case for a particular nuclear shape typically characterized by a quadrupole $Q_{20}$ and an octupole deformation $Q_{30}$, only the ground states of a nucleus with that shape, and the wave function of the nucleus was represented as a linear superposition of these fixed shape ground states. In nuclear literature the GCM was considered over the years in two flavors, with time-independent “generator wave function” function $f(Q)$ and with a time-dependent “generator wave function” $f(Q,t)$ respectively

but in both cases with static (generalized) Slater determinants $\Phi(\xi_{1},\ldots,\xi_{A}|Q)$ Ring and Schuck (2004); Verriere and Regnier (2020); Schunck and Regnier (2022); Schunck (2022), as solutions of chosen by “educated guesses” of a restricted set of nuclear shape constraint Hartree-Fock (HF) or Hartree-Fock-Bogoliubov (HFB) solutions. These two choices are neither well nor uniquely defined and consequently the “sum” over “nuclear shapes” parametrized by the multidimensional variable $Q$ lead to a uncontrolled approximation of the many-fermion wave function $\Psi(\xi_{1},\ldots,\xi_{A})$ or $\Psi(\xi_{1},\ldots,\xi_{A},t)$, see Ref. sup . A correct expansion in terms of either static or time-dependent (generalized) Slater determinants should include all possible many-fermion excitations at a “given shape,” not only the “ground state” at a given “shape.” At the same time by including all excited states at all “nuclear shapes” will make the entire set overcomplete. In the literature, as far as I know, there is no study of the accuracy of the GCM or of the accuracy of the similar framework based on molecular orbitals. At best, studies demonstrate that there is “agreement” with a reduced set of observables, when a reasonable basis set of wave functions, chosen by educated guesses, is used.

Over the many decades of GCM use in nuclear physics, the belief, or rather the hope, has been that the total wave functions will follow an adiabatic evolution from one instantaneous “ground state shape” to the next neighboring instantaneous “ground state shape,” which led to the development of the adiabatic time-dependent Hartree-Fock theory Ring and Schuck (2004). In chemistry it has been known for decades that this naive Born-Oppenheimer approximation Born and Oppenheimer (1927) is violated Tully (1990). There are also other possible theoretical frameworks, see Refs. Requist and Gross (2016); Agostini and Gross (2020); Bulgac et al. (2019a); Bulgac (2022) and earlier references therein. In the case of nuclear LACM the nuclear shape has to evolve typically in low energy nuclear dynamics in such a manner that the local Fermi momentum distribution remains approximately spherical Bertsch (1980); Barranco et al. (1988, 1990); Bertsch and Flocard (1991); Bertsch and Bulgac (1997), as otherwise the volume contribution to the energy of the nucleus will change. Nuclear systems are saturating many fermion systems and in the low energy dynamics the contribution to the volume energy and the average density should remain essentially constant Meitner and Frisch (1939); Hill and Wheeler (1953); Bertsch (1980). Only the Coulomb and surface isoscalar and isovector contributions to the total nucleus energy can vary considerably, in agreement with the brilliant insight Meitner and Frisch (1939) had in the case of nuclear fission. The energy contribution due to the emergence of pairing correlations is always a small contribution, but however crucial. In the Bethe-Weiszäcker mass formula the odd-even correction term is about $3$ times smaller than the root mean square error of the binding energy of any nucleus with an atomic mass larger than $A\approx 50$. At a well defined “nuclear shape” a many nucleon system has a very rich spectrum of excited states, qualitatively similar to the various states of electron systems, or similar to the full spectrum of a nucleus described in a large shell model basis. While changing its shape a nucleus does not have to hop from one ground state of a given shape to another ground state of a different shape only. This “conspicuous” deficiency was attempted to be “fixed” in Ref. Müther et al. (1977); Bernard et al. (2011); Chen and Egido (2017), in an approach where only a patently insufficient number of excited states was taken into account. As far as I am aware this (uncontrollable) extension of GCM Bernard et al. (2011) was never shown to lead to a satisfactory description of LACM, and moreover, it is not clear whether such an extension is numerically feasible or practically meaningful if all relevant quasiparticle excitations are taken into account.

It is well established theoretically that LACM of a fissioning nucleus beyond the outer fission barrier is a strongly dissipative non-equilibrium process Bulgac et al. (2016, 2019b, 2020); Bender and et al. (2020) and the adiabaticity invoked by Wheeler and collaborators Hill and Wheeler (1953); Griffin and Wheeler (1957) is an unphysical assumption, which however is still widely used in GCM “microscopic” approaches to nuclear fission up to the scission configuration Verriere and Regnier (2020); Schunck and Regnier (2022); Schunck (2022). In theoretical studies of LACM in the first and second potential well, and in particular in the case of spontaneous fission Balantekin and Takigawa (1998); Sadhukahn (2020) the attitude is that the collective motion is adiabatic. The argument brought forward in the case of spontaneous fission, which is to a large extent an under the barrier penetration process, is that due to pairing effects the adiabatic approximation is valid. This goes against the solid theoretical arguments presented by Caldeira and Leggett (1983) and widely accepted in condensed matter physics, that the coupling to internal excitations (in nuclear language parlance, coupling to excitations above “the instantaneous ground shape” and when "polarization effects" are neglected) leads to longer tunneling times. This is in qualitative agreement with experimental observations of the spontaneous life-times of odd-mass and odd-odd nuclei Vandenbosch and Huizenga (1973), and also with our recent findings Bulgac (2024), where in particular we demonstrate that the Pauli blocking approximation is not valid in non-equilibrium processes and therefore Cooper pairs are broken during LACM.

In 1983 Reinhard et al. (1983) suggested to replace the static fixed set of ground states a nucleus evolves through during a LACM with the solutions of a time-dependent HF problem. In this framework the total time dependent nuclear wave function of the nucleus acquires a more complex structure, it is in general an infinite sum/integral over many (generalized) time-dependent Slater determinants $\Phi(\xi_{1},\ldots,\xi_{A}|Q,t)$. Allegedly, this prescription may describe a dissipative non-equilibrium process such as nuclear fission or heavy-ion reaction. In Ref. Reinhard et al. (1983) it was implicitly assumed that various TDHF trajectories $\Phi(\xi_{1},\ldots,\xi_{A}|Q,t)$ are started simultaneously and span a sufficiently large set of initial nuclear shapes, described by the shape (multidimensional) parameter $Q$. The “generator wave function” $f(Q,t)$ is a solution of the time-dependent Hill-Wheeler equation

where the nucleon coordinates $\xi_{1},\ldots,\xi_{A}$ have been suppressed and the matrix elements ${\cal N}(Q,Q^{\prime})=\langle\Phi(Q,t)|\Phi(Q^{\prime},t)\rangle$ and $\langle\Phi(Q,t)|H-H_{MF}|\Phi(Q^{\prime},t)\rangle$ are evaluated by integrating over the coordinates $\xi_{1},\ldots,\xi_{A}$. Here $H$ and $H_{MF}$ stand for the many-body and mean field Hamiltonians, and $i\hbar\partial_{t}\Phi(Q,t)=H_{MF}\Phi(Q,t)$.

In order to understand the incompleteness of the arguments presented by Reinhard et al. (1983) I have to make a detour to introducing in quantum many-body theory the collective and intrinsic degrees of freedom, a process described a long time ago by Feynman and Vernon (1963), from which either a classical  Grangé et al. (1983); Weidenmüller and Zhang (1984) or quantum Fokker-Planck equation Bulgac et al. (1996, 1998); Kusnezov et al. (1999) for the “collective degrees of freedom” can be derived. In the GCM framework one typically interprets the (generalized) Slater determinants labels $Q$ as a set of collective variables, thus practically introducing a poor man’s Feynman-Vernon separation between the “collective degrees of freedom $Q$”, which need to be requantized, and the "intrinsic degrees of freedom $\xi_{1,\ldots,A}$." In GCM this is achieved by adopting the Gaussian overlap approximation (GOA) Ring and Schuck (2004); Verriere and Regnier (2020); Schunck and Regnier (2022); Schunck (2022) of the norm and Hamiltonian overlaps and “deriving” a Schrödinger-like equation for the “collective degrees of freedom $Q$.” In the most theoretically advanced implementation of the TDDFT extended to superfluid systems Bulgac et al. (2016, 2019b, 2020), as in any previous time-dependent mean field approach the “collective degrees of freedom $Q$” are simply labels for the initial conditions in case of time-dependent problems. Mixing different (TD)DFT trajectories $Q$ requires knowledge of the function either $f(Q)$ or $f(Q,t)$, which can be obtained only by solving the corresponding static or time-dependent GCM equations or using other phenomenological or microscopic approaches Wilkins et al. (1976); Brosa et al. (1990); Lemaître et al. (2015); Ishizuka et al. (2017); Sierk (2017); Albertsson et al. (2020); Ivanyuk et al. (2024); Sadhukhan et al. (2016, 2017); Sadhukhan (2020).

A very instructive example of a GCM-like application, though it was never considered as one, is the treatment of pairing correlations, and specifically the treatment of a bound electron pair by Cooper (1956). The pairing Hamiltonian in second quantization is

and the electron pair wave function is $|\phi\rangle=\sum_{k}C_{k}a^{\dagger}_{k}a^{\dagger}_{\overline{k}}|0\rangle$, with a sum over two-fermion Hartree-Fock states in time reversed states $k$ and $\overline{k}$, and here the coefficients $C_{k}$ are obtained as solution of a GCM-like equation, in which the norm overlap is diagonal, while the Hamiltonian overlap is not and leads to hops between different “generator coordinate” states $a^{\dagger}_{k}a^{\dagger}_{\overline{k}}|0\rangle$, where $|0\rangle$ is the vacuum state. In GCM calculations the basis states are not typically orthogonal to each other and the norm overlap as a result is non-diagonal, which leads to the usual technical problems. One can extend trivially Cooper (1956)’s treatment of a single pair to several pairs, but when the number of pairs $N\geq 5$ one should switch to the Hartree-Fock-Bogoliubov (HFB) approximation.

There were a few attempts in recent literature Regnier and Lacroix (2019); Hasegawa et al. (2020); Marević et al. (2023, 2024); Li et al. (2023, 2024) to implement the framework suggested by Reinhard et al. (1983), but the results obtained so far have not shown how to obtain the expected mass and charge distributions of fission fragments (FFs) or how to describe an expected effective mixing between different trajectories Li et al. (2024), or with a better accuracy than in a quite wide variety of other approaches Wilkins et al. (1976); Brosa et al. (1990); Lemaître et al. (2015); Ishizuka et al. (2017); Sierk (2017); Albertsson et al. (2020); Ivanyuk et al. (2024); Sadhukhan et al. (2016, 2017); Sadhukhan (2020).

There are several reasons why the framework outlined by Reinhard et al. (1983) is not going to succeed in describing the non-equilibrium dynamics of a many-fermion systems, since the initial conditions are not well defined. As shown in Refs. Bulgac et al. (2019b, 2020), fission trajectories started along the rim of the outer fission barrier, after a very short time appear to focus into a rather narrow bundle. Unfortunately these trajectories reach a particular set of FF separations at different times, as the saddle-to-scission times vary quite a bit, depending on the initial values of $Q$, and as a result the norm and Hamiltonian overlaps evaluated of a given time $t$, counted from the time each specific trajectory was started, leads to practically vanishing values of the Hamiltonian overlaps and therefore to no “expected quantum mixing”’ between different trajectories parametrized by the “classical variables $Q$” Li et al. (2024). This is very similar to the Thomas Young two-slit or the Mach-Zehnder experiments, with either light waves or quantum objects, when particles hit the screen and interfere after traveling for different times along different trajectories and hit a spot on the screen after different times spent from the emission. This is the main element overlooked by Reinhard et al. (1983). In induced fission an incident low energy neutron beam excites the system in its ground state potential well and different “classical trajectories,” described by the “collective variables $Q$,” reach the rim of the outer fission barrier at different times and eventually the FF detector, the “screen,” where the interference occurs.

The solution to this apparently complicated problem is rather simple, see also the supplement sup . There a number of relevant references are listed Peierls and Yoccoz (1957); Peierls and Thouless (1962); Ryssens et al. (2015); Bulgac (2013); Bulgac et al. (2018, 2021, 2022); Scamps et al. (2023); Zagrebaev et al. (2013); Sekizawa and Yabana (2013); Berry (1984); Bulgac (1988, 2024); Barrett et al. (2013); Scharnweber et al. (1971); Bulgac (1989); Taylor (1909); Merli et al. (1976); Arndt et al. (1999); Hanbury Brown and Twiss (1956); Boal et al. (1990); Bell (1987), a relation to the old works of Peierls and Yoccoz (1957) and Peierls and Thouless (1962) is presented and discussed how very strong correlations, absent in mean field, are generated. Instead of the norm and Hamiltonian overlap matrix elements in Eq. (3) one should evaluate the new generalized norm and the Hamiltonian overlap matrix elements and solve the extended GCM (eGCM) equations

using the variational approach with the eGCM many-body wave function $\sumint_{Q,t}|\Phi(Q,t)\rangle f(Q,t)$, which is akin to Feynman path integral representation sup . One can introduce the separation between FFs along each TDDFT trajectory Bulgac (2019); Bulgac et al. (2019b); Bulgac (2013) as a function of time $s(t)$. The new norm ${\cal N}(Q,s|Q^{\prime},s^{\prime})$ and the Hamiltonian ${\cal H}(Q,s|Q^{\prime},s^{\prime})$ overlaps are non-vanishing if $s(t)$ is similar in magnitude to $s(t^{\prime})$, in the case of fission of an axially symmetric nucleus, when the FFs are in close proximity of each other sup . Note that TDDFT is used to generate a physically sound choice set of many-body wave functions $|\Phi(Q,t)\rangle$ (analogous to the set of harmonic single-particle wave functions used in shell model calculations, and which is desirably complete sup for a controlled approximation) and they do not have to be related with the Hamiltonian $H$ used in Eq. (Enhanced Generator Coordinate Method: eGCM), unlike in (TD)GCM applications Verriere and Regnier (2020); Schunck and Regnier (2022); Schunck (2022). eGCM is capable of describing interference patterns in the detector, even though particles arriving by different paths and at different times $t\neq t^{\prime}$, a fact known otherwise for a very long time sup ; Taylor (1909); Merli et al. (1976); Arndt et al. (1999).

It is instructive to develop an intuition concerning the meaning and importance of the parameters $(Q,t)$ and important to discuss the feasibility of implementing eGCM. In an analogy with with a non-uniform diffraction grating with a finite number of slits, one can associate $Q$ with the label of a specific slit, in this case including not only the $(Q_{20},Q_{30})$ position of the initial point on the rim of the outer barrier, but if desired also the Euler angles specifying of the fissioning nucleus needed to perform an angular momentum projection.

In the case of generalized Slater determinants used in the case of induced fission one has to include all allowed quasiparticle states in constructing a time-dependent “trajectory,” see Ref.  Bulgac et al. (2024a), which for a typical situation in a box $30^{2}\times 60$ fm³ and a lattice constant $l=1$, which is a very large number $4\times 30^{2}\times 60$ = 216,000, unlike in all known to us real-time numerical simulations of fission used by other authors in literature so far. Fortunately, there is a solution to this problem, at each specific set of values $(Q,t)$ one should introduce the canonical wave functions and in that case the size of the needed number of quasiparticle wave functions dramatically drops to a several hundreds, see Fig. 1 and Refs. Scamps et al. (2023); Bulgac et al. (2024a), and the needed norm and Hamiltonian overlaps can be evaluated rather rapidly.

I will assume that one performs induced fission simulations in a typical ($64$ fm)³ box, if angular momentum projection is considered. Otherwise a $32^{2}\times 64$ fm³ simulation box is appropriate. For evaluating the first and second spatial derivatives the use of FFT, which leads to machine precision and using powers of 2 is the best choice for spatial dimensions Jin et al. (202); Bulgac et al. (2011). The lattice constant $l=1$ fm corresponds to a maximum momentum cutoff in one cartesian direction $p_{max}=\hbar\pi/l\approx$ 600 MeV/c, which is of the order of the maximum momentum cutoff considered in chiral effective field theory of nucleon interactions. It would be sufficient to use a number $N_{Q}\leq 15$ for the set of quadrupole and octupole deformations $(Q_{20},Q_{30})$ for an axially symmetric even-even compound nucleus along the rim of the outer fission barrier, as was done in Refs. Bulgac et al. (2016, 2019b, 2020). Sometimes it maybe convenient to parameterize a given trajectory in terms of the separation between FFs, starting with a separation $s\approx 7$ fm, when the neck is emerging, until $s\approx 22$ fm, when the FF shapes are relaxed and spatially rather well separated, and use the relation between the time $t$ along the trajectory and the FFs separation $s(t)$, where $t$ is the running time along a given trajectory. Retaining for the eGCM a $\Delta s=0.25$ fm one would need $N_{s}=60$ different FFs separations along a given trajectory. As it was amply demonstrated recently, see again Fig. 1 and Refs. Scamps et al. (2023); Bulgac et al. (2024b, c), it is sufficient to use not more than hopefully $N_{\rm{cwfs}}<200,\ldots,300$ canonical wave functions with the highest occupation probabilities, for the proton and neutron subsystems respectively.

If one might try instead to eschew the determination of the canonical wave functions and use instead quasiparticle wave functions with the highest occupation probabilities at each FF separation $s$, or use a static or time-dependent Bardeen-Cooper-Schrieffer approximation, the errors are large, essentially irrespective of how many states are included. A further simplification of the eGCM framework can be achieved by adopting an appropriate generalized Gaussian overlap approximation (GOA) for the norm kernel ${\cal N}(Q,t|Q^{\prime},t^{\prime})$ and derive a Schrödinger-like equation for $f(Q,t)$. The GOA is however a further approximation, prone to large errors Verrière et al. (2017), and it was prompted only to find a simple simulacrum to a “collective Hamiltonian quadratic in collective coordinates” to the Feynman-Vernon influence functional, which emerges from a path integral formulation of the many-body problem. As it was shown in Ref. Bulgac (1988), in the case of nonlocal equations a reduction to a second order partial differential equation is not always possible, as the system has a behavior similar to a birefringent medium.

In the case of fissioning even-even nuclei with axial symmetry one can build a rather simple angular momentum projection procedure of the nucleus as only the Euler angles $\gamma$ and $\beta$ are needed and one can use the icosahedral group with only 60 different angles to project on the total spin $J=0$ and 2, using the positions of the 60 carbon atoms in a $C_{60}$ buckyball for different orientations of the “collective” variables $(Q,Q^{\prime})$. This approach is orders of magnitude more economic that the usual angular momentum projecting techniques Ring and Schuck (2004); Lu et al. (2022); Scamps et al. (2023). This leads to a conservative estimate of the dimension of the norm overlap matrices (separately for protons and neutrons) $N_{eGCM}=N_{\rm{angles}}\times N_{\rm{cwfs}}\times N_{s}\times N_{Q}=60\times 200\times 60\times 15=7,200,000$, which is a number significantly smaller than the largest dimensions of the shell-model calculations Johnson (2018); Zbikowski and Johnson (2023). If an angular momentum projection is not performed the dimension of norm overlap matrices is significantly smaller $N_{eGCM}=N_{\rm{cwfs}}\times N_{s}\times N_{Q}=200\times 60\times 15=180,000$. Since only the largest eigenvalues of the norm overlap matrices are needed, full diagonalization is not necessary. As it was recently shown in Ref. Bulgac et al. (2024b), the differences between the number-projected and number-unprojected number densities as a function of time at any point in space are at most at the level 1$\%$ in the case of induced fission treated in TDDFT with pairing correlations included, during the entire time-evolution from the top of the outer fission barrier to complete fission fragments (FFs) separation. Naturally, this new eGCM framework is equally applicable to other cases of nuclear LACM, in particular to collisions of heavy ions. Since the emergence of powerful supercomputers during the last two decades or so, the numerical implementation of eGCM appears to be doable with many existing and rather modest computer platforms.

The most contentious and difficult issue in using GCM in nuclear LACM is the fact that, in DFT and in TDDFT in particular, a true Hamiltonian does not exist and the evaluation of the Hamiltonian overlap within GCM is not well defined. The Hamiltonian overlaps ${\cal H}(Q,t|Q^{\prime},t^{\prime})$ in Eq. (Enhanced Generator Coordinate Method: eGCM) should and can be be evaluated using a true many-body Hamiltonian with two- and three-body interactions, when the evaluation of the norm and Hamiltonian overlaps with generalized mean field many-body wave functions is a well defined procedure Ring and Schuck (2004); Bertsch and Robledo (2012); Carlsson and Rotureau (2021), which has no ambiguity and eGCM can be used with generator wave functions from TDDFT as outlined above and in the supplement sup .

In conclusion, I have outlined an extended GCM framework, dubbed eGCM, which is free of a number of undefined steps in previous formulations, which is expected to correctly describe interference and entanglement between different “classical” fission trajectories, parametrized by the initial shapes $(Q_{20},Q_{30})$ on the rim of the outer fission barrier and different Euler angles $\gamma$ and $\beta$. Since rather soon after scission many TDDFT fission trajectories follow a very similar path Bulgac et al. (2019b, 2020), different trajectories originating at different initial points $Q$ will find themselves in relatively close proximity of each other, albeit each parametrized with a different “time” $t$, which translates into to an actual FFs separation $s(t)$. Having many different trajectories in close proximity of each other will result into a strong mixing, a fact which seemed almost impossible to realize until now Li et al. (2024). Since LACM dissipation is now organically incorporated in TDDFT extended to pairing correlations, there seem to be no other limitation to a full microscopic approach to fission. The only nagging element is that by construction GCM and equally the method of delocalized orbitals in chemistry and condensed matter physics are still uncontrollable approximations, even though the framework is based on sound physical assumptions.

Acknowledgements

I thank I. Stetcu, I. Abdurrahman, and M. Kafker for discussions, L. Troy, I. Abdurrahman, and M. Kafker for reading and suggesting a number of improvements, I. Abdurrahman for preparing the data for Fig. 1, and also D. Vretenar for raising number of questions on an earlier version of the manuscript. The funding from the Office of Science, Grant No. DE-FG02-97ER41014 and also the partial support provided by NNSA cooperative Agreement DE-NA0003841 is greatly appreciated. This research used resources of the Oak Ridge Leadership Computing Facility, which is a U.S. DOE Office of Science User Facility supported under Contract No. DE-AC05-00OR22725.