Lastly, at high to relativistic energies which can be referred to as the participant-spectator regime, the short-ranged two-body interaction overshadows the mean-field contribution and dictates the main evolution. For each of those energy ranges, there are different dedicated models that are welladapted to the specific situation. For instance, the low energy mechanisms are efficiently portrayed within the time-dependent Hartree-Fock (TDHF) framework, while the large-amplitude fluctuations at Fermi energies are more suitably described in stochastic mean-field approaches (Boltzmann-Langevin equation) or in molecular dynamics models. It has been a longstanding challenge for nuclear many-body theories to properly address the competing types of instabilities at the threshold between Fermi and low energies within a unified picture. At the state-of-the-art, models at Fermi energies tend to lack the mean-field phenomenology to some extent, from isospin drifts to collective modes, and, conversely, models approaching low energy miss large-amplitude fluctuation and lack mechanisms where profound transformation, or even splits, occur in the system. Ultimately, a theoretical approach able to cover the transition from Fermi to low energies in one single comprehensive description, would be a well suited framework to study the evolution of fragments and clusters as a function of time and density.

We are adapting a stochastic TDHF formulation by starting from the many-body Schroedinger equation, since this is the fundermental physics which determines the collective behaviour of any quantum system at these energies. The main idea of this model is the decomposition of the non-local nucleonic wave functions into a set of moving basis functions. As a consequence, the mean-field properties are preserved since the nucleonic wave function is not constrained to be localized and furthermore the system is prepared into a set of moving basis functions in order to follow a similar scheme as in the analogous semi-classical Boltzmann-Langevin approach.

I dedicate my dissertation work to my family and many friends. A special feeling of gratitude to my loving parents, Nguyen Thi Lan Phuong and Dinh Van Duc whose words of encouragement have lead me to who I am today. My brother Dinh Trong Tin who has always believed in me and kept me going with his optimism. I also dedicate this dissertation to my close friends who helped me get through the difficult times during the pandemic.

First of all, I would like to thank my supervisor, Paolo Napolitani, for not only his academic guidance, but also for his kindness and profound knowledge about history and philosophy which has helped me to widen my perspective on life. I am grateful for his overwhelming support since the very beginning when we first met in Orsay. Particularly during the difficult times of the confinement due to the COVID-19 pandemic, we stayed frequently in contact which has eased my nerves and lightened me up.

He stayed patient when I was trying to solve a problem but also guided me in the right direction when I got stuck.

Further, I greatly appreciate all the time when I was invited to dinner at home with his family and friends, especially the time I was invited to Christmas which I will never forget. I also highly value that he has provided me with many opportunities to travel and to attend conferences, and that he would go the extra mile to connect me to people. On that note, I would like to thank Maria Colonna and Stefan Typel for fruitful discussions and suggestion and for their time and efforts in writing the thesis report for me. In particular, I had the honor to meet Maria twice in Catania and I greatly appreciate the hospitality that she provided. In addition, it was a pleasure to have a detailed discussion with her where she asked many critical questions. Her vast expertise about various types of nuclear transport models and her clear illustration of her thought process in tackling a problems are an inspiration to me.

I am also grateful to my other committee members, Michael Urban and Abdou Chbihi for accepting my invitation. In particular, I want to thank Michael for taking his time for discussions and advice throughout the years since it was easier for me to open up to someone in the lab who also spoke german.

I shall take this opportunity to send my gratitude and best wishes to my colleagues, specifically my good friends Melih Ozcelik, Carlo Flore, Liza Yedelkina and Stefano Burrello.

i Chapter 1 Introduction

Scope and context

Heavy-ion collisions (HIC) are our best terrestrial laboratory probe for the nuclear equation of state (EOS) [Baldo2012, Baran2005]. The EOS is such a fundamental object that it does not only characterize finite nuclei, but also a large variety of astrophysical scenarios where dense matter is involved [Steiner2005, Oertel2017, Burgio2021, Fantina2022]. Heavy-ion collisions affect the density landscape of a nuclear system [Ring1980], from collective modes to disordered perturbations. In particular, a broad range of situations arise, where neutrons and protons oscillate in phase or out of phase. In both cases, the outcome could often be a profound transformation of the system structure.

Sporadically, since in nuclear matter the nuclear interaction leads to clusterization, the exit channel could also be a catastrophic process, where the whole system disassembles into nuclear clusters and fragments [Chomaz2004]. In contrast to astrophysical conditions, this correlation between heavy-ion collisions and the EOS is highly indirect because of two main reasons. First of all, especially in the early stage of the process, timescales are comparable with the relaxation time of the nuclear interaction so that chaotic behaviour and volume instabilities may arise. Secondly, when moving from nuclear matter conditions to finite systems, also the surface of the system comes into play imposing non-trivial geometries and surface instabilities. To handle such complex phenomenology in a theoretical framework requires to combine one-body (mean field) and beyond-one-body features which have so far been largely approximated in state-of-the-art models or kept uncorrelated. These are the description of the non-local character of wave function and quantum fluctuations deriving from their interference (one body), or the large-amplitude fluctuations which are triggered by dissipative effects (beyond one body). These aspects of the nuclear process are translated into experimental observables wich are becoming accessible in upcoming measurements [Sorensen2023], such as particle-particle correlations and the chronology of particle emission.

Along this line, these experiments on heavy-ion collisions are expected to bring a completely new microscopic view on nuclear dynamics, based on particle-particle correlations and chronology-related observables. This is thanks to the higher granularity of new detectors, higher kinematic resolution and improved identification power. As a consequence, it will be possible to extract the equation of state of nuclear matter from heavy-ion collisions profiting from more microscopic observables, different to previous studies, where thermodynamical observables had to be privileged. As well, high resolution correlation observables will allow one to gain more insight on the process of formation of nuclear fragments and nuclear clusters, rather than the study of the final-state properties. From this prospective, density-functional theories and microscopic models will be essential both to prepare new experiments and to impose even more strict constraints in connecting theory and observables.

For the scope of this thesis, the focus is aimed towards the theoretical description of the dynamical process of heavy-ion collisions, specifically at low to Fermi energies, typically in the range of a few MeV to hundreds of MeV per nucleon, with the purpose of constructing a well suited new microscopic modeling framework. The interest of addressing nuclear reactions in this specific energy range arises because they are characterized by a wide variety of phenomena from collective modes to dissipative processes due to the complex interplay between nuclear structure, nucleon-nucleon interactions, and many-body effects.

From early experiments and models of heavy-ion collisions to nuclear dynamics: a brief history

The study of heavy ion collisions has a long and fascinating history that dates back to the early 20th century [Weinberg1984]. In 1919, Ernest Rutherford performed the first nuclear reaction experiment by bombarding nitrogen gas with alpha particles, which resulted in the production of oxygen and protons [Rutherford1919]. This discovery paved the way for the study of nuclear structure and reactions.

In the following years, one of the first major discoveries in this field was the observation of nuclear fission by Lise Meitner, Fritz Strassmann and Otto Hahn in 1938 [Meitner1938, Meitner1939]. In the 1930s experimental techniques for studying induced nuclear reactions started to be developed and researchers began to explore the use of accelerators to produce high-energy beams for nuclear physics research [Bohr1936]. The 1950s saw a significant advancement in the study of heavy-ion collisions, as researchers began to use particle accelerators to produce high-energy beams of heavy ions. Initially, with this novel technique, heavy nuclei could be bombarded with light particles (neutrons, protons, deuterons) and exit channels combining heavy residues and several light particles could be observed for the first time. This finding led to the earliest and one of the most popular description of a nuclear reaction [Serber1947], where a fast stage of excitation is followed by a second stage of statistical decay [Weisskopf1937] from a thermalized compound nucleus [Thomas1968].

Starting from the 1970s, research in heavy-ion collisions continued to grow, as scientists got more and more interested into the dynamics of nuclear reactions leading to perturbations of a nuclear system. One of the most significant developments in this field was the study of deep-inelastic collisions [Volkov1978, Beck1996, Sanders1999], which provided important insights into the structure of atomic nuclei. This process, as well as fission, inspired the earliest models to describe fragment formation in the spirit of the liquid-drop model and the saddle-point interpretation [Frolich1973, Griffin1976, Brosa1990, Businaro55a, Businaro55b]. Later on, also intermediate energies [Nervik1954] and relativistic energies [RobbGrover1962] could be explored, leading to the discovery that a nuclear system under violent perturbation can also disintegrate into intermediate-mass fragments (IMF), generally lighter than fission products and heavier than light charged particles. This finding led to intense pioneering investigations carried on throughout the 1980s and 1990s on IMF emission from highly excited nuclear sources [Kaufman76, Warwick82, Jakobsson1982, Hirsch84, Andronenko86, Barz86, Korteling90, Kotov95, Hsi97, Avdeyev98, Viola2006]. In the same period, pioneering experiments with heavy-ion beams from Fermi-energy to relativistic domains led to the the discovery that heavy-ion collisions could result into a multifragmentation process in many IMF [Hufner84, Cole2000, Chomaz2006].

Along this line, more and more dissipative conditions could be accessed using increasingly sophisticated experimental techniques, especially around Fermi energies, paving the way to the study of the nuclear EOS. Even though the earliest theoretical approaches were a sort of saddle-point interpretation extended to highly excited systems [Moretto1992, Moretto1988, Lopez94], later on, new descriptions were gradually involving many-body theories and searching analogies with phase transitions [Schmidt2001, Hock2009, Kim2004, Steinheimer2012, Moretto2011]. These interpretations led to a flourishing exploration of statistical approaches and hydrodynamic descriptions [Botvina85b, Botvina90, Bondorf95, Raduta1997, Gross2001]. Expecting that a connection could be established between the nuclear matter equation of state and the arising of signals of liquid-gas phase transitions [Jacquaman1983, Muller1995], two decades of experiments on heavy-ion collisions were addressed to tracking signals of caloric curves [Bowman1991, Pochodzalla1995, Dagostino2000, Schmidt2002, Natowitz2002, Souza2006, Pichon2006, Borderie2008, Bonnet2009, Lehaut2009, Moretto2011, Zheng2012, Borderie2013]. In parallel, an intense study of the instabilities which lead to the phase transition phe-nomenology were carried on, ranging from thermodynamic theories of nuclear matter and finite nuclei [Binder1984, Gulminelli1999, Chomaz2000, Chomaz2002, Gulminelli2003, Chomaz2003, Gross2011] to density-fuctional theories condensed into transport models [Chomaz2004, Baran2005, Colonna2017, Zhang2018, Colonna2021]. The success of transport models in investigating the EOS led to the search of more and more microscopic observables related to mean-field features, like neutron and proton currents and beyond mean-field features, like dissipative regimes and nuclear cluster production. Moreover, since transport models do not depend on equilibrium assumptions, they led to a more comprehensive study of nuclear dynamics, in presence of instabilities, quantum fluctuations and large-amplitude fluctuations related to nucleonic degrees of freedoms. All these prospectives are inspiring upcoming experiments and research on the nuclear EOS [Li2014, Sorensen2023] based on nuclear dynamics.

In particular, nuclear dynamics investigates the evolution in time, and therefore as a function of density, of a fermionic complex from its formation until its possible fragmentation throughout successive transformations and combinations. The richest variety of evolution paths is achieved when largeamplitude fluctuations come into play. They result from violent perturbations like phase transitions or frustration phenomena, typically induced by external actions, like heavy-ion collisions. There exists a large history of modeling approaches for nuclear dynamics, originally mostly based on statistical approximations and later evolving to more and more microscopic approaches. In some circumstances, when the timescale of the process largely exceeds the relaxation time of the nuclear interaction, statistical approximations are fully justified and very efficient. This is typically the case of astrophysical processes, but statistical descriptions also fit relatively slow processes in heavy-ion collisions, like the decay by fission or particle emission of an excited nuclear complex. Quite generally, however, the initial stages of heavy-ion collisions evolve in so short timescales that chaotic behaviours dominate. In this case, explanatory theories for finite systems should not start from equilibrium assumptions. This is the case of transport approaches which we address in this thesis. The observables which could be studied experimentally vary as a function of the incident energy, especially around Fermi energy. Below Fermi energy, deep inelastic regimes are mainly associated with collective modes. At low energies, up to 15 AMeV, the probing wavelength is larger than the target nucleus which means that mean-field effects dominate over two-body nucleon-nucleon collisions. In green it is indicated that collisionless models such as the TDHF models describe very well this energy regime. At high energies, around 1 AGeV, the probing wavelength is of the size of the nucleonic interaction which means that collisional effects dominate the dynamical process, which is indicated in red as the participant-spectator regime. Lastly in blue, around Fermi energy, the probing wavelength is of the size of the target nucleus which shows that both mean field as well as two-body correlations have to be taken into account, which means that beyond-mean-field models have to be used to describe this regime accurately. See text.

Transport models for heavy-ion collisions, from collective modes to dissipative regimes

In this case, collective effects dominate and a corresponding modelling framework relies on mean-field features while two-body nucleon-nucleon collisions are highly suppressed by Pauli-blocking factors. Well suited to low energies are therefore ('collisionless') one-body theories, like time-dependent Hartree-Fock (TDHF) approaches [Ring1980, Goeke1983, Suraud1995, Maruhn2010, Simenel2012, Umar2017], or the semi-classical analogue Vlasov approach [Tang1981].

The opposite situation characterize largest incident energies, ranging from about 200 MeV per nucleon to relativistic energies, where the probing wavelength is of the order of the nuclear-interaction range. This situation distinguishes the intermediate-energy regime or, at relativistic energies, the participantspectator regime, where two-body nucleon-nucleon collisions gradually dominate over mean-field effects. Well suited to this energy regime are Molecular Dynamics (MD) approaches, either Anti-symmetrized Molecular Dynamics (AMD [Ono1992, Kanada1995]), Fermionic Molecular Dynamics (FMD [Feldmeier1990, Feldmeier1995]) or semi-classical Quantum Molecular Dynamics (QMD [Aichelin1986, Aichelin1991, Bohnet1991, Peilert1989, Peilert1991, Hartnack1989, Maruyama1998, Chikazumi2001, Papa2001]), where nucleon-nucleon correlations are explicitly treated while the one-body description is approximated. Due to a dominant thermodynamical contribution in the nuclear process which heats up the system, also statistical models are successfully applied to this energy regime.

Intermediate between the two above situations, the domain around Fermi energy is identified by a probing wavelength of the size of heavy ions. The Fermi-energy regime is the place of dissipative heavy-ion collisions where mechanical contributions (like deformations) are dominant and where the action of large-amplitude fluctuations results in the largest variety of mechanisms and exit channels, from fusion and deep inelastic to breakup into fragments and clusters. The dissipative character of the process requires to combine both mean-field and two-body effects. For this purpose, semi-classical transport approaches have been rather successful, both based on a Boltzmann one-body treatment or on Molecular Dynamics.

Based on the Boltzmann equation, one class of models is the Boltzmann-Uehling-Uhlenbeck (BUU) or Boltzmann-Nordheim-Vlasov (BNV) transport model [Cassing1990, Bonasera1994], which treats the nucleons as classical particles and follows their trajectories in phase space, taking into account their evolution within the mean-field potential and introducing dissipation through the additional contribution on nucleon-nucleon collisions, sampled via the test-particle method [Bertsch1988]. Alternatively, Molecular Dynamics models introduce two-body correlations by treating the nucleonic degrees of freedom explicitly. The AMD approach uses a variational approach to solve the many-body Schrödinger equation, treating the nucleons as correlated quantum particles rather than semi-classical particles.

The additional contribution from large-amplitude fluctuations which characterize the Fermi-energy domain requires to include also stochastic terms. In the Boltzmann sector, they could be treated in the spirit of the Brownian motion, in approaches based on the Boltzmann-Langevin (BL) equation [Ayik1990, Reinhard1992, Reinhard1992bis, Chomaz2004, Colonna1994a]. In approaches like the Boltzmann-Langevin-One-Body (BLOB) model, Beyond-mean-field extensions are introduced via a fluctuation contribution in phase space which is consistent with the dispersion relation of the unstable modes [Napolitani2013, Napolitani2017, Colonna2017]. This stochastic treatment in a one-body framework allows to describe fluctuations, bifurcations, the formation of fragments and nuclear clusters from the inhomogeneities emerging in the density landscape. Alternatively, in the MD sector, fluctuations are directly linked to its explicit many-body correlation treatment. In the more advanced approach of AMD, a quantum branching processes [Ono1999] is used as a random source in analogy to a Langevin term and has the effect to treat dissipation and fragment formation. In addition, three-body contributions are added in the nucleon-nucleon collisions term in order to describe the formation of nuclear cluster [Ono2016, Ono2019].

Limitations and common assumptions in current transport approaches

Despite the success of the above kinetic theories, in particular approaches like AMD and BLOB, there are still many challenges in simulating heavy-ion collisions in a broad range of incident energies. By construction, because of their intrinsic approximations, transport models based on the Boltzmann equation are semi-classical approaches that disregard quantum effects, while the AMD model provides a limited description of mean-field effects. In particular, all the above mentioned models present difficulties when applied to encompass different energy regimes, for instance when the aim is to describe low energies, Fermi energies and intermediate energies within the same unified framework. For example, at the threshold between Fermi and low energies, different types of instabilities compete leading to an intricate picture. We think that this limitation could be related to how non-locality is approximated both in MD and Boltzmann theories. While this approximation affects both theoretical frameworks, it comes about in different forms according to specific model implementations. .

MD models approximate non-locality by localizing the nucleonic degrees of freedom in order to keep nucleonic packets compact. This favors fragment and cluster formation via a more effective handling of two-body and three-body nucleon-nucleon correlations. As a drawback, localized nucleonic degrees of freedom affect significantly the collective behaviour of the system or, more generally, the description of the mean-field properties. This is evidently the case in semi-classical QMD models, but FMD or AMD suffer from this drawback as well [Feldmeier2000]. In fact, even though the anti-symmetrization, explicitly included by calculating the Slater determinant of the N-body system, leads to the direct treatment of Pauli-principle, it does not restore the mean-field phenomenology. As it will be justified in the following chapters, even though the localized nucleonic wave functions are anti-symmetrized, the available function space is still restricted to the span of localized basis functions. Thus, MD models are more suited to describe phenomena at intermediate to high energies, where short-range interactions dominate over mean-field contributions.

Mean-field equations, when dealing with large-amplitude dynamics, are usually translated into a semi-classical form in order to treat complex beyond-mean-field terms, like two-body collisions and fluctuating contributions. As a drawback, the semi-classical evolution tends to drive the initial Fermi distribution into a Boltzmann distribution at large times, with the consequence that stability is reduced when low energies are described. In addition, the way residual contributions are usually implemented impose some approximation of non-locality. Hereafter, we illustrate briefly how such approximations come about.

At low energy, in conditions of small-amplitude fluctuations, the mean-field paths which describe the dynamics of the system could be treated as correlated channels in approaches like the time-dependent generator coordinate method (TDGCM) [Reinhard1987, Goutte2005], where states are propagated coherently, or they can be described through a variational approach à la Balian-Vénéroni [Balian1981, Balian1992, Simenel2012]. Beyond this single-particle picture, more excited systems undergo largeamplitude dynamics, where dissipative effects from two-body nucleon-nucleon collisions should be additionally considered. In this conditions, a nuclear system may experience unstable regimes and explore several degrees of freedom leading to chaotic regimes and bifurcations; the formation of IMF is a signature of such behaviour. We can progress from a general picture of the hierarchy of correlations that are required to build up and extend a kinetic theory, such as the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy for density matrices [Bogoliubov1946, Born1946, Kirkwood1946, Balescu1976]. Kinetic equation, either of the BUU type, or based on dissipative extensions of the TDHF approach [Wong1978, Wong1979, Lacroix2004, Besse2020], can be associated to truncations in the BBGKY hierarchy up to two-body correlations [Cassing1992, Campa2014], while mean-field approaches result from only the leading one-body contribution. Large-amplitude regimes which may lead to bifurcations are related to higher order correlation, even beyond kinetic equations (i.e. beyond the second-order truncation of the BBGKY hierarchy). However, since treating higher-order correlations specifically is technically too intricate, it was proposed to adopt a coherent description in a statistical framework, leading to the formulation of stochastic TDHF (STDHF) [Reinhard1992, Suraud2014, Slama2015]. In order to achieve a semi-classical scheme, a Wigner-transform is required to describe the evolution of the phase-space density rather than the original quantum density operator. The resulting form is the Boltzmann Langevin (BL) Equation [Ayik1988], which is the semi-classical analogue of STDHF [Reinhard1992, Abe1996, Reinhard1992bis]. It includes the fluctuation 𝛿𝐼 [𝑓 ] of the average correlations. This procedure should ensure that Pauli blocking is fully satisfied in the collision but the geometry of the packets is not tracked in the mean field and it should therefore be defined at each time step on the basis of the efficiency of the scattering. In order to better agitate the phase space, it turns out that the 'semi-classical' wave-packets involved in the collision should be compact in phase space.

collision integral ⟨𝐼 coll [𝑓 ]⟩, 𝜕𝑓 1 𝜕𝑡 + ì 𝑝 1 𝑚 • ì ∇𝑓 1 -ì ∇𝑈 • ì ∇ 𝑝 1 𝑓 1 = ⟨𝐼 coll [𝑓 ]⟩ + 𝛿𝐼 [𝑓 ] . ( 1
This condition of compactness is however not consistent with the non-local character of a quantum mean-field.

1.5 Purpose: merging mean-field dynamics and nucleonic degrees of freedom Therefore, it is timely to develop a new theoretical model that can take into account both quantum effects, like the delocalized character of wave functions and their interference, as well as the resulting quantum fluctuations, in addition to an accurate sampling of the mean-field evolution. In this respect, progressing from the TDHF theory is a promising start to overcome these limitations. The TDHF approach solves the time-dependent Schrödinger equation for the many-body wave function in coordinate space. This model treats the nucleons as delocalized quantum particles and includes the mean-field potential. It can provide a clear description of both the quantum effects and mean-field effects. However, the TDHF theory still has some limitations, particularly because the extension to large-amplitude, beyond-mean-field effects is more complicated than in a semi-classical Boltzmann model. For this reason it is usually applied to nuclear structure at low energies.

To overcome these limitations, a new quantum transport model based on TDHF can be developed by defining delocalized wave functions and a suited decomposition into a dynamical basis. Progressing from such a construction, we set up a framework which allows to adopt some successful methods from BL approaches, like stochastic collision terms, and by MD approaches, like the explicit handling of nucleon-nucleon correlations, whilst keeping the quantum properties of the TDHF theory. Therefore, an attempt is addressed to combine and generalize features which characterize both the time-dependent Hartree-Fock theories and the anti-symmetrized molecular models.

The purpose of this thesis is to focus on developing the base of a theoretical framework which covers the essential quantum-mechanical features but also opens an avenue to possibly include large-amplitude fluctuations. We present a quantum-mechanical transport model based on the theory of TDHF which transitions to a generalized molecular dynamics model which keeps de-localized fermionic wave functions for a better description of the mean field. It is the quantum analogue to semi-classical models where test particles are used to describe the density distribution, however, contrary to semi-classical models, the test particles are not independent of each other but we keep track of the information of nucleonic degrees of freedom. That means, some test particles belong to the same nucleonic wave function. As a result, the mean-field evolution equations are modified such that quantum effects are directly included in the equation of motion, which allows us to explore low energies as well. Finally, we want to provide a numerical tool which improves on semi-classical models in the sense that we evolve time-dependent test particles in a quantum mean-field but additionally provide the information of de-localized nucleonic wave functions. The goal is to give a new foundation which handles explicitly the collective behaviour of low energy phenomena but is able to be extended by beyond-mean-field effects due to the knowledge of nucleonic degrees of freedom. However, the main focus in the following will be setting up such a theoretical approach starting from the underlying equation of motion and leaving the extension with collisional dynamics to be an important task in the future.

Outline

In Chapter 2 we will give a formal introduction and a detailed formulation of the quantum meanfield evolution equations of the dynamical bases which build up the de-localized wave packets. Since the focus of the thesis is to develop a rigorous quantum-mechanical many-body transport scheme, we present carefully each important step in solving the equation of motion. It starts at introducing the parametrization of the dynamical basis, where we also discuss important properties. Afterwards, we calculate the the essential gradients that appear in the evolution equation. Finally, it is also necessary to discuss how we initialize the nuclei properly to start an evolution. We use a self-consistency routine based on the harmonic oscillator and finish with a cooling procedure based on the imaginary-time evolution. We follow up in Chapter 3 with a discussion of the history of the numerical code that we developed. We discuss how it started as a model called DYWAN [Jouault1998] which already had the same idea that we follow to construct the model many years ago. However, as it turns out, that model went through a decoherence approximation which transformed it to a typical semi-classical model. Therefore, it is important to show in this chapter how the code works in detail. We start with the detailed explanation of the harmonic oscillator routine which includes the fitting process of the nuclear potential with an harmonic oscillator potential. After the routine converges we present the decomposition of the harmonic oscillator states into dynamical Gaussians and how three-dimensional states are constructed. Lastly, we explain for two cases the cooling procedure. We will discover how additional Gaussians add interference contributions to the evolution of other Gaussians, which shows the direct influence of quantum effects in the equation of motion, which are not present in semi-classical models. Finally, we show how two nuclei are initialized and how invariant transformations are defined in our parametrization. In Chapter 4 we show some tests which validates the foundation of the scheme to open up a new avenue for beyond-mean-field effects. We present the result of the cooling and the time-evolution of nuclei. It turns out that the cooling is necessary to avoid the initial structure which was given by the harmonic oscillator.

Chapter 2

Theoretical description of the model

|Ψ⟩ = A (|𝜑 1 ⟩ ⊗ . . . ⊗ |𝜑 𝐴 ⟩) . (2.3)
To account for the lack of additional Slaters as in eq. 2.2, the problem shifts to solving for the appropriate set of single-particle wave functions |𝜑 𝑖 ⟩ to describe the system. In TDHF models, the single-particle wave functions are determined by imposing orthonormality. Alternatively, in molecular dynamics [Feldmeier1995, Dote1997, Ono1998, Hasnaoui2008] a specific parametrization of the onebody wave function, usually of gaussian type, is introduced. Both approaches have their advantages and disadvantages. For example, in TDHF the orthonormality of the single-particle states lead to a simplification of the one-body density operator

ρ = 𝐴 Tr 2...𝐴 |Ψ⟩ ⟨Ψ| = 𝐴 ∑ 𝑖 𝑗 |𝜑 𝑖 ⟩ ⟨𝜑 𝑗 | 𝑂 -1 𝑖 𝑗 (2.4) = 𝐴 ∑ 𝑖 |𝜑 𝑖 ⟩ ⟨𝜑 𝑖 | (2.5)
where eq. 2.4 is the general solution of the density operator of a single-Slater determinant with 𝑂 𝑖 𝑗 being the overlap matrix of the one-particle states. In the special case of TDHF the overlap matrix reduces to 𝑂 𝑖 𝑗 = 𝛿 𝑖 𝑗 which is why the density operator simplifies to eq. 2.5. Furthermore, the single-particle wave functions are not restricted to a specific functional form resulting in a very good description of the mean-field properties. On the other hand, molecular dynamics models treat individual channels independently due to localization. That means, if a nucleon has been emitted, the residue nucleus is made up by 𝐴 -1 nucleons contrary to TDHF where only a part of a single-particle state will go out of the nucleus and the other part will remain in the residue nucleus [Ono2019]. In our model, we will follow the TDHF approach in order to keep a good description of the mean field but then we will switch to a parametrization of Gaussians in order to treat different channels accordingly.

Derivation of the mean-field equations of motion

In this section, we will present all the important steps that lead to the model that we are using and discuss in detail the impact on it. First of all, we start, as in any microscopic model, with the Lagrangian of the system and perform a variational principle [Balian1981] L = Tr (𝑖 h𝜕 𝑡 -ĥ) ρ (2.6)

= 𝐴 ∑ 𝑖 ⟨𝜑 𝑖 |𝑖 h𝜕 𝑡 -ĥ|𝜑 𝑖 ⟩ = 𝐴 ∑ 𝑖 ⟨𝜑 𝑖 |𝑖 h𝜕 𝑡 |𝜑 𝑖 ⟩ -E tot (2.7)
where ĥ is the single-particle Hamiltonian and E tot is the total energy of the system. As mentioned in the previous section, we follow the philosophy of TDHF approaches, where we look for sets of orthogonal nucleonic wave functions, leading us to eq. 2.7 using the simplified density operator in eq. 2.5. In order to find the single-particle wave functions, we choose to express them in a linear combination of parameterized functions 𝑔[𝑞 𝜈 ]

|𝜑 𝑖 ⟩ = 𝑁 𝑖 ∑ 𝑛 𝑔 𝑖 𝑛 ≡ 𝐼 ∑ 𝑛 |𝑔 𝑛 ⟩ (2.8)
where the number of basis functions 𝑁 𝑖 ≡ 𝐼 may vary for every nucleonic wave function. Furthermore, we introduce a convention that the sum over 𝐼 implies that the sum will go through the set of parameterized functions which belong to the specific nucleon 𝑖, which is why we omit the index in |𝑔 𝑛 ⟩ for the sake of clarity. However, we have to keep in mind that the drawback of simplifying the Lagrangian in eq. 2.7 is that the single-particle wave functions are additionally constrained by the orthogonality condition. Hence, those wave functions cannot access the full function space that is spanned by the basis functions but are restricted to a subspace of the manifold.

Finally, the equations of motion are given by the variational principle using the Euler-Lagrange equa- (2.11)

In eq. 2.10 we can get an insight on how the nucleonic wave functions evolve. In fact, the onebody Schrödinger equation has to be fulfilled along all directions of the tangent space that means all the parametric derivatives 𝜕 𝑞 𝜇 of every basis function 𝑔 𝑚 belonging to that nucleonic wave function.

Therefore, there are in total 𝐼 • 𝑃 different equations where 𝐼 is the number of gaussians and 𝑃 the number of parameters. If we rewrite the time-derivative of the single-particle wave function 𝜕 𝑡 𝜑 𝑖 as the sum over the time-derivative of the parameters 𝑞 𝑛 𝜈 ≡ 𝑞 𝜈 using the chain rule, we obtain eq. 2.11. Again, we will apply a convention that the greek letter 𝜈 implies that the parameter belongs to the specific gaussian 𝑛 which is why we drop the extra superscript. As a result, we find 𝐼 • 𝑃 unknown variables 𝑞 𝜈 for which we have the exact amount of equations to solve them. The only things that we need to know are the Hesse matrix elements 𝜕 𝑞 𝜈 𝑔 𝑛 𝜕 𝑞 𝜇 𝑔 𝑚 and the derivatives of the energy expectation values 𝜕 𝑞 𝜇 E tot . Thus, at this point it is necessary to choose the parameterization for the basis function in order to compute those matrix elements.

The dynamical basis function and their properties

First of all, we have to note that the evolution equations in eq. 2.11 are driven by the imaginary part of the Hesse matrix elements. In order to have a unique solution of the time-evolution of the parameters, that matrix has to be invertible. Since the scalar product is defined as

⟨𝑔 𝑛 |𝑔 𝑚 ⟩ = 𝑔 𝑛 (ì 𝑥)𝑔 𝑚 (ì 𝑥)d 3 𝑥 = ⟨𝑔 𝑚 |𝑔 𝑛 ⟩ (2.12) =⇒ ℑ 𝜕 𝜈 𝑔 𝑛 𝜕 𝜇 𝑔 𝑚 = ℑ 𝜕 𝜇 𝑔 𝑚 𝜕 𝜈 𝑔 𝑛 = -ℑ 𝜕 𝜇 𝑔 𝑚 𝜕 𝜈 𝑔 𝑛 (2.13)
the evolution matrix ℑ 𝜕 𝜈 𝑔 𝑛 𝜕 𝜇 𝑔 𝑚 will be anti-symmetric. As a consequence, that matrix has to be of even dimension in order to be invertible. We recall that the dimension is given by 𝐼 • 𝑃 and since the number of gaussians 𝐼 can take any natural number we follow that the number of parameters 𝑃 has to be even. This is of no surprise, because it is known that in analytical mechanics you need pairs of conjugate variables such as position and momentum to describe the system. With that in mind we will choose the following parameterization of the dynamical basis function:

𝑔 𝑚 (𝑥) = 1 2𝜋 𝜒 𝑚 1/4 exp -𝜉 𝑚 (𝑥 -𝑥 𝑚 ) 2 2 + 𝑖𝑘 𝑚 (𝑥 -𝑥 𝑚 )) (2.14) with 𝜉 𝑚 ≡ 1 2𝜒 𝑚 -2𝑖𝛾 𝑚 , 𝛾 𝑚 ≡ 𝜎 𝑚 2𝜒 𝑚 𝑔 𝑚 (ì 𝑥) = 𝑐 𝑚 𝑒 𝑖𝜃 𝑚 𝑔 𝑥 𝑚 (𝑥)𝑔 𝑦 𝑚 (𝑦)𝑔 𝑧 𝑚 (𝑧).
(2.15)

As we can see, we use a standard one-dimensional gaussian-type parameterization in eq. 2.14 with the time-dependent parameters {𝑥 𝑚 , 𝑘 𝑚 , 𝜒 𝑚 , 𝛾 𝑚 }. Additionally we multiply a complex weight 𝑐 𝑚 𝑒 𝑖𝜃 𝑚 in eq. 2.15 after we combine all gaussians of each three dimensions together. In total there are fourteen parameters to describe one gaussian basis function. Concerning the physical meaning behind the parameters, let us look at some properties of this parametrization:

|𝑔 𝑚 | 2 (𝑥) = 1 2𝜋 𝜒 𝑚 1/2 exp - (𝑥 -𝑥 𝑚 ) 2 2𝜒 𝑚 (2.16) ĝ𝑚 (𝑘) = 1 √ 2𝜋 𝑔 𝑚 (𝑥)𝑒 -𝑖𝑘𝑥 d𝑥 = 1 2𝜋 𝜒 𝑚 𝜉 2 𝑚 1/4 exp - (𝑘 -𝑘 𝑚 ) 2 2𝜉 𝑚 -𝑖𝑘𝑥 𝑚 (2.17) | ĝ𝑚 | 2 (𝑘) = 1 2𝜋𝜙 𝑚 1/2 exp - (𝑘 -𝑘 𝑚 ) 2 2𝜙 𝑚 , 𝜙 𝑚 ≡ 𝜒 𝑚 |𝜉 𝑚 | 2 .
(2.18)

In the first eq. 2.16 we can see that 𝜒 𝑚 represents the variance of the spatial probability distribution and 𝑥 𝑚 its position, also shown in fig. 2.1. Next, we show the Fourier-transform and its momentum probability distribution in eq. 2.17 and eq. 2.18. We notice that 𝑘 𝑚 corresponds to the angular wave number and 𝜙 𝑚 ≡ 𝜒 𝑚 |𝜉 𝑚 | 2 to the variance of the momentum probability distribution. Finally, if we check the product of the variances

𝜒 𝑚 𝜙 𝑚 = 𝜒 2 𝑚 |𝜉 𝑚 | 2 = 1 4 + 4𝜒 2 𝑚 𝛾 2 𝑚 = 1 4 + 𝜎 2 𝑚 ≥ 1 4 (2.19)
it shows that for 𝛾 𝑚 = 0, or equivalently 𝜎 𝑚 = 0, the gaussian has minimum uncertainty and for any other value 𝛾 𝑚 ̸ = 0 the uncertainty can become arbitrarily large. That means 𝛾 𝑚 is a parameter which indicates how much more than (h/2𝜋) 3 the gaussian occupies a phase-space volume. Furthermore, we can rewrite the parameterization of a one-dimensional gaussian basis function into

𝑔 𝑚 (𝑥) = 1 2𝜋 𝜒 𝑚 1/4 exp          -𝜉 𝑚 𝑥 -(𝑥 𝑚 + 𝑖 𝑘 𝑚 𝜉 𝑚 ) 2 2 - 𝑘 2 𝑚 2𝜉 𝑚          . (2.20)
where we rearranged the square in the exponent of the gaussian in eq. 2.14. In this form we can easily recognize that the center of the gaussian is located in the complex plane at (𝑥 𝑚 + 𝑖 𝑘 𝑚 𝜉 𝑚 ). Note that the positions of the gaussians in the complex plane does not directly correlate to the position (𝑥 𝑚 , 𝑘 𝑚 ) in phase-space since 𝜉 𝑚 is complex, thus there is an additional contribution to the real part in (𝑥 𝑚 +𝑖 𝑘 𝑚 𝜉 𝑚 ). Nevertheless, that fact helps us to visualize the scalar product between two gaussians. As shown in eq. 2.12 it is defined as the integral of the product function 𝑔 𝑛 (ì 𝑥)𝑔 𝑚 (ì 𝑥) along the real axes. Since we take the complex conjugate of the first gaussian 𝑔 𝑛 the center of that gaussian is affected accordingly.

To be specific, the center of 𝑔 𝑛 is located in the complex plane at (𝑥 𝑛 -𝑖 𝑘 𝑛

𝜉 𝑛

). Moreover, the product of two gaussians can be re-expressed into a new gaussian 𝑔 𝑛𝑚 whose new center is positioned in-between the centers of the two defining gaussians at (𝑥 𝑛 -𝑖 𝑘 𝑛 𝜉 𝑛

) and (𝑥 𝑚 + 𝑖 𝑘 𝑚 𝜉 𝑚 ). In detail we find

𝑔 𝑛𝑚 (𝑥) = 𝑔 𝑛 (𝑥)𝑔 𝑚 (𝑥) = 1 4𝜋 2 𝜒 𝑛 𝜒 𝑚 1/4 exp - (𝑥 -𝑋 𝑛𝑚 ) 2 2𝜒 𝑛𝑚 + 𝑋 2 𝑛𝑚 2𝜒 𝑛𝑚 + 𝐶 𝑛𝑚 (2.21)
with 𝜒 𝑛𝑚 = 1

𝜉 𝑛 + 𝜉 𝑚 (2.22) 𝑋 𝑛𝑚 = 1 𝜉 𝑛 + 𝜉 𝑚 𝜉 𝑛 (𝑥 𝑛 -𝑖 𝑘 𝑛 𝜉 𝑛 ) + 𝜉 𝑚 (𝑥 𝑚 + 𝑖 𝑘 𝑚 𝜉 𝑚 ) = 1 𝜉 𝑛 + 𝜉 𝑚 𝜉 𝑛 𝑥 𝑛 + 𝜉 𝑚 𝑥 𝑚 + 𝑖 (𝑘 𝑚 -𝑘 𝑛 ) (2.23) 𝐶 𝑛𝑚 = -𝜉 𝑛 𝑥 2 𝑛 2 -𝜉 𝑚 𝑥 2 𝑚 2 + 𝑖 (𝑘 𝑛 𝑥 𝑛 -𝑘 𝑚 𝑥 𝑚 ) (2.24)
where 𝜒 𝑛𝑚 is the variance and 𝑋 𝑛𝑚 is the position of the product function. We recognize in eq. 2.23 that the position is determined by the positions of the two gaussians 𝑔 𝑛 and 𝑔 𝑚 weighted by the inverse variances 𝜉 𝑛 and 𝜉 𝑚 . A schematic visualization of this process is shown in fig. 2.2. Each of these double indexed variables are also self-adjoint in those indexes

𝜒 𝑛𝑚 = 𝜒 𝑚𝑛 , 𝑋 𝑛𝑚 = 𝑋 𝑚𝑛 , 𝐶 𝑛𝑚 = 𝐶 𝑚𝑛 .
(2.25)

Finally, from eq. 2.21 we find an expression for the scalar product

⟨𝑔 𝑛 |𝑔 𝑚 ⟩ 𝑖 = 𝜒 2 𝑛𝑚,𝑖 𝜒 𝑛,𝑖 𝜒 𝑚,𝑖 1/4 exp 𝑋 2 𝑛𝑚,𝑖 2𝜒 𝑛𝑚,𝑖 + 𝐶 𝑛𝑚,𝑖 (2.26) ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ = 𝑐 𝑛 𝑐 𝑚 𝑒 𝑖 (𝜃 𝑚 -𝜃 𝑛 ) ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ 𝑥 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ 𝑦 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ 𝑧 (2.27)
where ⟨•|•⟩ 𝑖 denotes the one-dimensional scalar product in the dimension 𝑖 = {𝑥, 𝑦, 𝑧}. The threedimensional scalar product is the product of all the one-dimensional scalar products. The exponential can be further broken down with the help of the following relationship

𝑋 𝑛𝑚 (𝑥 𝑛 + 𝑥 𝑚 ) = 1 𝜉 𝑛 + 𝜉 𝑚 𝜉 𝑛 𝑥 2 𝑛 + 𝜉 𝑚 𝑥 2 𝑚 + (𝜉 𝑛 + 𝜉 𝑚 )𝑥 𝑛 𝑥 𝑚 + 𝑖 (𝑘 𝑚 -𝑘 𝑛 )(𝑥 𝑛 + 𝑥 𝑚 ) (2.28) ⇒ 2𝜒 𝑛𝑚 𝐶 𝑛𝑚 = 1 𝜉 𝑛 + 𝜉 𝑚 -𝜉 𝑛 𝑥 2 𝑛 -𝜉 𝑚 𝑥 2 𝑚 + 2𝑖 (𝑘 𝑛 𝑥 𝑛 -𝑘 𝑚 𝑥 𝑚 ) = -𝑋 𝑛𝑚 (𝑥 𝑛 + 𝑥 𝑚 ) + 𝑥 𝑛 𝑥 𝑚 + 𝑖 (𝑘 𝑚 + 𝑘 𝑛 )(𝑥 𝑛 -𝑥 𝑚 ) 𝜉 𝑛 + 𝜉 𝑚 (2.29) Figure 2
.2: Schematic picture of the derivation of the position of the product function 𝑔 𝑛𝑚 . The red cross represents the central position of the gaussian 𝑔 𝑚 in the complex plane and the black cross the central position of the gaussian 𝑔 𝑛 . The blue cross is the position of the complex conjugate 𝑔 𝑛 . The product function 𝑔 𝑛𝑚 , in green, is the average position between red and blue, which are weighted by the complex inverse variances 𝜉 𝑚 and 𝜉 𝑛 which is why the green position does not necessarily have to be located directly on the line in-between red and blue. so that we finally obtain for the exponent of the scalar product

𝑋 2 𝑛𝑚 2𝜒 𝑛𝑚 + 𝐶 𝑛𝑚 = (𝑋 𝑛𝑚 -𝑥 𝑛 )(𝑋 𝑛𝑚 -𝑥 𝑚 ) 2𝜒 𝑛𝑚 + 𝑖 (𝑘 𝑚 + 𝑘 𝑛 )(𝑥 𝑛 -𝑥 𝑚 ) 2 .
(2.30)

It is instructive to take a closer look on that expression. Especially because the product (𝑋 𝑛𝑚 -𝑥 𝑛 )(𝑋 𝑛𝑚 -𝑥 𝑚 ) will appear a lot in the next subsection where we calculate the expectation values of derivatives. Therefore, we first express each factor individually in terms of the differences (𝑥 𝑚 -𝑥 𝑛 ) and

(𝑘 𝑚 -𝑘 𝑛 ) (𝑋 𝑛𝑚 -𝑥 𝑛 ) = 1 𝜉 𝑛 + 𝜉 𝑚 𝜉 𝑚 (𝑥 𝑚 -𝑥 𝑛 ) + 𝑖 (𝑘 𝑚 -𝑘 𝑛 ) (2.31) (𝑋 𝑛𝑚 -𝑥 𝑚 ) = 1 𝜉 𝑛 + 𝜉 𝑚 -𝜉 𝑛 (𝑥 𝑚 -𝑥 𝑛 ) + 𝑖 (𝑘 𝑚 -𝑘 𝑛 ) (2.32)
and we find that they have a very symmetric expression. The product of those two equation results in

(𝑋 𝑛𝑚 -𝑥 𝑛 )(𝑋 𝑛𝑚 -𝑥 𝑚 ) = - 1 (𝜉 𝑛 + 𝜉 𝑚 ) 2 𝜉 𝑛 𝜉 𝑚 (𝑥 𝑚 -𝑥 𝑛 ) 2 + (𝑘 𝑚 -𝑘 𝑛 ) 2 -𝑖 (𝜉 𝑚 -𝜉 𝑛 )(𝑥 𝑚 -𝑥 𝑛 )(𝑘 𝑚 -𝑘 𝑛 )
which resembles a distance in a metric space. This is not surprising since the origin of this term is the scalar product between two gaussians. Therefore the magnitude of the scalar product should be proportional to the distance between those gaussians. As it is known, the magnitude of a complex number written as 𝑒 𝑧 is given by the real part of the exponent z. Thus, the real part of eq. 2.30 only comes from the first addend (𝑋 𝑛𝑚 -𝑥 𝑛 )(𝑋 𝑛𝑚 -𝑥 𝑚 ) since the second addend is purely imaginary.

Jacobian and Hesse matrix elements

After we introduced the parameterization of the basis function of our choice in the previous subsection, we are now able to compute the evolution matrix elements in eq. 2.11. Therefore, we need the derivatives with respect to every parameter in

P 𝑚 ≡ { ì 𝑘 𝑚 , ì 𝛾 𝑚 , ì 𝑥 𝑚 , ì 𝜒 𝑚 , 𝑐 𝑚 , 𝜃 𝑚 } ì ∇ 𝑘 𝑚 ,𝑖 𝑔 𝑚 = 𝑖 (ì 𝑥 𝑖 -ì 𝑥 𝑚,𝑖 )𝑔 𝑚 (2.33) ì ∇ 𝛾 𝑚 ,𝑖 𝑔 𝑚 = 𝑖 (ì 𝑥 𝑖 -ì 𝑥 𝑚,𝑖 ) 2 𝑔 𝑚 (2.34) ì ∇ 𝑥 𝑚 ,𝑖 𝑔 𝑚 = -𝑖 ( ì 𝑘 𝑚,𝑖 + ì 𝜉 𝑚,𝑖 ì ∇ 𝑘 𝑚 ,𝑖 )𝑔 𝑚 (2.35) ì ∇ 𝜒 𝑚 ,𝑖 𝑔 𝑚 = - 1 4 ì 𝜒 𝑚,𝑖 (1 + 2𝑖 ì 𝜉 𝑚,𝑖 ì ∇ 𝛾 𝑚 ,𝑖 )𝑔 𝑚 (2.36) 𝜕 𝑐 𝑚 𝑔 𝑚 = 1 𝑐 𝑚 𝑔 𝑚 , 𝜕 𝜃 𝑚 𝑔 𝑚 = 𝑖𝑔 𝑚 (2.37)
where 𝑖 = 1, 2, 3 is the index for the dimension. The Hesse matrix ì ∇ 𝑞,𝑖 𝑔 𝑛 ì ∇ 𝑝,𝑗 𝑔 𝑚 with 𝑞 ∈ P 𝑛 and 𝑝 ∈ P 𝑚 is defined as the scalar product between derivatives of two gaussians. There are two cases to consider, first when the dimensions are different, that is 𝑖 ̸ = 𝑗 and second when the derivatives are in the same dimension 𝑖 = 𝑗. In the first case with 𝑖 ̸ = 𝑗 the scalar product of the Hesse matrix factorizes

into ì ∇ 𝑞,𝑖 𝑔 𝑛 ì ∇ 𝑝,𝑗 𝑔 𝑚 = ì ∇ 𝑞,𝑖 𝑔 𝑛 𝑔 𝑚 𝑖 𝑔 𝑛 ì ∇ 𝑝,𝑗 𝑔 𝑚 𝑗 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ 𝑘 (2.38) = ì ∇ 𝑞,𝑖 𝑔 𝑛 𝑔 𝑚 𝑔 𝑛 ì ∇ 𝑝,𝑗 𝑔 𝑚 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ (2.39)
where 𝑘 ̸ = 𝑖, 𝑗 is the third dimension. As a result we have to calculate the expectation values of the Jacobians 𝑔 𝑛 ì ∇ 𝑝,𝑖 𝑔 𝑚 :

𝑔 𝑛 ì ∇ 𝑘 𝑚 ,𝑖 𝑔 𝑚 = ⟨𝑔 𝑛 |𝑖 (ì 𝑥 𝑖 -ì 𝑥 𝑚,𝑖 )|𝑔 𝑚 ⟩ = ⟨𝑔 𝑛 |𝑖 (ì 𝑥 𝑖 -ì 𝑋 𝑛𝑚,𝑖 + ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑚,𝑖 )|𝑔 𝑚 ⟩ = 𝑖 ( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑚,𝑖 ) ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ (2.40) 𝑔 𝑛 ì ∇ 𝛾 𝑚 ,𝑖 𝑔 𝑚 = ⟨𝑔 𝑛 |𝑖 (ì 𝑥 𝑖 -ì 𝑥 𝑚,𝑖 ) 2 |𝑔 𝑚 ⟩ = 𝑖 ì 𝜒 𝑛𝑚,𝑖 + ( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑚,𝑖 ) 2 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ (2.41) 𝑔 𝑛 ì ∇ 𝑥 𝑚 ,𝑖 𝑔 𝑚 = -𝑖 ⟨𝑔 𝑛 |( ì 𝑘 𝑚,𝑖 + ì 𝜉 𝑚,𝑖 ì ∇ 𝑘 𝑚 ,𝑖 )|𝑔 𝑚 ⟩ (2.42) 𝑔 𝑛 ì ∇ 𝜒 𝑚 ,𝑖 𝑔 𝑚 = - 1 4 ì 𝜒 𝑚,𝑖 ⟨𝑔 𝑛 |(1 + 2𝑖 ì 𝜉 𝑚,𝑖 ì ∇ 𝛾 𝑚 ,𝑖 )|𝑔 𝑚 ⟩ (2.43) 𝑔 𝑛 𝜕 𝑐 𝑚 𝑔 𝑚 = 1 𝑐 𝑚 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ , 𝑔 𝑛 𝜕 𝜃 𝑚 𝑔 𝑚 = 𝑖 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ .
(2.44)

In the second case 𝑖 = 𝑗 the derivatives act in the same dimension, thus there are additional terms related to the interplay between the derivatives:

ì ∇ 𝑘 𝑛 ,𝑖 𝑔 𝑛 ì ∇ 𝑘 𝑚 ,𝑖 𝑔 𝑚 = ⟨𝑔 𝑛 |(ì 𝑥 𝑖 -ì 𝑥 𝑛,𝑖 )(ì 𝑥 𝑖 -ì 𝑥 𝑚,𝑖 )|𝑔 𝑚 ⟩ = ì 𝜒 𝑛𝑚,𝑖 + ( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑚,𝑖 )( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑛,𝑖 ) ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ = ì 𝜒 𝑛𝑚,𝑖 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ + ì ∇ 𝑘 𝑛 ,𝑖 𝑔 𝑛 𝑔 𝑚 𝑔 𝑛 ì ∇ 𝑘 𝑚 ,𝑖 𝑔 𝑚 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ (2.45) ì ∇ 𝛾 𝑛 ,𝑖 𝑔 𝑛 ì ∇ 𝛾 𝑚 ,𝑖 𝑔 𝑚 = ⟨𝑔 𝑛 |(ì 𝑥 𝑖 -ì 𝑥 𝑛,𝑖 ) 2 (ì 𝑥 𝑖 -ì 𝑥 𝑚,𝑖 ) 2 |𝑔 𝑚 ⟩ = 2 ì 𝜒 𝑛𝑚,𝑖 ì 𝜒 𝑛𝑚,𝑖 + 4 ì 𝜒 𝑛𝑚,𝑖 ( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑚,𝑖 )( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑛,𝑖 ) + ì 𝜒 𝑛𝑚,𝑖 + ( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑛,𝑖 ) 2 ì 𝜒 𝑛𝑚,𝑖 + ( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑚,𝑖 ) 2 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ = 2 ì 𝜒 𝑛𝑚,𝑖 ì 𝜒 𝑛𝑚,𝑖 + 4 ì 𝜒 𝑛𝑚,𝑖 ( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑚,𝑖 )( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑛,𝑖 ) ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ + ì ∇ 𝛾 𝑛 ,𝑖 𝑔 𝑛 𝑔 𝑚 𝑔 𝑛 ì ∇ 𝛾 𝑚 ,𝑖 𝑔 𝑚 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ (2.46) ì ∇ 𝛾 𝑛 ,𝑖 𝑔 𝑛 ì ∇ 𝑘 𝑚 ,𝑖 𝑔 𝑚 = ⟨𝑔 𝑛 |(ì 𝑥 𝑖 -ì 𝑥 𝑛,𝑖 ) 2 (ì 𝑥 𝑖 -ì 𝑥 𝑚,𝑖 )|𝑔 𝑚 ⟩ = 2 ì 𝜒 𝑛𝑚,𝑖 ( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑛,𝑖 ) + ì 𝜒 𝑛𝑚,𝑖 + ( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑛,𝑖 ) 2 ( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑚,𝑖 ) ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ = 2 ì 𝜒 𝑛𝑚,𝑖 ( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑛,𝑖 ) ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ + ì ∇ 𝛾 𝑛 ,𝑖 𝑔 𝑛 𝑔 𝑚 𝑔 𝑛 ì ∇ 𝑘 𝑚 ,𝑖 𝑔 𝑚 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ (2.47) ì ∇ 𝑥 𝑛 ,𝑖 𝑔 𝑛 ì ∇ 𝑝,𝑖 𝑔 𝑚 = 𝑖 ( ì 𝑘 𝑛,𝑖 + ì 𝜉 𝑛,𝑖 ì ∇ 𝑘 𝑛 ,𝑖 )𝑔 𝑛 ì ∇ 𝑝,𝑖 𝑔 𝑚 (2.48) ì ∇ 𝜒 𝑛 ,𝑖 𝑔 𝑛 ì ∇ 𝑝,𝑖 𝑔 𝑚 = - 1 4 ì 𝜒 𝑛,𝑖 (1 + 2𝑖 ì 𝜉 𝑛,𝑖 ì ∇ 𝛾 𝑛 ,𝑖 )𝑔 𝑛 ì ∇ 𝑝,𝑖 𝑔 𝑚 (2.49) 𝜕 𝑐 𝑛 𝑔 𝑛 ì ∇ 𝑝,𝑖 𝑔 𝑚 = 1 𝑐 𝑛 𝑔 𝑛 ì ∇ 𝑝,𝑖 𝑔 𝑚 , 𝜕 𝜃 𝑛 𝑔 𝑛 ì ∇ 𝑝,𝑖 𝑔 𝑚 = -𝑖 𝑔 𝑛 ì ∇ 𝑝,𝑖 𝑔 𝑚 .
(2.50)

The last three lines from eq. 2.48 to eq. 2.50 are formulas to compute the Hesse matrix element with an arbitrary parameter 𝑝 ∈ P 𝑚 . That is possible due to the fact that the derivative with respect to { ì 𝑥 𝑚 , ì 𝜒 𝑚 , 𝑐 𝑚 , 𝜃 𝑚 } can be expressed in terms of a constant plus derivatives with respect to { ì 𝑘 𝑚 , ì 𝛾 𝑚 }. Thus, essentially only the first three lines from eq. 2.45 to eq. 2.47 are necessary to build up the entirety of the Hesse matrix. In addition, the first three equations have been written in a way where the independent part has been separated to clearly distinguish the terms that are different to the non-diagonal case in eq. 2.39. Thus, the Hesse matrix elements can also be read as following

ì ∇ 𝑞,𝑖 𝑔 𝑛 ì ∇ 𝑝,𝑗 𝑔 𝑚 = ì ∇ 𝑞,𝑖 𝑔 𝑛 𝑔 𝑚 𝑔 𝑛 ì ∇ 𝑝,𝑗 𝑔 𝑚 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ + 𝛿 𝑖 𝑗 𝑀 𝑞 𝑖 ,𝑝 𝑖 (2.51)
where the first term is the independent part which factorizes into the product of the corresponding Jacobians and the second term is an additional mixing term which appears in case when the dimensions 𝑖 and 𝑗 coincide. For completeness, we summarize the mixing terms of the Hesse matrix elements

𝑀 𝑘 𝑛,𝑖 ,𝑘 𝑚,𝑖 = ì 𝜒 𝑛𝑚,𝑖 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ (2.52) 𝑀 𝛾 𝑛,𝑖 ,𝛾 𝑚,𝑖 = 2 ì 𝜒 𝑛𝑚,𝑖 ì 𝜒 𝑛𝑚,𝑖 + 4 ì 𝜒 𝑛𝑚,𝑖 ( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑚,𝑖 )( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑛,𝑖 ) ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ (2.53) 𝑀 𝛾 𝑛,𝑖 ,𝑘 𝑚,𝑖 = 2 ì 𝜒 𝑛𝑚,𝑖 ( ì 𝑋 𝑛𝑚,𝑖 -ì 𝑥 𝑛,𝑖 ) ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ (2.54) 𝑀 𝑥 𝑛,𝑖 ,𝑝 𝑖 = 𝑖 𝜉 𝑛,𝑖 𝑀 𝑘 𝑛,𝑖 ,𝑝 𝑖 (2.55) 𝑀 𝜒 𝑛,𝑖 ,𝑝 𝑖 = 𝑖 𝜉 𝑛,𝑖 2 ì 𝜒 𝑛,𝑖 𝑀 𝛾 𝑛,𝑖 ,𝑝 𝑖 (2.56) 𝑀 𝑐 𝑛 ,𝑝 𝑖 = 0 (2.57) 𝑀 𝜃 𝑛 ,𝑝 𝑖 = 0 (2.58)
where 𝑝 ∈ P 𝑚 . All in all, this concludes the derivation of the evolution matrix in eq. 2.11. However, to complete the full equation we have to compute the right-hand side as well. Those terms are related to derivatives of the expectation value of the hamiltonian. Hence, we have to introduce the nuclear potential that we will be using in our model.

Nuclear interaction and the energy expectation values

In mean-field approaches the nuclear interaction is often an effective, phenomenological interaction interactions, since they have shown a great success in mean-field calculations. The main difference between these two types of interactions is that the Gogny interaction has finite range whereas the Skyrme interaction is a contact force with zero range. In our framework, we choose a Skyrme-type interaction for the reason that the point-like interaction allows us to have a simple analytical expression of the energy functional in terms of the nucleon densities. In particular, the nuclear energy-density that we are using for our model is given by

H Skyrme (𝜌, 𝜁 ) = H 0 + H 3 + H surf (2.59) H 0 = 3 8 𝑡 0 𝜌 2 - 1 2 1 2 𝑡 0 𝑥 0 + 1 2 𝜁 2
(2.60)

H 3 = 1 16 𝑡 3 𝜌 𝜎+2 - 1 2 1 12 𝑡 3 𝑥 3 + 1 2 𝜌 𝜎 𝜁 2 (2.61) H surf = 1 16 9 4 𝑡 1 -𝑡 2 𝑥 2 + 5 4 ì ∇𝜌 2 - 1 32 3𝑡 1 𝑥 1 + 1 2 + 𝑡 2 𝑥 2 + 1 2 ì ∇𝜁 2 (2.62)
where 𝜌 = 𝜌 n +𝜌 p is the isoscalar density, 𝜁 = 𝜌 n -𝜌 p is the isovector density and {𝑡 0 , 𝑡 1 , 𝑡 2 , 𝑡 3 , 𝑥 0 , 𝑥 1 , 𝑥 2 , 𝑥 3 , 𝜎 } the parameters of the Skyrme interaction. The neutron and proton densities 𝜌 n,p are defined by the configuration space representation of the density operator which only include the nucleonic wave functions of the relevant isospin

𝜌 n,p = 𝑁 n,p ∑ 𝑖 ⟨ì 𝑥 |𝜑 𝑖 ⟩ ⟨𝜑 𝑖 | ì 𝑥⟩ = 𝑁 n,p ∑ 𝑖 𝜑 𝑖 (ì 𝑥)𝜑 𝑖 (ì 𝑥) . (2.63)
In this equation, it is implied that the sum runs through the subset of wave functions of the chosen isospin. Further, H 0 is related to the local density contribution, H 3 gives the many-body contribution and H surf involves the gradient of the densities and thus is associated to the surface term. We note that we are using an approximation of the Skyrme parameterization where we neglect for now the contributions coming from the effective mass and the spin-orbit effects. It is known that the spin-orbit effects are crucial to obtain the magic numbers beyond 𝑁 , 𝑍 > 20 [Ring1980]. However, the simplified potential can still reproduce approximately the main bulk nuclear ground state properties such as energy, radii and surface. Therefore, we decided for numerical reasons, to choose an approximated effective nuclear force for the dynamical evolution of the system.

The potential energy related to the nuclear interaction then reads

E Skyrme = ∭ H Skyrme d 3 𝑥 (2.64) 𝑈 n,p Skyrme ≡ 𝛿H Skyrme 𝛿𝜌 n,p (2.65)
where we introduce at this moment the effective Skyrme potential 𝑈 n,p Skyrme for neutrons and protons, respectively. Since the energy density is only a polynomial in 𝜌 n,p , it is straight-forward to perform the functional derivative

𝑈 n,p Skyrme = 3 4 𝑡 0 𝜌 -𝜏 n,p 1 2 𝑡 0 𝑥 0 + 1 2 𝜁 + 𝜎 + 2 16 𝑡 3 𝜌 𝜎+1 - 𝜎 2 1 12 𝑡 3 𝑥 3 + 1 2 𝜌 𝜎 -1 𝜁 2 -𝜏 n,p 1 12 𝑡 3 𝑥 3 + 1 2 𝜌 𝜎 -1 𝜁 - 1 8 9 4 𝑡 1 -𝑡 2 𝑥 2 + 5 4 ∆𝜌 + 𝜏 n,p 1 16 3𝑡 1 𝑥 1 + 1 2 + 𝑡 2 𝑥 2 + 1 2 ∆𝜁 (2.66) 𝜏 n = 𝜕𝜁 𝜕 n = 1 , 𝜏 p = 𝜕𝜁 𝜕 p = -1 (2.67)
where the effective neutron potential differs from the effective proton potential by a sign 𝜏 n,p which appears in the terms where the derivative acted on the isovector densities. In addition to the nuclear energy, we have to add the kinetic as well as the Coulomb energy in order to have the full expression of the total energy. If we want to implement the Coulomb interaction into our mean-field theory we have to remember that the Coulomb force is a long-range interaction, hence we have to include an exchange term in addition to the direct term. The direct term and exchange term are given in the Slater approximation [Slater1951] by

E Coul,direct = 1 2 ∭ 𝜌 p (ì 𝑥)𝑉 Coul (ì 𝑥) d 3 𝑥 (2.68) E Coul,exchange = - 𝑒 2 4𝜋𝜖 0 ∭ 3 4 𝜌 p 3 𝜋 𝜌 p 1/3 d 3 𝑥 (2.69) 𝑉 Coul (ì 𝑥) = 𝑒 2 4𝜋𝜖 0 ∭ 𝜌 p ( ì 𝑥 ′ ) | ì 𝑥 -ì 𝑥 ′ | d 3 𝑥 ′ (2.70)
where 𝜌 𝑝 is the proton density. To simplify our calculations we construct an effective Coulomb energy density from which we can define an effective Coulomb potential

H Coul = 1 2 𝜌 p (ì 𝑥)𝑉 Coul (ì 𝑥) - 𝑒 2 4𝜋𝜖 0 3 4 𝜌 p 3 𝜋 𝜌 p 1/3
(2.71)

𝑈 Coul (ì 𝑥) ≡ 𝛿H Coul 𝛿𝜌 p = 𝑉 Coul (ì 𝑥) - 𝑒 2 4𝜋𝜖 0 3 𝜋 𝜌 p (ì 𝑥) 1/3
(2.72)

E Coul = ∭ H Coul d 3 𝑥 . (2.73)
Eventually, the kinetic energy consists of the kinetic contribution of each nucleonic wave function.

Those contributions are calculated by the expectation value of the kinetic energy operator

⟨𝜑 𝑖 | T|𝜑 𝑖 ⟩ = ⟨𝜑 𝑖 | h2 ì 𝒌 2 2𝑚 |𝜑 𝑖 ⟩ (2.74) = - h 2𝑚 ∑ 𝑛𝑚 ⟨𝑔 𝑛 |( ì ∇) 2 |𝑔 𝑚 ⟩ = - h 2𝑚 ∑ 𝑛𝑚 ⟨𝑔 𝑛 |𝜕 2 𝑥 + 𝜕 2 𝑦 + 𝜕 2 𝑧 |𝑔 𝑚 ⟩ (2.75) = h 2𝑚 ∑ 𝑛𝑚 ⟨𝜕 𝑥 𝑔 𝑛 |𝜕 𝑥 𝑔 𝑚 ⟩ + 𝜕 𝑦 𝑔 𝑛 𝜕 𝑦 𝑔 𝑚 + ⟨𝜕 𝑧 𝑔 𝑛 |𝜕 𝑧 𝑔 𝑚 ⟩ (2.76) = h 2𝑚 ∑ 𝑛𝑚 𝜕 𝑥 𝑛 𝑔 𝑛 𝜕 𝑥 𝑚 𝑔 𝑚 + 𝜕 𝑦 𝑛 𝑔 𝑛 𝜕 𝑦 𝑚 𝑔 𝑚 + 𝜕 𝑧 𝑛 𝑔 𝑛 𝜕 𝑧 𝑚 𝑔 𝑚 (2.77)
where we applied the definition of the momentum operator ì 𝒌 = -𝑖 ì ∇ in eq. 2.74, then we performed an integration by parts to move one derivative to the left side of the braket and lastly we used the fact that 𝜕 𝑥 𝑔 𝑚 = -𝜕 𝑥 𝑚 𝑔 𝑚 in eq. 2.76. The final result is that we can express the kinetic energy contribution of a nucleonic wave function by the sum over Hesse matrix elements, which we have already calculated before. All in all, we can combine all energy contributions to obtain the total energy expectation value

E tot = ∑ 𝑖 ⟨𝜑 𝑖 | ĥ|𝜑 𝑖 ⟩ (2.78) = E Skyrme + E Coul + ∑ 𝑖 ⟨𝜑 𝑖 | T |𝜑 𝑖 ⟩ . (2.79)
The next step is to calculate its derivatives with respect to the gaussian parameters. For the potential part of the total energy we will make use of the effective potentials.

⟨𝑔 𝑛 |𝑈 pot |𝑔 𝑚 ⟩ = ∭ 𝑈 pot (ì 𝑥) 𝑔 𝑛𝑚 (ì 𝑥) d 3 𝑥 = ∭ √ 2𝜋 𝜒 𝑛𝑚,𝑖 √ 2𝜋 𝜒 𝑛𝑚,𝑖 √ 2𝜋 𝜒 𝑛𝑚,𝑗 √ 2𝜋 𝜒 𝑛𝑚,𝑗 √ 2𝜋 𝜒 𝑛𝑚,𝑘 √ 2𝜋 𝜒 𝑛𝑚,𝑘 𝑈 pot (ì 𝑥) 𝑔 𝑛𝑚 (ì 𝑥) d 3 𝑥 ⟨𝑔 𝑛 |𝑈 pot |𝑔 𝑚 ⟩ = ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ ∭ 1 √ 𝜋 3 𝑈 pot ( √ 2 ì 𝜒 𝑛𝑚,𝑖 ì 𝑥 ′ 𝑖 + ì 𝑋 𝑛𝑚,𝑖 ) 𝑒 -ì 𝑥 ′ 2 d 3 𝑥 ′ ≡ ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ 𝐼 0 (2.85) ⟨𝑔 𝑛 |(𝑥 𝑖 -𝑋 𝑛𝑚,𝑖 ) 𝑈 pot |𝑔 𝑚 ⟩ = ∭ (𝑥 𝑖 -𝑋 𝑛𝑚,𝑖 ) 𝑈 pot (ì 𝑥) 𝑔 𝑛𝑚 (ì 𝑥) d 3 𝑥 = ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ ∭ 1 √ 𝜋 3 √ 2𝜒 𝑛𝑚,𝑖 𝑥 ′ 𝑖 𝑈 pot ( √ 2 ì 𝜒 𝑛𝑚,𝑖 ì 𝑥 ′ 𝑖 + ì 𝑋 𝑛𝑚,𝑖 ) 𝑒 -ì 𝑥 ′ 2 d 3 𝑥 ′ ≡ ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ 𝐼 1,𝑖 (2.86) ⟨𝑔 𝑛 |(𝑥 𝑖 -𝑋 𝑛𝑚,𝑖 ) 2 𝑈 pot |𝑔 𝑚 ⟩ = ∭ (𝑥 𝑖 -𝑋 𝑛𝑚,𝑖 ) 2 𝑈 pot (ì 𝑥) 𝑔 𝑛𝑚 (ì 𝑥) d 3 𝑥 = ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ ∭ 1 √ 𝜋 3 2𝜒 𝑛𝑚,𝑖 𝑥 ′ 𝑖 2 𝑈 pot ( √ 2 ì 𝜒 𝑛𝑚,𝑖 ì 𝑥 ′ 𝑖 + ì 𝑋 𝑛𝑚,𝑖 ) 𝑒 -ì 𝑥 ′ 2 d 3 𝑥 ′ ≡ ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ 𝐼 2,𝑖 (2.87)
Hence, we define the integrals 𝐼 𝑛,𝑖 as

𝐼 𝑛,𝑖 = 1 √ 𝜋 3 ∭ √ 2𝜒 𝑛𝑚,𝑖 𝑥 𝑖 𝑛 𝑈 pot ( √ 2 ì 𝜒 𝑛𝑚,𝑖 ì 𝑥 𝑖 + ì 𝑋 𝑛𝑚,𝑖 ) 𝑒 -ì 𝑥 2 d 3 𝑥 (2.88)
where 𝑛 defines the order and 𝑖 defines the dimension in which the order of the moment is in. As a consequence, we can list each derivative of the potential energy

ì ∇ 𝑘 𝑚 ,𝑖 E pot = ∑ 𝑛 ℜ 𝑔 𝑛 𝑈 pot ì ∇ 𝑘 𝑚 ,𝑖 𝑔 𝑚 = ∑ 𝑛 ℜ 𝑖 𝐼 1,𝑖 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ + 𝐼 0 𝑔 𝑛 ì ∇ 𝑘 𝑚 ,𝑖 𝑔 𝑚 (2.89) ì ∇ 𝛾 𝑚 ,𝑖 E pot = ∑ 𝑛 ℜ 𝑔 𝑛 𝑈 pot ì ∇ 𝛾 𝑚 ,𝑖 𝑔 𝑚 = ∑ 𝑛 ℜ 𝑖 𝐼 2,𝑖 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ + 2𝐼 1,𝑖 𝑔 𝑛 ì ∇ 𝑘 𝑚 ,𝑖 𝑔 𝑚 + 𝐼 0 ⟨𝑔 𝑛 | ì ∇ 𝛾 𝑚 ,𝑖 -𝑖 𝜒 𝑛𝑚,𝑖 |𝑔 𝑚 ⟩ (2.90) ì ∇ 𝑥 𝑚 ,𝑖 E pot = -𝑖 ( ì 𝑘 𝑚,𝑖 + ì 𝜉 𝑚.𝑖 ì ∇ 𝑘 𝑚 ,𝑖 )E pot (2.91) ì ∇ 𝜒 𝑚 ,𝑖 E pot = - 1 4 ì 𝜒 𝑚,𝑖 (1 + 2𝑖 ì 𝜉 𝑚,𝑖 ì ∇ 𝛾 𝑚 ,𝑖 )E pot (2.92) 𝜕 𝑐 𝑚 E pot = 1 𝑐 𝑚 E pot (2.93) 𝜕 𝜃 𝑚 E pot = 𝑖 E pot (2.94)
At last, we need to take the derivative of the kinetic energy which can be written as

𝜕 𝑖 ⟨𝜑 𝑖 | T |𝜑 𝑖 ⟩ 𝜕𝑞 𝑚 = 𝜕 𝑞 𝑚 ⟨𝜑 𝑖 | T |𝜑 𝑖 ⟩ = ℜ 𝜑 𝑖 T 𝜕 𝑞 𝑚 𝜑 𝑖 (2.95) = - h 2𝑚 ℜ 𝜑 𝑖 ( ì ∇) 2 𝜕 𝑞 𝑚 𝜑 𝑖 = - h 2𝑚 ℜ 𝜑 𝑖 𝜕 2 𝑥 + 𝜕 2 𝑦 + 𝜕 2 𝑧 𝜕 𝑞 𝑚 𝜑 𝑖 .
(2.96)

In eq. 2.96 there are two types of expectation values to calculate. First case is when the spatial derivative is acting in a different dimension than the derivative with respect to the parameter 𝑞 𝑚 . As a result, the derivatives will factorize into

𝜑 𝑖 𝜕 2 𝑥 𝜕 𝑞 𝑚,𝑦 𝜑 𝑖 = ∑ 𝑛 𝑔 𝑛 𝜕 2 𝑥 𝜕 𝑞 𝑚,𝑦 𝑔 𝑚 = ∑ 𝑛 ⟨𝑔 𝑛 |𝜕 2 𝑥 |𝑔 𝑚 ⟩ 𝑔 𝑛 𝜕 𝑞 𝑚,𝑦 𝑔 𝑚 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ (2.97) = - ∑ 𝑛 𝜕 𝑥 𝑛 𝑔 𝑛 𝜕 𝑥 𝑚 𝑔 𝑚 𝑔 𝑛 𝜕 𝑞 𝑚,𝑦 𝑔 𝑚 ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ . (2.98)
The first factor of the numerator can be expressed again as the negative Hesse matrix element of the position derivatives similar to eq. 2.77 and the second factor is the expectation value of the Jacobian as in eq. 2.40 to eq. 2.44. In the second case, the spatial derivative acts in the same dimension as the parametric derivative, hence the derivatives do not factorize. We will move the double derivative 𝜕 2 𝑥 to the left-hand side of the braket using integration by parts and insert the explicit expression of 𝜕 2 𝑥 𝑔 𝑛 into the expectation value

𝜕 2 𝑥 𝑔 𝑛 = [ -(𝜉 𝑛 + 𝑘 2 𝑛 ) + 𝜉 2 𝑛 (𝑥 -𝑥 𝑛 ) 2 -2𝑖 𝜉 𝑛 𝑘 𝑛 (𝑥 -𝑥 𝑛 )]𝑔 𝑛 = [ -(𝜉 𝑛 + 𝑘 2 𝑛 ) -𝑖 𝜉 2 𝑛 𝜕 𝛾 𝑛 -2𝜉 𝑛 𝑘 𝑛 𝜕 𝑘 𝑛 ]𝑔 𝑛 (2.99) ⇒ 𝜑 𝑖 𝜕 2 𝑥 𝜕 𝑞 𝑚,𝑥 𝜑 𝑖 = ∑ 𝑛 𝑔 𝑛 𝜕 2 𝑥 𝜕 𝑞 𝑚,𝑥 𝑔 𝑚 = ∑ 𝑛 𝜕 2 𝑥 𝑔 𝑛 𝜕 𝑞 𝑚,𝑥 𝑔 𝑚 = ∑ 𝑛 [ -(𝜉 𝑛 + 𝑘 2 𝑛 ) -𝑖 𝜉 2 𝑛 𝜕 𝛾 𝑛 -2𝜉 𝑛 𝑘 𝑛 𝜕 𝑘 𝑛 ]𝑔 𝑛 𝜕 𝑞 𝑚,𝑥 𝑔 𝑚 = ∑ 𝑛 -(𝜉 𝑛 + 𝑘 2 𝑛 ) 𝑔 𝑛 𝜕 𝑞 𝑚,𝑥 𝑔 𝑚 + 𝑖 𝜉 𝑛 2 𝜕 𝛾 𝑛 𝑔 𝑛 𝜕 𝑞 𝑚,𝑥 𝑔 𝑚 -2𝜉 𝑛 𝑘 𝑛 𝜕 𝑘 𝑛 𝑔 𝑛 𝜕 𝑞 𝑚,𝑥 𝑔 𝑚 (2.100)
where the last equation eq. 2.100 is solely described in terms of the calculated Hesse matrix elements.

Finally, we have derived in this subsection as well as in the previous subsection 2.1.2 all the essential pieces that are required to solve the Euler-Lagrange equations in eq. 2.11, namely each Hesse matrix element along with every derivative of the total energy with respect to the parameters. Before we continue with the evolution, we have to discuss in the next section, how we construct an appropriate initialization of the nucleus, in order to start the evolution.

Initialization of a nucleus

We recall from the beginning of the chapter that the function space of a nucleonic wave function is spanned by a sum over any number of the parameterized gaussian basis functions. Hence, if there were no constraints the nucleonic wave function could take any functional form. However, in the construction of the density operator in eq. 2.5 we made the assumption that the system of singleparticle wave functions have to be orthogonal. Therefore, the function space reduces for each nucleonic wave function that we add to a system, since it has to be orthogonal to each other wave function already in that set. This property is important to fulfill, especially later during the time-evolution, since it ensures the fermionic nature of the system, that means the Pauli-principle. In non-orthogonal quantum models, such as in anti-symmetrized molecular dynamics models, one needs to include the interference terms among single-particle wave functions into the definition of the density operator, as shown in eq. 2.4, to account for the Pauli-principle. We will present in the following a self-consistent procedure that finds a suitable set of harmonic wave functions, which are naturally orthogonal to each other.

Self-consistent harmonic oscillator routine

The goal of the initialization process is to find a set of wave functions which is at equilibrium in its mean-field potential. For that reason, we want to restrict the function space to the set of harmonic oscillator wave functions, as they are naturally orthogonal to each other, hence we do not have to worry about an orthogonalization procedure. To find harmonic oscillator states, we need to know the meanfield potential of nucleons. We choose for the initial effective nuclear potential the density-independent

Woods-Saxon potential [Woods1954] 𝑈 initial WS (𝑟 ) = 𝑉 0 1 + exp 𝑟 -𝑅 𝑑 (2.101)
where 𝑉 0 = -53.0 MeV represents the potential well depth, 𝑑 = 0.66 fm shows the surface thickness and 𝑅 = 𝑟 0 𝐴 1/3 is the nuclear radius where 𝑟 0 = 1.3 fm and 𝐴 is the mass number. Additionally, we have to include the Coulomb potential generated by the protons of the nucleus. Assuming an uniform distribution, it can be expressed as

𝑈 initial Coul (𝑟 ) =        𝑒 2 4𝜋𝜖 0 1 𝑟 , 𝑟 > 𝑅 𝑒 2 4𝜋𝜖 0 1.5-0.5( 𝑟 𝑅 ) 2 𝑅 , 𝑟 ≤ 𝑅 , for protons.
(2.102) 1. The first step is fitting the effective total potentials by a three-dimensional harmonic oscillator potential

𝑈 eff 𝑛,𝑝 ≈ 𝑈 𝑛,𝑝 0 + 1 2 𝑚𝜔 2 𝑛,𝑝 ì 𝑥 2 . (2.105)
From the fit we obtain the parameters 𝑈 𝑛,𝑝 0 the depth of the effective potentials and 𝜔 𝑛,𝑝 the angular frequency of the harmonic oscillators. There will be one set of parameters for the neutrons and one for the protons, since we fit two different potentials. The full details of the fitting procedure will be discussed in the next section where we present the numerical implementation of the model. For now, the important fact is that the angular frequencies allow us to define a family of energy eigenstates of the quantum harmonic oscillator. In one dimension, they are given in terms of the Hermite functions 𝐻 𝑛 :

𝜑 𝑛 (𝑥) = 1 √ 2 𝑛 𝑛! 𝑚𝜔 𝜋h 1/4 𝑒 -𝑚𝜔𝑥 2 2h 𝐻 𝑛 √ 𝑚𝜔 h 𝑥 , 𝑛 = 0, 1, 2, . . . (2.106) 𝐻 𝑛 (𝑧) = (-1) 𝑛 𝑒 𝑧 2 𝑑 𝑛 𝑑𝑧 𝑛 𝑒 -𝑧 2 .
(2.107)

The three dimensional wave function can be constructed from a particular set of quantum numbers {𝑛 𝑥 , 𝑛 𝑦 , 𝑛 𝑧 } by the simple product of the corresponding one-dimensional eigenstates, since the Hamiltonian separates into independent one-dimensional harmonic oscillators

𝜑 (𝑛 𝑥 ,𝑛 𝑦 ,𝑛 𝑧 ) (ì 𝑥) = 𝜑 𝑛 𝑥 (𝑥)𝜑 𝑛 𝑦 (𝑦)𝜑 𝑛 𝑧 (𝑧) (2.108) 𝐸 𝑁 = h𝜔 𝑛 𝑥 + 𝑛 𝑦 + 𝑛 𝑧 + 3 2 ≡ h𝜔 𝑁 + 3 2 (2.109)
where 𝐸 𝑁 is the N-th energy level of the harmonic oscillator.

2. The second step of the routine is to build up the nucleus level by level starting from the ground state 𝑁 = 0. After that we go through the next excited energy level 𝑁 = 1 and assign for each combination of {𝑛 𝑥 , 𝑛 𝑦 , 𝑛 𝑧 } the corresponding harmonic wave function 𝜑 (𝑛 𝑥 ,𝑛 𝑦 ,𝑛 𝑧 ) such that

𝑛 𝑥 + 𝑛 𝑦 + 𝑛 𝑧 = 1.
For the first excited energy level, there are three different combination. In general, the degeneracy of a level N without including the spin degeneracy is given by

𝑔 𝑁 = (𝑁 + 1)(𝑁 + 2) 2 .
(2.110)

To take into account the spin multiplicity we double the number of states per level. When we complete an energy level 𝑁 with nucleons we continue with the next higher energy level 𝑁 + 1 as long as there are nucleons left to be defined. The process ends when there are no more nucleons to be constructed. The final result is that we end up with a family of harmonic wave functions which have been assembled from the bottom to the top.

3. Next, we can calculate the neutron and proton densities, since we have now defined the explicit distribution of each individual nucleonic wave function

𝜌 𝑛 (ì 𝑥) = 𝑁 ∑ 𝑖 𝜑 𝑛,𝑖 (ì 𝑥) 2 (2.111) 𝜌 𝑝 (ì 𝑥) = 𝑍 ∑ 𝑖 𝜑 𝑝,𝑖 (ì 𝑥) 2 (2.112)
where 𝑁 is the number of neutrons 𝑍 is the charge number. Consequently, we are able to calculate the density dependent effective potential 𝑈 eff pot

𝑈 eff pot = 𝑈 eff Skyrme + 𝑈 eff Coul (2.113)
where the effective Skyrme potential is given in eq. 2.66 and the effective Coulomb Potential in eq. 2.72. Note that the expression of the Skyrme potential includes terms with the second derivative of the density. We can perform the derivative analytically using the following relation for Hermite functions

d d𝑥 𝐻 𝑛 (𝑥) = 2𝑛 𝐻 𝑛-1 (𝑥) .
(2.114)

Thus, the derivatives of the harmonic oscillator states can be expressed by a combination of lower levels

d d𝑥 𝜑 𝑛 (𝑥) = -𝛽 2 𝑥 𝜑 𝑛 (𝑥) + 𝛽 √ 2𝑛 𝜑 𝑛-1 (𝑥) (2.115) d 2 d𝑥 2 𝜑 𝑛 (𝑥) = 𝛽 2 (𝛽 2 𝑥 2 -1) 𝜑 𝑛 (𝑥) -2𝛽 2 √ 2𝑛𝑥 𝜑 𝑛-1 (𝑥) + √ 4𝑛(𝑛 -1) 𝜑 𝑛-2 (𝑥) (2.116) d 2 d𝑥 2 𝜑 2 𝑛 (𝑥) = 2𝜑 𝑛 (𝑥) d 2 d𝑥 2 𝜑 𝑛 (𝑥) + 2 d d𝑥 𝜑 𝑛 (𝑥) 2 (2.117)
where the harmonic beta value is defined as

𝛽 = √ 𝑚𝜔 h . (2.118)
Finally, with the newly obtained effective potential we can repeat the self-consistent routine from step 1, where we fit an harmonic potential to the effective one. The iteration ends, when the total energy defined as the sum over the energy levels of the harmonic oscillator converges to a given precision 𝜀 > 0

𝐸 (𝑛) tot = 𝐴 ∑ 𝑖 𝐸 (𝑛) 𝑖 (2.119) 𝐸 (𝑛) tot -𝐸 (𝑛-1) tot ≤ 𝜀 (2.120)
where the superscript 𝑛 denotes the n-th iteration of the routine.

Cooling procedure

So far we have constructed a set of orthogonal wave function which are self-consistent to the effective potential that they create. However, we have restricted the wave functions so far to only harmonic oscillator eigenstates which are only a subset of all orthogonal states. Nevertheless, they provide a good initial guess of the nucleus to start a minimization procedure. For that reason we will represent each harmonic wave function as a superposition of our dynamical gaussian basis function. The basis function are dynamical in the sense that the parameters of the gaussians are time-dependent, hence the components can be moved similar to the gaussians in Molecular Dynamics models. The difference is that we use multiple gaussians to represent the de-localized nucleonic wave functions rather than just one. In order to approximate the harmonic eigenstates as close as possible, we fit each level with a different number of gaussian depending on the excitation. For example, the ground state of the harmonic oscillator is simply a gaussian, therefore we only have to fit with one gaussian. However, the second excitation in one dimension, see fig. 2.5, has two maxima and one minimum, thus we fit the wave function with three gaussians, one for every extrema. We will use the same principle for every other level as well, where we identify the number of extrema and fit that many gaussians. After we have changed our system to sets of dynamical gaussians per nucleon we can start our cooling procedure.

The minimization is based on the imaginary time evolution [Ono2004, Kanada1995] where we find that the nucleonic wave functions will evolve along the negative gradient of the total energy in imaginary time. Thus by using the same matrix elements as we have calculated in the previous subsection, we can use the prescription of the imaginary time-evolution equations in eq. 2.122

to evolve the dynamical gaussians to a further minimum which is closer to the real solution than the minimum which we obtained by restricting only to the sub-manifold of harmonic oscillator functions.

Conservation laws during the dynamical evolution

It is important to check, that the evolution equations indeed satisfy the conservation of different physical quantities. For instance, it is crucial that nucleonic wave functions stay orthogonal at all times since the full derivation of the evolution equation was based on that condition. Hence we perform the time-derivative of the scalar product between two states

𝜕 𝑡 𝜑 𝑖 𝜑 𝑗 = 𝜕 𝑡 𝜑 𝑖 𝜑 𝑗 + 𝜑 𝑖 𝜕 𝑡 𝜑 𝑗 = ĥ 𝑖 h 𝜑 𝑖 𝜑 𝑗 + 𝜑 𝑖 ĥ 𝑖 h 𝜑 𝑗 = - 1 𝑖 h ( 𝜑 𝑖 ĥ † 𝜑 𝑗 -𝜑 𝑖 ĥ 𝜑 𝑗 ) = 0 (2.123)
where we used in the last step, that the Hamiltonian is hermitian ĥ † = ĥ. This shows that the scalar product remains constant over time

𝜑 𝑗 𝜑 𝑖 (𝑡) = 𝜑 𝑗 𝜑 𝑖 (0) = 𝑂 𝑖 𝑗 .
(2.124)

Therefore the conservation of orthogonality is ensured. In addition, it also means that the normalization ⟨𝜑 𝑖 |𝜑 𝑖 ⟩ stays constant, which implies that the mass number of the system is conserved as well

𝐴 ∑ 𝑖 ⟨𝜑 𝑖 |𝜑 𝑖 ⟩ (𝑡) = 𝐴 = const. . (2.125)
The conservation of energy can be derived in the same way where we check the time-derivative of the energy expectation value instead of the scalar product 𝜕 𝑡 𝜑 𝑖 ĥ 𝜑 𝑗 = 𝜕 𝑡 𝜑 𝑖 ĥ 𝜑 𝑗 + 𝜑 𝑖 ĥ 𝜕 𝑡 𝜑 𝑗 = ĥ 𝑖 h 𝜑 𝑖 ĥ 𝜑 𝑗 + 𝜑 𝑖 ĥ ĥ

𝑖 h 𝜑 𝑗 = - 1 𝑖 h ( 𝜑 𝑖 ĥ2 𝜑 𝑗 -𝜑 𝑖 ĥ2 𝜑 𝑗 ) = 0 . (2.126)
In fact, it is also possible to see energy and normalization conservation directly in the equation of motion. Let us therefore look at the equation related to the derivative with respect to the phase

𝜃 𝑚 2hℑ 𝜕𝜑 𝑖 𝜕𝑡 𝜕𝑔 𝑚 𝜕𝜃 𝑚 = 𝜕E tot 𝜕𝜃 𝑚 (2.127) 2hℜ 𝜕𝜑 𝑖 𝜕𝑡 𝑔 𝑚 = 2ℜ 𝜑 𝑖 ĥ 𝜕𝑔 𝑚 𝜕𝜃 𝑚 = 2ℑ ⟨𝜑 𝑖 | ĥ|𝑔 𝑚 ⟩ (2.128) ⇐⇒ ∑ 𝑚 2hℜ 𝜕𝜑 𝑖 𝜕𝑡 𝑔 𝑚 = hℜ 𝜕𝜑 𝑖 𝜕𝑡 𝜑 𝑖 = ℑ ⟨𝜑 𝑖 | ĥ|𝜑 𝑖 ⟩ = 0 (2.129) =⇒ 𝜕 𝑡 ⟨𝜑 𝑖 |𝜑 𝑖 ⟩ = 0 (2.130)
where in the second line we used 𝜕 𝜃 𝑔 𝑚 = 𝑖 𝑔 𝑚 and that the imaginary unit transforms the Imaginary operator ℑ to the real operator ℜ and vice versa. In the third line we summed up all equations of eq. 2.128 to recover the total nucleonic wave function 𝜑 𝑖 . Then, we recall that the hamiltonian is hermitian, thus the imaginary part of the expectation value is zero. Finally we find that the normalization of the nucleonic wave function 𝜑 𝑖 is conservered over time. That means the introduction of the phase in the definition of the Gaussian allowed us to explicitly handle normalization conservation in the equation of motion. Similarly we find for the derivation with respect to the coefficient

𝑐 𝑚 2hℑ 𝜕𝜑 𝑖 𝜕𝑡 𝜕𝑔 𝑚 𝜕𝑐 𝑚 = 𝜕E tot 𝜕𝑐 𝑚 (2.131) 2h 𝑐 𝑚 ℑ 𝜕𝜑 𝑖 𝜕𝑡 𝑔 𝑚 = 2ℜ 𝜑 𝑖 ĥ 𝜕𝑔 𝑚 𝜕𝑐 𝑚 = 2 𝑐 𝑚 ℜ ⟨𝜑 𝑖 | ĥ|𝑔 𝑚 ⟩ (2.132) ⇐⇒ hℑ 𝜕𝜑 𝑖 𝜕𝑡 𝜑 𝑖 = ℜ ⟨𝜑 𝑖 | ĥ|𝜑 𝑖 ⟩ (2.133) ℜ 𝑖 h 𝜕𝜑 𝑖 𝜕𝑡 𝜑 𝑖 = ℜ ⟨𝜑 𝑖 | ĥ|𝜑 𝑖 ⟩ (2.134)
which ensures in our equation of motion that the evolution follows explicitly the Schrödinger equation.

This property was important to show the energy conservation in eq. 2.126.

Summary

In this chapter we have derived in full detail the theoretical foundation of our model. We started with a single-Slater approximation and added additionally an orthogonality condition to the set of singleparticle wave functions. This allowed us to drop the mixing terms in the density operator in eq. 2.5.

As a result the Lagrangian separated to a sum of Lagrangians that only included the behaviour of one nucleonic wave function. Therefore, the equation of motion, derived by the Euler-Lagrange equations, are only driven by the parameters that describe that particular single-particle wave function.

Specifically, the nucleonic wave functions are decomposed into a superposition of gaussians 𝑔 𝑚 each characterized by a position ì 𝑥 𝑚 , angular wavenumber ì 𝑘 𝑚 , spatial variance ì 𝜒 𝑚 , phase-space uncertainty ì 𝛾 𝑚 and complex weight 𝑐 𝑚 𝑒 𝑖𝜃 𝑚 parameter. For that choice of functional form, we continued to derive the Jacobian and Hesse matrix elements, which appear in the coefficients of the evolution matrix. Furthermore, we introduced the effective Skyrme-type potential and the Coulomb potential, since we had to calculate energy expectation values, which also governs the time evolution. Afterwards, we showed how we initialize the nucleus in our system. We start from an initial guess of the nuclear potential, the Woods-Saxon potential, and started a self-consistent routine, which approximates the nuclear potential with an harmonic oscillator in order to obtain a family of orthogonal energy eigenstates. Those eigenstates form a density distribution which in turn defines the previously mentioned effective Skyrme and Coulomb potential. The routine repeats iterates in this loop, where the new potential is fitted again, until the total energy of the harmonic oscillator converges. At last, we eventually decompose the de-localized harmonic wave functions into gaussians, where the number of gaussians depend on the level of excitation of the harmonic oscillator, and apply to the gaussian basis function a cooling procedure, based on the imaginary time evolution. Thus, we should be able to find a static solution which is close to the true minimum of the effective potential. Finally, we showed as well that the orthogonality condition, the mass number and the total energy of the system is conserved throughout the evolution at all times. Specifically, the introduction of the variable coefficient and the phase introduced into our equation of motion directly the conservation equation, which have to be solved at the same time as the evolution equations for the other parameters.

Chapter 3 Numerical Implementation of the model

In the previous chapter we have introduced a new approach to solve the mean-field evolution. Those mean-field models are considered in a first approximation and build a foundation for beyond-mean-field approaches [Simenel2010, Lacroix2004]. The main motivation was to develop a theory which generalizes the BUU-like models [Aichelin1985, Bertsch1988, Cassing1990, Bonasera1994] in a quantum framework. Thus, we remained in the TDHF description instead of its semi-classical analogue. However, different to traditional TDHF models [Balian1992, Ayik2008] where the Hartree-Fock states are described by the evolution of coefficients on a fixed grid, we represent the single-particle states by a set of parameterized basis functions. It is similar to the test-particle method in semi-classical models, only that we do not describe the phase-space distribution by a linear combination of dynamical functions but rather express the de-localized nucleonic wave functions by a quantum superposition of Gaussians. In a sense, the single-particle states become de-localized through the fact that we combine different localized basis functions together. Moreover, we can find an interesting observation if we take the limit to only one Gaussian per nucleon. As a result the single-particle wave function will be localized and not orthogonal to each other anymore. Thus, the general form of the density in eq. 2.4 with the inverse overlap matrix 𝑂 𝑖 𝑗 has to be used to ensure Pauli-principle. Therefore, we find that through those two changes we can make a connection from our approach to molecular dynamics models for fermions such as FMD or AMD [Feldmeier1990, Ono1992]. Essentially, we transitioned from a TDHF theory to a generalized molecular dynamics model with de-localized wave functions. We will discuss the connection between different types of transport models in the next section in more detail. It shows that our approach has many elements of different models, which is why we have to present the numerical framework of our model in the following section to clearly identify the advantages of our model. The common feature of all models is the many-body Schrödinger equation which determines the microscopic evolution. The general solution is given by a superposition of Slater determinants. However, practical models start from the single-Slater approximation. Then, there are two paths depending on whether the single-particle states are orthogonal or not. Our model decomposes Hartree-Fock states into dynamical Gaussians. We can draw the connection to F/AMD if we go to the limit of one Gaussian per nucleon. Finally BUU-models are connected through the Wigner transformation.

Connection between microscopic models

In this section we will follow the schematic diagram shown in fig. 3.1 which summarizes the crucial steps to draw the important connection between different nuclear transport models. As we can see each microscopic model at low to intermediate energies starts from the 𝐴-body Schrödinger equation. The general solution would be given by a superposition of Slater-determinants that are constructed from an arbitrary complete set of single-particle wave function. In practice, most of the existing models start from the single-Slater approximation however. Then, two common types of models can be constructed. On the one hand, TDHF-like models which keep non-local, orthonormal states [Maruhn2013, Sagert2016] and on the other hand molecular dynamics models for fermions that constrain the states to local, non-orthogonal gaussian-type functions. As we have mentioned, in order to determine the Hartree-Fock states, many models decompose the wave function 𝜑 𝑖 onto a fixed set basis function 𝑔 𝑗 and adjust the coefficients 𝑐 𝑖 𝑗 accordingly

|𝜑 𝑖 ⟩ = 𝑁 ∑ 𝑐 𝑖 𝑗 𝑔 𝑗 . (3.1)
It is clear that the decomposition becomes more accurate as the number of basis functions 𝑁 increases. The minimum number of basis functions is also given by the number of nucleons 𝐴, since we need at least that many coefficients to define 𝐴 linear independent single-particle states. An example from computational chemistry is the so-called linear combination atomic orbital (LCAO) approximation [Sekino2013]. The assumption is that the number of molecular orbitals |𝜑 𝑖 ⟩ is equal to the number of atomic orbitals 𝑔 𝑗 included in the linear expansion. The coefficients 𝑐 𝑖 𝑗 are the weights of the contributions of the atomic orbitals to the molecular orbital. We extended this method so that the basis functions 𝑔 𝑗 are also time-dependent and thus dynamical. In this case, we can then draw the connection to molecular dynamics for fermions when we go to the lowest possible number 𝑁 = 𝐴 of dynamical gaussian-type basis functions and additionally reintroduce the general definition of the density operator with the overlap matrix. Finally, we can also connect any quantum model to the semi-classical approaches such as BUU-models by performing a Wigner transformation which will map the quantum operators to its phase-space equivalent. Especially, if we introduce a decoherence assumption between the dynamical Gaussians such that the density operator further simplifies into

ρ = 𝑁 Gau ∑ 𝑖 |𝑐 𝑖 | 2 |𝑔 𝑖 ⟩ ⟨𝑔 𝑖 | (3.2) the Wigner transformation [Weyl1927] 𝑓 (𝑥, 𝑘) = 1 𝜋 ∞ -∞ d𝑠 𝑒 -2𝑖𝑘𝑠 ⟨𝑥 + 𝑠 | ρ|𝑥 -𝑠⟩ (3.3)
acts on each single density operator |𝑔 𝑖 ⟩ ⟨𝑔 𝑖 | individually since it is a linear operation. Note that the Wigner transformation of a gaussian-type density operator will turn to a Gaussian g𝑖 in phase-space with the corresponding position 𝑥 𝑖 , momentum 𝑘 𝑖 and variances 𝜒 𝑖 and 𝜙 𝑖 . All in all, we recover

ρ = 𝑁 Gau ∑ 𝑖 |𝑐 𝑖 | 2 |𝑔 𝑖 ⟩ ⟨𝑔 𝑖 | Wigner --------------→ transformation 𝑓 (𝑥, 𝑘) = 𝑁 Gau ∑ 𝑖 𝜔 𝑖 g𝑖 (𝑥, 𝑘) (3.4)
which is the decomposition in semi-classical transport models of the phase-space distribution into gaussian-type basis functions where 𝜔 𝑖 is the phase-space weight and 𝑁 Gau is the number of Gaussians.

History of the numerical framework

The numerical code that we have been working with has been initially developed by Francois Sebille and Bruno Jouault in the mid 90s [Jouault1996, Jouault1998, DeLaMota2001, Eudes2013]. The model is called DYWAN which stands for Dynamical Wavelets in Nuclei. Therefore, the idea to decompose the nucleonic wave functions into a dynamical basis was already formed during that time. The scheme of using wavelets as the basis functions was a clever and innovative theory, since wavelets can be constructed in such a way that they are inherently orthogonal to each other [Ali2016]. Today, this technique has found great success in computational physics and nuclear chemistry [Anderson2019, Fosso2013, Hamada2011, Yanai2004]. However, technically they do not form a dynamical basis but rather a static basis which can be assimilated with a TDHF implementation on a lattice. As a result, the first version of DYWAN has made some strong approximations. First, the wavelets have been reduced to the most simple case of a wavelet, which is a Gaussian. Furthermore, the decoherence assumption has been adopted [DeLaMota2019]. The latter approximation resulted to the problem that the model essentially moved from a quantum mean-field model to a semi-classical, BUU-like model where quantum effects like Pauli-principle are treated in terms of density distributions. The peculiarity of DYWAN is that unlike any other Boltzmann model, the test particles had a time-dependent width parameter and a time-constant weight which was different for each gaussian. The fact, that the weights were different originated from the initial treatment of de-localized wave-functions. However, due to the decoherence the development of DYWAN evolved into a model which has attempted to, but did not completely implement the de-localization of wave functions. Years later, the same model was applied in periodic conditions for stellar matter [Sebille2011] and more recently revisited by Gregoire

Besse for his PhD [Besse2017]. He made improvements on the computational side of the model, such as updating the effective Skyrme interaction and initializing the phase-space density more accurately.

But the model still included the decoherence assumption which has kept the calculation closer to a semi-classical approach rather than a quantum mean-field theory. The difference was that the initialization appeared to be more formal than in semi-classical approaches which initialize with a basic Thomas-Fermi model [Thomas1927], but in the end a Wigner transformation was applied to the quantum operators nevertheless. Hence, the mean-field evolution equation did not differ with respect to the Vlasov equation [Vlasov1938, Vlasov1968]. Later on a newer version implemented two-body collisions without however considering the nucleonic degree of freedom [Besse2020]. Therefore, by simply adding new features to the code, it ultimately would not improve the fundamental dynamics because the approach still behaves under the semi-classical approximation. Therefore, we decided to rework the approach so that it does not pass through a Wigner transformation anymore in order to prevent the transition to the decoherence assumption in the model. Our goal is to reduce the level of approximation and develop a new framework which keeps the detail of a quantum mean-field TDHF description, but establishes additionally a foundation to include beyond mean-field extensions of largeamplitude fluctuations. Consequently, we worked out the current theoretical approach which is why we went so thoroughly through the derivation of the model in the previous chapter. We want to make sure that there are no strong approximations which would spoil the quantum behavior of our theory as it has been for so long in the past. As a result, we kept the numerical framework that was started by Francois Sebille and which was developed over many years but restructured the code so that the original idea of de-localized wave functions are well handled. Instead, we use the general structure that has been provided and modified the crucial parts that were associated to the mean-field evolution. To be specific, the code is separated into two parts, the first one calculates the ground state of a nucleus and stores the list of gaussian parameters into an initialization data file. The second part reads the data from the initialization files as well as information about the kinematics and starts a dynamical evolution of the nuclear reaction under those conditions. In the following we will present the numerical implementation of both routines.

Implementation of the statics

The goal of the first part of the code is to initialize a nucleus according to the steps we discussed in subsection 2.2. Note, that numerically we solve the problem in one dimension first and construct the three-dimensional nucleus subsequently. Thus, we assume that the effective mean-field potential is radially symmetric and that we are able to fit it with a radially symmetric harmonic potential. In the following we discuss how we fit an harmonic oscillator to the effective potentials of the neutrons and protons and show how we decompose the associated harmonic oscillator energy eigenstates into our dynamical gaussian-type basis function in three dimensions.

Harmonic oscillator fitting procedure

We recall that the first step of the self-consistent routine is to define the two initial potentials 𝑈 initial n and 𝑈 initial p in eq. 2.103 and eq. 2.104 where the only difference is the additional Coulomb potential for the proton potential. In the code, we are working on a fine lattice from -20 fm to 20 fm with a lattice spacing of 0.04 fm. The potentials are given on that lattice. The next procedure is to fit those potentials by an harmonic oscillator. The fit is not performed immediately but we prepare an approximated potential first to improve the fitting process. To be specific, we extend the sides of the We see that the approximation extends the sides of the effective potential. In-between, the smoothing is exact since the graph does not need to be approximated. The Woods-Saxon potential has already a positive curvature at every point within that interval. Therefore, the curve overlaps with the initial potential as we can see in the left plot. Finally, we perform a simple minimization with the fitting function

Û = 𝑈 0 + Ω𝑥 2 (3.5)
where 𝑈 0 and Ω are the parameters of the fit. On the right-hand side, we see the same process after the first iteration. If we remember, the harmonic oscillator defines a set of eigenstates which in turn gives the spatial distribution through the sum of their complex squares. Thus, we can determine the density-dependent effective Skyrme potential for the second iteration. The Skyrme potential has a more complicated structure than the initial Woods-Saxon potential, hence the smoothing approximation is needed. For the Skyrme parameterization there are many different sets of parameters. Among the most popular parameterization that are used for the TDHF code SKY3D are the SkM [Bartel1982] and the SLy4 parameters [Chabanat1997]. In our calculations we used mainly the SV-bas parameterization [Klupfel2009] and reduced it to the 𝑡 0 -and 𝑡 3 -term in order to compare the results with previous calculations. Note that we do not refit the 𝑡 0 -and 𝑡 3 -parameters, although we do not include the 𝑡 1 -and 𝑡 2 -contributions to the effective mass. Therefore, the properties of infinite nuclear matter are slightly changed, namely the saturation density, the energy per nucleon at saturation and the symmetry energy become 𝜌 0 ≈ 0.188 fm -3 , 𝐸/𝐴 = -18.5 MeV and 𝑎 sym = 43.5 MeV, respectively, while the full SV-bas parameterization gives 𝜌 0 = 0.1596 fm -3 , 𝐸/𝐴 = -15.9 MeV and 𝑎 sym = 30 MeV [Klupfel2009].

We can see in the right plot that the repulsive center of the effective potential lifts the approximation which results to an overall higher potential well. In fig. the effective potential such that the harmonic fit function gives a different angular frequency for proton wave functions compared to neutron wave functions. Also we repeat the harmonic fits for the initial and second iterations for a heavier nucleus 124 Sn in fig. 3.4 and fig. 3.5 . In that case we observe that the initial potential is much wider since the nucleus is bigger. Therefore, the harmonic oscillator fit has a deeper initial potential well. In fig. 3.4 the potential well of the harmonic oscillator fit in the This becomes even more apparent if we look at the total evolution of the effective potential during the self-consistency routine in fig. 3.6 for the 48 Ca and fig. 3.7 for the 124 Sn initialization. In both figures the left plot shows the full evolution of the effective neutron potential and the right plot the evolution of the proton potential. We represent again in solid black the initial potential for the respective cases.

The dashed black lines are potentials at intermediate steps of the precedure and finally the solid blue line shows the final effective potential, when the self-consistency routine converged to the condition in eq. 2.120 where we choose 𝜀 = 0.5 MeV. Comparing the two figures, we deduce that the restriction of the Hartree-Fock states to the subspace of harmonic oscillator wave functions is better suited for heavier nuclei. The reason for that is, that the structural effect are gradually washed out as the nuclei gets heavier. Thus, the system can be sufficiently well described on a restricted sub-manifold as long as the details of the individual components do not play a significant role. Conversely, for lighter nuclei the structural properties of harmonic oscillator states certainly differ from the actual Hartree-Fock states.

That is why, we add another routine to the initialization which helps us to get closer to that solution by cooling the states. Before that, we first need to discuss how we decompose the one-dimensional harmonic wave functions into parameterized Gaussians and use them construct the three-dimensional nucleus.

Decomposition into dynamical Gaussians

From the previous subsection, we obtain from the fit of the harmonic oscillator potential the angular frequency 𝜔 which in turn defines us the harmonic beta value 𝛽 = √ 𝑚𝜔 h , where 𝑚 = 931.5 MeV is the atomic mass, to define the harmonic oscillator states 𝜑 𝑛 (𝑥) according to eq. 2.106. As mentioned in subsection 2.2.2 we decompose the harmonic oscillator states by fitting Gaussians equal to the number of peaks of the wave function. Therefore, the fitting function φ𝑛 (𝑥) reads

φ𝑛 (𝑥) = 𝑁 𝑛 ∑ 𝑖 𝑔 𝑖 (𝑥) (3.6)
where 𝑁 𝑛 = 𝑛 + 1 is the number of peaks of the harmonic oscillator state 𝑛 and 𝑔 𝑖 are the Gaussian paramaterization defined in eq. 2.14 with the parameters P 𝑖 ≡ {𝑘 𝑖 , 𝛾 𝑖 , 𝑥 𝑖 , 𝜒 𝑖 , 𝑐 𝑖 , 𝜃 𝑖 }. As an initial guess to fit the wave function we assign the position of the peaks to each Gaussian and determine the full width at half maximum (FWHM) for each peak individually as well. Using the relation

FWHM 𝑖 = 4 √ ln 2 𝜒 𝑖 =⇒ 𝜒 𝑖 = 1 4 ln 2 FWHM 𝑖 2 2 (3.7)
we can also initialize the variances using the full width at half maximum of each peak. Lastly, the coefficients 𝑐 𝑖 are given by the height of the peaks and all the remaining parameters 𝑘 𝑖 ,𝛾 𝑖 and 𝜃 𝑖 are ini-tialized by zero. All in all, we can start a minimization procedure. We choose a simple Gauss-Newton algorithm to solve non-linear least squares problem. The results of the fitting procedure are shown in fig. 3.8. As we can see in the first plot for 𝑛 = 0, the ground state of the harmonic oscillator is exactly described by a Gaussian. For the higher levels 𝑛 ≥ 1 the fitting procedure adjusts the position and variance of the Gaussians to match the harmonic oscillator states. In the case of 𝑛 = 1 we see that the position of the Gaussians are slightly shifted towards each other and the width as well as the amplitude is increased to account for the fact that left Gaussian has a negative amplitude while the right Gaussian has a positive amplitude. Thus, the two Gaussians annihilate each other in the overlap region. It is even more clear for 𝑛 = 3, where the two outer Gaussians have an amplitude close to the peaks amplitude but the central Gaussian needs to be increased to counteract the additive contribution of the other two. For the sake of clarity, we only show the decomposition for 𝑛 = 3 but not 𝑛 = 4, 5 since we do not obtain any new information out of them. From the plot for 𝑛 = 4 we notice that the central Gaussians are more narrow than the outer Gaussians which can be explained by the Heisenberg-uncertainty principle. Different to the ground state, these are Gaussians that belong to a higher level, so the variance in momentum space is larger than for a lower level such as the ground state. As a consequence, the spatial variance is lower if we assume minimal uncertainty.

Lastly, the decomposition is not fully completed yet. We treated the problem only in one dimension so far. However, the solution of the three-dimensional harmonic wave function factorizes, see eq. 2.108, where the wave function can be differently excited in each of the three dimensions. The total energy level is then given by the sum of the energy levels in each dimension. As discussed in the second step of the self-consistency routine in subsection 2.2.1, we build the entire nucleus using a bottom-up approach. Let us explain the method to assign the levels to the nucleons by the explicit example of the 124 Sn nucleus. The nucleus consists of 74 neutrons and 50 protons. With the degeneracy 𝑔 𝑁 in eq. 2.110 we can associate the maximum energy level for both isospin types. The accumulative sum of the degeneracy is given by

𝑁 max ∑ 𝑁 =0 𝑔 𝑁 = (𝑁 max + 1)(𝑁 max + 2)(𝑁 max + 3) 6 . (3.8)
Therefore, we can create a table for different maximum energy levels 𝑁 max . The tab. 3.1 shows the number of different three-dimensional harmonic oscillator states for that level and the number of total states up to that level. Furthermore, we double this number to include the spin degeneracy as well.

The reason is that we initialize the same state twice, where one corresponds to spin up and the other to spin down. Finally, we can build up the nucleus. We start with the 74 neutrons. From tab. 3.1 we

𝑁

𝑔 𝑁 and its double. The sum provides the number of different three-dimensional harmonic oscillator states that can be construct up to that level. We double the number of states due to the spin degeneracy.

Gaussians is then given by the different combinations of Gaussians of those wave functions

𝜑 (𝑛 𝑥 ,𝑛 𝑦 ,𝑛 𝑧 ) (ì 𝑥) = 𝜑 𝑛 𝑥 (𝑥) 𝜑 𝑛 𝑦 (𝑦) 𝜑 𝑛 𝑧 (𝑧) (3.9) = 𝑁 𝑥 ∑ 𝑖 𝑁 𝑦 ∑ 𝑗 𝑁 𝑧 ∑ 𝑘 𝑔 𝑛 𝑥 ,𝑖 (𝑥) 𝑔 𝑛 𝑦 ,𝑗 (𝑦) 𝑔 𝑛 𝑧 ,𝑘 (𝑧) (3.10) = Ñ ∑ 𝑖 𝑔 𝑖 (ì 𝑥) (3.11)
where 𝑁 𝑥 , 𝑁 𝑦 and 𝑁 𝑧 are the respective number of Gaussians and Ñ its product. The last line indicates that we decomposed the wave function into three-dimensional gaussian-type functions. Moreover, we make another copy of this set of Gaussians and attribute a different spin to it, in order to account for the spin degeneracy. The excess nucleons that have not been assigned yet, in this example the 74 -70 = 4 neutrons, are placed in the upper level 𝑁 max = 5 which is not completely filled. In the routine, they are attributed to a combination (𝑛 𝑥 , 𝑛 𝑦 , 𝑛 𝑧 ) such that 𝑛 𝑥 + 𝑛 𝑦 + 𝑛 𝑧 = 𝑁 max according to a specific sorting algorithm, which fills first the states which have the maximum number of gaussians (in the current example you start with (1,2,2) and go through (2,1,2), (2,2,1) and only then move to (1,1,3) and so on). The same decomposition are performed for the 50 protons, where we keep in mind that the harmonic wave functions are different since the potentials were different. Again, for 50 protons we see that we can fill up the levels up to 𝑁 = 3 and the remaining 50 -40 = 10 protons are distributed among the states of the maximum level 𝑁 max = 4.

Implementation of the cooling process

The last routine in the first part of the numerical code is related to the cooling procedure described in subsection 2.2.2. The evolution equations eq. 2.120 have been derived from the imaginary time evolution method [McArdle2019], which transformed our equation of motion in eq. 2.10 by replacing the imaginary part of the Hesse matrix elements by their real part and the derivatives of the energy expectation value by its negative. To understand the effects of the evolution it is instructive to investigate the single-Gaussian and the two-Gaussian evolution. We will see that a wave function that is localized to a single Gaussian, evolves like an independent particle in a common mean-field potential whereas a wave function that consists of multiple Gaussians carries additional interference effects in the evolution equation preventing them to overlap.

Cooling of a single-Gaussian wave function

We need to calculate the Hesse matrix elements ì ∇ 𝑞,𝑖 𝑔 1 ì ∇ 𝑝,𝑖 𝑔 1 . In the special case of a single Gaussian, the formulas presented in subsection 2.1.2 reduce to very simple expressions. The reason is that the product function 𝑔 11 has the properties of the single Gaussian

𝜒 11,𝑖 = 1 2 ℜ𝜉 1,𝑖 = 𝜒 1,𝑖 (3.12) 𝑋 11,𝑖 = 𝑥 1,𝑖 (3.13) 𝐶 11,𝑖 = - 𝑥 2 1,𝑖 2 𝜒 1,𝑖 (3.14)
thus the scalar product simplifies to

⟨𝑔 1 |𝑔 1 ⟩ 𝑖 = 1 =⇒ ⟨𝑔 1 |𝑔 1 ⟩ = 𝑐 2 1 . (3.15)
For the Jacobians we also find very simple expressions

𝑔 1 ì ∇ 𝑘 1 ,𝑖 𝑔 1 = 0 (3.16) 𝑔 1 ì ∇ 𝛾 1 ,𝑖 𝑔 1 = 𝑖 ì 𝜒 1,𝑖 ⟨𝑔 1 |𝑔 1 ⟩ (3.17) 𝑔 1 ì ∇ 𝑥 1 ,𝑖 𝑔 1 = -𝑖 ì 𝑘 1,𝑖 ⟨𝑔 1 |𝑔 1 ⟩ = 0 (3.18) 𝑔 1 ì ∇ 𝜒 1 ,𝑖 𝑔 1 = - 1 4 ì 𝜒 1,𝑖 ⟨𝑔 1 |(1 + 2𝑖 ì 𝜉 1,𝑖 ì ∇ 𝛾 1 ,𝑖 )|𝑔 1 ⟩ = 0 (3.19) 𝑔 1 𝜕 𝑐 1 𝑔 1 = 1 𝑐 1 ⟨𝑔 1 |𝑔 1 ⟩ , 𝑔 1 𝜕 𝜃 1 𝑔 1 = 𝑖 ⟨𝑔 1 |𝑔 1 ⟩ (3.20)
where we used our specific initialization ì 𝑘 1 = 0 and ì 𝛾 1 = 0, thus ì 𝜉 1,𝑖 = 1 2 𝜒 1,𝑖 . Notice only derivatives with respect to ì 𝛾 1 , 𝑐 1 or 𝜃 1 are non-vanishing. We also recall from eq. 2.51 that the matrix elements can be written in terms of an independent part and a mixing term

( Ĥ𝑞,𝑝 ) 𝑖 𝑗 ≡ ì ∇ 𝑞,𝑖 𝑔 1 ì ∇ 𝑝,𝑗 𝑔 1 = ì ∇ 𝑞,𝑖 𝑔 1 𝑔 1 𝑔 1 ì ∇ 𝑝,𝑗 𝑔 1 ⟨𝑔 1 |𝑔 1 ⟩ + 𝛿 𝑖 𝑗 𝑀 𝑞 𝑖 ,𝑝 𝑖 . (3.21)
The independent part consists of the product of two Jacobians. Hence, only the matrix elements that include any combination between ì 𝛾 1 , 𝑐 1 and 𝜃 1 have a non-zero independent-part contribution. The mixing terms are calculated individually. We obtain

𝑀 𝑘 1,𝑖 ,𝑘 1,𝑖 = ì 𝜒 1,𝑖 ⟨𝑔 1 |𝑔 1 ⟩ (3.22) 𝑀 𝛾 1,𝑖 ,𝛾 1,𝑖 = 2 ì 𝜒 2 1,𝑖 ⟨𝑔 1 |𝑔 1 ⟩ (3.23) 𝑀 𝑥 1,𝑖 ,𝑘 1,𝑖 = 𝑖 𝜉 1,𝑖 ì 𝜒 1,𝑖 ⟨𝑔 1 |𝑔 1 ⟩ = 𝑖 1 2 ⟨𝑔 1 |𝑔 1 ⟩ (3.24) 𝑀 𝑥 1,𝑖 ,𝑥 1,𝑖 = 𝜉 1,𝑖 2 𝜒 1,𝑖 1 2 ⟨𝑔 1 |𝑔 1 ⟩ = 1 4 𝜒 1,𝑖 ⟨𝑔 1 |𝑔 1 ⟩ (3.25) 𝑀 𝜒 1,𝑖 ,𝛾 1,𝑖 = 𝑖 𝜉 1,𝑖 2 ì 𝜒 1,𝑖 2 ì 𝜒 2 1,𝑖 ⟨𝑔 1 |𝑔 1 ⟩ = 𝑖 1 2 ⟨𝑔 1 |𝑔 1 ⟩ (3.26) 𝑀 𝜒 1,𝑖 ,𝜒 1,𝑖 = 𝜉 1,𝑖 2 ì 𝜒 1,𝑖 1 2 ⟨𝑔 1 |𝑔 1 ⟩ = 1 2 1 2 ì 𝜒 1,𝑖 2 ⟨𝑔 1 |𝑔 1 ⟩ (3.27)
𝑀 𝑥 1,𝑖 ,𝛾 1,𝑖 = 0 , 𝑀 𝛾 1,𝑖 ,𝑘 1,𝑖 = 0 (3.28)

𝑀 𝜒 𝑛,𝑖 ,𝑘 1,𝑖 = 0 , 𝑀 𝑥 1,𝑖 ,𝜒 1,𝑖 = 0 (3.29)

𝑀 𝑐 𝑛 ,𝑝 𝑖 = 0 , 𝑀 𝜃 𝑛 ,𝑝 𝑖 = 0 .

(3.30)

We deduce, the only real contributions come from Hesse matrix elements with symmetric derivatives, that means of the form ì ∇ 𝑝,𝑖 𝑔 1 ì ∇ 𝑝,𝑖 𝑔 1 except for one matrix element

𝜕 𝜃,𝑖 𝑔 1 ì ∇ 𝛾 1 ,𝑖 𝑔 1 = ì 𝜒 1,𝑖 ⟨𝑔 1 |𝑔 1 ⟩ . (3.31)
Hence, we can visualize the matrix as following

ĤCooling = Ĥ𝑘,𝑘 0 0 0 ì 0 ì 0 0 Ĥ𝛾,𝛾 0 0 ì 0 ì H 𝜃,𝛾 0 0 Ĥ𝑥,𝑥 0 ì 0 ì 0 0 0 0 Ĥ𝜒,𝜒 ì 0 ì 0 ì 0 ⊤ ì 0 ⊤ ì 0 ⊤ ì 0 ⊤ H 𝑐,𝑐 0 ì 0 ⊤ ì H ⊤ 𝜃,𝛾 ì 0 ⊤ ì 0 ⊤ 0 H 𝜃,𝜃 (3.32)
where ĤCooling = ℜ ì ∇ 𝑞 𝑔 1 ì ∇ 𝑝 𝑔 1 is the cooling evolution matrix, and its Ĥ𝑝,𝑝 the sub-matrices of the Hesse matrix where 𝑝 ∈ P 1 is a parameter of the Gaussian. The real part of the sub-matrices have the simple expressions

Ĥ𝑘,𝑘 ℜ = 𝜒 1,𝑥 0 0 0 𝜒 1,𝑦 0 0 0 𝜒 1,𝑧 ⟨𝑔 1 |𝑔 1 ⟩ (3.33) Ĥ𝛾,𝛾 ℜ = 3 𝜒 2 1,𝑥 𝜒 1,𝑦 𝜒 1,𝑥 𝜒 1,𝑧 𝜒 1,𝑥 𝜒 1,𝑥 𝜒 1,𝑦 3 𝜒 2 1,𝑦 𝜒 1,𝑧 𝜒 1,𝑦 𝜒 1,𝑥 𝜒 1,𝑧 𝜒 1,𝑦 𝜒 1,𝑧 3 𝜒 2 1,𝑧 ⟨𝑔 1 |𝑔 1 ⟩ (3.34) Ĥ𝑥,𝑥 ℜ = 1 2 1 2 𝜒 1,𝑥 0 0 0 1 2 𝜒 1,𝑦 0 0 0 1 2 𝜒 1,𝑧 ⟨𝑔 1 |𝑔 1 ⟩ (3.35) Ĥ𝜒,𝜒 ℜ = 1 2 1 2 𝜒 1,𝑥 2 0 0 0 1 2 𝜒 1,𝑦 2 0 0 0 1 2 𝜒 1,𝑧 2 ⟨𝑔 1 |𝑔 1 ⟩ (3.36) ì H 𝜃,𝛾 ℜ = ì 𝜒 1 ⟨𝑔 1 |𝑔 1 ⟩ , H 𝑐,𝑐 ℜ = 1 𝑐 2 1 ⟨𝑔 1 |𝑔 1 ⟩ , H 𝜃,𝜃 ℜ = ⟨𝑔 1 |𝑔 1 ⟩ . (3.37)
Secondly, we need to calculate the energy derivatives. The kinetic energy of a single-Gaussian wave function takes a simple form

E kin = ⟨𝑔 1 | T|𝑔 1 ⟩ = h 2𝑚 ( 𝜕 𝑥 1 𝑔 1 𝜕 𝑥 1 𝑔 1 + 𝜕 𝑦 1 𝑔 1 𝜕 𝑦 1 𝑔 1 + 𝜕 𝑧 1 𝑔 1 𝜕 𝑧 1 𝑔 1 ) (3.38) = h 2𝑚 ( ì 𝑘 2 1 + 𝜙 1,𝑥 + 𝜙 1,𝑦 + 𝜙 1,𝑧 )𝑐 2 1 = h 2𝑚 ( 1 4𝜒 1,𝑥 + 1 4𝜒 1,𝑦 + 1 4𝜒 1,𝑧 )𝑐 2 1 (3.39)
which only depends on ì 𝜒 1 and 𝑐 1 . Remember that 𝜙 1 = 1 4𝜒 1 is the momentum variance. That means, the kinetic energy is given by the momentum variance of the Gaussian. The energy contributions

h ì 𝑘 2 1
2𝑚 related to overall momentum vanishes, because ì 𝑘 1 = ì 0 was initialized. We summarize for the derivatives of the kinetic energy

ì ∇ 𝜒 1 ,𝑖 E kin = - h 2𝑚 ⟨𝑔 1 |𝑔 1 ⟩ 4 𝜒 2 1,𝑖 (3.40) 𝜕 𝑐 1 E kin = h 2𝑚 ⟨𝑔 1 |𝑔 1 ⟩ 2𝑐 1 𝜒 1,𝑖 (3.41)
𝜕 𝑝 E kin = 0 (3.42)

for 𝑝 ∈ P 1 \{ ì 𝜒 1 , 𝑐 1 }. Lastly, in order to determine the derivatives of the potential energy, we recall that the derivatives can be expressed as following

𝜕 𝑞 E pot = ∭ 𝑈 Skyrme 𝜕𝑔 11 𝜕𝑞 d 3 𝑥 (3.43)
according to eq. 2.81 with 𝑞 ∈ P 1 . As a result, the derivatives depend on the derivatives of the product

function 𝑔 11 𝑔 11,𝑖 (𝑥 𝑖 ) = 1 2𝜋 𝜒 1,𝑖 1/2 exp - (𝑥 𝑖 -𝑥 1,𝑖 ) 2 2𝜒 1,𝑖 (3.44) 𝑔 11 (ì 𝑥) = 𝑐 2 1 𝑔 11,𝑥 (𝑥) 𝑔 11,𝑦 (𝑦) 𝑔 11,𝑧 (𝑧) (3.45)
which only depends on ì 𝜒 1 , ì 𝑥 1 . Therefore, there is only a contribution for the position and variance parameter derivatives of the potential energy. All in all, the evolution equations are condensed to

2h Ĥ𝑘,𝑘 0 0 0 ì 0 ì 0 0 Ĥ𝛾,𝛾 0 0 ì 0 ì H 𝜃,𝛾 0 0 Ĥ𝑥,𝑥 0 ì 0 ì 0 0 0 0 Ĥ𝜒,𝜒 ì 0 ì 0 ì 0 ⊤ ì 0 ⊤ ì 0 ⊤ ì 0 ⊤ H 𝑐,𝑐 0 ì 0 ⊤ ì H ⊤ 𝜃,𝛾 ì 0 ⊤ ì 0 ⊤ 0 H 𝜃,𝜃 • ì 𝑘 1 ì 𝛾 1 ì 𝑥 1 ì 𝜒 1 𝑐 1 𝜃 1 = - ì 0 ì 0 ì ∇ 𝑥 1 E pot ì ∇ 𝜒 1 E kin + ì ∇ 𝜒 1 E pot 𝜕 𝑐 1 E kin + 𝜕 𝑐 1 E pot 0 . (3.46)
By inverting the evolution matrix ĤCooling we obtain the time-derivatives of the gaussian parameters.

We find

ì 𝑥 1,𝑖 = - 2 𝜒 1,𝑖 𝑐 2 1 h ì ∇ 𝑥 1 ,𝑖 E pot (3.47) ì 𝜒 1,𝑖 = - 4 𝜒 2 1,𝑖 𝑐 2 1 h ì ∇ 𝜒 1 ,𝑖 E tot (3.48) 𝑐 1 = -𝑐 1 E tot (3.49) ì 𝑘 1 = ì 0 , ì 𝛾 1 = ì 0 , 𝜃 1 = 0 . (3.50)
These equations are really insightful. They show that in the cooling procedure, the momentum does not change when we initialize them with zero. The same holds for the 𝛾 parameter related to the minimum uncertainty and the phase 𝜃 1 . The only parameters that change are the position and the variance, which follow the negative gradients of the total energy with respect to the corresponding parameter. That means, the Gaussian is displaced step-by-step into the well of the effective potential.

Meanwhile, the variance adjusts itself to take into account the change in kinetic and potential energy.

If the spatial variance increases, then the momentum variance decreases which would lower the total energy. However, if the Gaussian flattens, then the potential well is lifted which would increase the total energy. Thus the procedure minimizes the width such that it fits correctly into the potential well. Lastly, we have to deal with the weight evolution 𝑐 1 . Obviously, if we continue to increase the weight as it is shown in eq. 3.49 the total energy behaves accordingly. The reason is that the imaginary-time evolution equations are not particle number conserving, which is also why they are not energy conserving. But that is what we aimed for, to minimize the initial conditions. We solve the problem, by restricting the evolution on the sub-manifold of the parameter-space which keeps the normalization of the wave function but varies the total energy. Numerically we ensure that in the code by re-normalizing the wave function each step of the evolution.

Cooling of a two-Gaussian wave function

In the previous example, we have seen that the single-Gaussian evolution is solely determined by the effective mean-field potential. Now we discuss, which additional effect appear if we include more Gaussians to the wave function. Therefore we build the cooling evolution matrix again. In detail, the matrix separates into single-Gaussian evolution matrices

ĤCooling = ĤCooling

where ĤCooling 11 and ĤCooling

22

are the same evolution matrices as in eq. 3.32 for Gaussian 𝑔 1 and 𝑔 2 .

The off-diagonal evolution matrices

( Ĥ21 ) 𝑞,𝑝 = ì ∇ 𝑞 𝑔 2 ì ∇ 𝑝 𝑔 1 (3.52) Ĥ21 = Ĥ † 12 (3.53)
are the interference contributions between the two Gaussians. In order simplify the expression as we did it for the single-Gaussian case, we notice that matrix elements ( Ĥ21 ) 𝑘,𝑥 , ( Ĥ21 ) 𝑘,𝜒 , ( Ĥ21 ) 𝑘,𝑐 and ( Ĥ21 ) 𝛾,𝑥 , ( Ĥ21 ) 𝛾,𝜒 , ( Ĥ21 ) 𝛾,𝑐 are purely imaginary in our case, where we initialize by ì 𝑘 = ì 0 and ì 𝛾 = ì 0.

Hence, the real part of the off-diagonal evolution matrix can be written as

Ĥ21 = ( Ĥ21 ) 𝑘,𝑘 ( Ĥ21 ) 𝛾,𝑘 0 0 ì 0 ( ì H 21 ) 𝜃,𝑘 ( Ĥ21 ) 𝑘,𝛾 ( Ĥ21 ) 𝛾,𝛾 0 0 ì 0 ( ì H 21 ) 𝜃,𝛾 0 0 ( Ĥ21 ) 𝑥,𝑥 ( Ĥ21 ) 𝜒,𝑥 ( ì H 21 ) 𝑐,𝑥 ì 0 0 0 ( Ĥ21 ) 𝑥,𝜒 ( Ĥ21 ) 𝜒,𝜒 ( ì H 21 ) 𝑐,𝜒 ì 0 ì 0 ⊤ ì 0 ⊤ ( ì H ⊤ 21 ) 𝑥,𝑐 ( ì H ⊤ 21 ) 𝜒,𝑐 (H 21 ) 𝑐,𝑐 0 ( ì H
⊤ 21 ) 𝑘,𝜃 ( ì H ⊤ 21 ) 𝛾,𝜃 ì 0 ⊤ ì 0 ⊤ 0 (H 21 ) 𝜃,𝜃 . (3.54)
As we can see, the evolution matrix separates again. Together with the fact that the total energy does not depend on the parameters ì 𝑘, ì 𝛾 and 𝜃 as well, we arrive to the same situation as in the previous case, that the momentum, 𝛾 and phase parameter do not change if we initialize them with zero from the start. We can split the evolution equations into two sub-systems of equations

2h ( Ĥ11 ) 𝑘,𝑘 0 ì 0 ( Ĥ21 ) 𝑘,𝑘 ( Ĥ21 ) 𝛾,𝑘 ( ì H 21 ) 𝜃,𝑘 0 ( Ĥ11 ) 𝛾,𝛾 ( ì H 11 ) 𝜃,𝛾 ( Ĥ21 ) 𝑘,𝛾 ( Ĥ21 ) 𝛾,𝛾 ( ì H 21 ) 𝜃,𝛾 ì 0 ⊤ ( ì H ⊤ 11 ) 𝛾,𝜃 (H 11 ) 𝜃,𝜃 ( ì H ⊤ 21 ) 𝑘,𝜃 ( ì H ⊤ 21 ) 𝛾,𝜃 (H 21 ) 𝜃,𝜃 ( Ĥ12 ) 𝑘,𝑘 ( Ĥ12 ) 𝛾,𝑘 ( ì H 12 ) 𝜃,𝑘 ( Ĥ22 ) 𝑘,𝑘 0 ì 0 ( Ĥ12 ) 𝑘,𝛾 ( Ĥ12 ) 𝛾,𝛾 ( ì H 12 ) 𝜃,𝛾 0 ( Ĥ22 ) 𝛾,𝛾 ( ì H 22 ) 𝜃,𝛾 ( ì H
⊤ 12 ) 𝑘,𝜃 ( ì H ⊤ 12 ) 𝛾,𝜃 (H 12 ) 𝜃,𝜃 ì 0 ⊤ ( ì H ⊤ 22 ) 𝛾,𝜃 (H 22 ) 𝜃,𝜃 • ì 𝑘 1 ì 𝛾 1 𝜃 1 ì 𝑘 2 ì 𝛾 2 𝜃 2 = - ì 0 ì 0 0 ì 0 ì 0 0 (3.55)
57 which shows that they remain constant during the cooling, while the second sub-system of equations

2h ( Ĥ11 ) 𝑥,𝑥 0 ì 0 ( Ĥ21 ) 𝑥,𝑥 ( Ĥ21 ) 𝜒,𝑥 ( ì H 21 ) 𝑐,𝑥 0 ( Ĥ11 ) 𝜒,𝜒 ì 0 ( Ĥ21 ) 𝑥,𝜒 ( Ĥ21 ) 𝜒,𝜒 ( ì H 21 ) 𝑐,𝜒 ì 0 ⊤ ì 0 ⊤ (H 11 ) 𝑐,𝑐 ( ì H ⊤ 21 ) 𝑥,𝑐 ( ì H ⊤ 21 ) 𝜒,𝑐 (H 21 ) 𝑐,𝑐 ( Ĥ12 ) 𝑥,𝑥 ( Ĥ12 ) 𝜒,𝑥 ( ì H 12 ) 𝑐,𝑥 ( Ĥ22 ) 𝑥,𝑥 0 ì 0 ( Ĥ12 ) 𝑥,𝜒 ( Ĥ12 ) 𝜒,𝜒 ( ì H 12 ) 𝑐,𝜒 0 ( Ĥ22 ) 𝜒,𝜒 ì 0 ( ì H ⊤ 12 ) 𝑥,𝑐 ( ì H ⊤ 12 ) 𝜒,𝑐 (H 12 ) 𝑐,𝑐 ì 0 ⊤ ì 0 ⊤ (H 22 ) 𝑐,𝑐 • ì 𝑥 1 ì 𝜒 1 𝑐 1 ì 𝑥 2 ì 𝜒 2 𝑐 2 = - ì ∇ 𝑥 1 E tot ì ∇ 𝜒 1 E tot 𝜕 𝑐 1 E tot ì ∇ 𝑥 2 E tot ì ∇ 𝜒 2 E tot 𝜕 𝑐 2 E tot (3.56)
reveal that the parameters are not only driven by the negative gradient of the total energy with respect to that parameter anymore but there are additional effects related to the other Gaussian. The reason is that the overlap between two Gaussians plays an important role in the evolution. If for example the overlap ⟨𝑔 2 |𝑔 1 ⟩ approaches zero then the matrices Ĥ21 and Ĥ12 would vanish since they are proportional to the overlap. Thus, the eq. 3.56 would diagonalize and we would have the situation of two independent Gaussians evolving in the mean-field potential. However, in the opposite case when the Gaussians are exactly the same, that is Ĥ21 = Ĥ12 = Ĥ11 = Ĥ22 , the matrix in eq. 3.56 is not invertible anymore and thus singular. Hence, the equation of motion avoid this singularity which is a quantum effect of the model. Specifically, Pauli-blocking is ensured since the energy levels cannot drop to lower levels in that way. The cooling process is eventually stopped when the total energy converges similar to eq. 2.120. Finally, all the properties of the Gaussians are saved into an initialization data file that are read for the second part of the numerical framework, which handles the time evolution of a heavy-ion collision.

Implementation of the dynamics

In second part of the code, we initialize two nuclei, which have been prepared from the implementation of the statics, on a bigger grid to simulate a heavy-ion collision. To be specific, the lattice is extended to range from -80 fm to 80 fm in the beam direction and -40 fm to 40 fm in both directions in the orthogonal plane. We will discuss briefly the initial condition of the collision kinematics, as well as the invariant transformation to boost the nuclei.

Initialization of the kinematics

The simulation of the collision is performed in the center-of-mass system. Furthermore, we need to input several parameters such as incident energy per nucleon 𝐸 inc , impact parameter 𝑏 and the relative distance 𝑅 between the two nuclei. Since the impact parameter provides the distance in the orthogonal plane, we determine together with 𝑅 the displacement of the two nuclei in configuration space. The total shift is given by

𝑍 1 = 𝑍 𝑅 𝐴 2 𝐴 1 + 𝐴 2 , 𝑍 2 = -𝑍 𝑅 𝐴 1 𝐴 1 + 𝐴 2 (3.57) 𝑋 1 = 𝑋 𝑅 𝐴 2 𝐴 1 + 𝐴 2 , 𝑋 2 = -𝑋 𝑅 𝐴 1 𝐴 1 + 𝐴 2 (3.58) with 𝑍 𝑅 = 𝑅 cos 𝜃 , 𝑋 𝑅 = 𝑅 sin 𝜃 (3.59)
where 𝐴 1 , 𝐴 2 are the atomic mass numbers and 𝑍 1 , 𝑍 2 the charge numbers of the respective nuclei 1 and 2. The rotation angle 𝜃 depends on the impact parameter b. In the routine there is an additional feature which takes into account the effects of Coulomb. The reason is that two nuclei are supposed to be initialized from infinity but we start at a given position. However, we note that in all our calculations the nuclei are initialized with a sufficiently large distance, which is why the effect of the Coulomb trajectory can be calculated dynamically during the approaching phase. Different to models where the mean-field is not as accurately described due to the localization of the wave packet the consideration of the Coulomb effects are more important. In that case the intrinsic nature of the semiclassical evolution equation drives the system to a Boltzmann distribution which is why the nuclei have to be initialized at close distance. Than it is important to account for the momentum shift due to the Coulomb trajectory. For our calculations we initialize at far enough distances, typically 20 fm, with minimal corrections of the Coulomb effects according to that routine. Similarly, the displacements of the Fermi spheres in momentum space are obtained from the incident energy per nucleon. The shifts in momentum space read

𝐾 𝑧,1 = 𝐾 𝑧,𝑅 𝐴 2 𝐴 1 + 𝐴 2 , 𝐾 𝑧,2 = -𝐾 𝑧,𝑅 𝐴 1 𝐴 1 + 𝐴 2 (3.60) 𝐾 𝑥,1 = 𝐾 𝑥,𝑅 𝐴 2 𝐴 1 + 𝐴 2 , 𝐾 𝑥,2 = -𝐾 𝑧,𝑅 𝐴 1 𝐴 1 + 𝐴 2 (3.61)
where ì 𝐾 𝑅 is the momentum shift, which are also given from the routine that handles the treatment of the Coulomb distance.

Invariant translation and boost transformation

The next step is to apply the translation to the position and momentum parameters of our wave functions. A transformation is defined by the preservation the scalar product. Since in the definition of our Gaussian basis function 𝑔 𝑚 , see eq. 2.14, the position only appears together with the position parameter as a difference (𝑥 -𝑥 𝑚 ) we conclude that

T𝑋 𝑥 𝑚 = 𝑥 𝑚 + 𝑋 -→ T𝑋 𝑔 𝑚 (𝑥) = 𝑔 𝑚 (𝑥 -𝑋 ) . (3.62)
where T𝑋 is the translation operator acting on the position parameter 𝑥 𝑚 by shifting it by 𝑋 . As a result, the spatial translation keeps the scalar product between Gaussians of the same nucleus invariant

T ì 𝑋 𝑔 𝑛 T ì 𝑋 𝑔 𝑚 = ∭ 𝑔 𝑛 (ì 𝑥 -ì 𝑋 )𝑔 𝑚 (ì 𝑥 -ì 𝑋 )d 3 𝑥 (3.63) = ∭ 𝑔 𝑛 (ì 𝑥 ′ )𝑔 𝑚 (ì 𝑥 ′ )d 3 𝑥 ′ = ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ . (3.64)
For a boost transformation T𝐾

T𝐾 𝑘 𝑚 = 𝑘 𝑚 + 𝐾 -→ T𝐾 𝑔 𝑚 (𝑥) = 𝑔 𝑚 (𝑥)𝑒 𝑖 𝐾(𝑥 -𝑥 𝑚 ) (3.65)
we notice that the Gaussian 𝑔 𝑚 picks up an addition phase shift instead of being translated. The constant phase shift 𝑒 -𝑖 𝐾𝑥 𝑚 can be absorbed into the phase parameter 𝜃 𝑚 by T𝜃 𝜃 𝑚 = 𝜃 𝑚 + 𝐾𝑥 𝑚 -→ T𝜃 T𝐾 𝑔 𝑚 (𝑥) = 𝑔 𝑚 (𝑥)𝑒 𝑖 𝐾𝑥 .

(3.66) Thus, we have constructed a combination of transformations for the momentum and the phase parameter which preserves the scalar product as well

T𝜃 T ì 𝐾 𝑔 𝑛 T𝜃 T ì 𝐾 𝑔 𝑚 = ∭ 𝑔 𝑛 (ì 𝑥)𝑒 -𝑖 ì 𝐾 ì 𝑥 𝑔 𝑚 (ì 𝑥)𝑒 𝑖 ì 𝐾 ì 𝑥 d 3 𝑥 = ⟨𝑔 𝑛 |𝑔 𝑚 ⟩ . (3.67)
Therefore we summarize that the Gaussian parameters of the first nucleus are transformed according to

ì 𝑥 𝑚 -→ ì 𝑥 𝑚 + ì 𝑋 1 (3.68) ì 𝑘 𝑚 -→ ì 𝑘 𝑚 + ì 𝐾 1 (3.69) 𝜃 𝑚 -→ 𝜃 𝑚 + ì 𝐾 1 • ì 𝑥 𝑚 (3.70)
where ì 𝑋 1 is the spatial displacement and ì 𝐾 1 the momentum displacement provided from the previous subsection. The same transformation rules hold for parameters of the second nucleus exchanging the displacements to ì 𝑋 2 and ì 𝐾 2 .

Time-evolution routine

For the evolution of the parameters we use a second-order predictor-corrector integration method.

Specifically, we first apply a two-step Adams-Bashforth method [Bashforth1883, Butcher2008], which takes into account the gradient of the previous time step

p(𝑛+1) = 𝑝 (𝑛-1) + ∆𝑡 2 (3 𝑝 (𝑛) -𝑝 (𝑛-1) ) (3.71)
where 𝑝 (𝑛) ∈ P is a parameter of a Gaussian at the n-th time iteration and p(𝑛+1) is the corresponding intermediate predictor parameter. In a second step, we calculate the gradients of the predictor to obtain the final parameter 𝑝 (𝑛+1) through a simple Euler method

𝑝 (𝑛+1) = 𝑝 (𝑛-1) + ∆𝑡 2 ( 𝑝 (𝑛) + p(𝑛) ) . (3.72)
In order to loop this integration method, we need the gradients of two consecutive time steps. Hence we initialize the time-evolution routine by a simplified version of the original method. For the first step we do not use a Adams-Bashforth method but just a forward time evolution for the first predictor and calculate the corrector as usual p(1) = 𝑝 (0) + ∆𝑡 𝑝 0 (3.73)

𝑝 (1) = 𝑝 (0) + ∆𝑡 2 ( 𝑝 (0) + p(1) ) (3.74)
where we save the gradients 𝑝 0 . With the new set of parameters 𝑝 (1) we calculate the new gradients 𝑝 (1) and start the time-evolution according to eq. 3.71 and eq. 3.72. For each iteration, we keep the gradients 𝑝 (𝑛) for the subsequent iteration. To get an idea on how the Gaussians behave we calculate explicitly again the evolution matrix of a single-Gaussian and two-Gaussian nucleonic wave-function. However, this time we investigate the time-evolution matrix in eq. 2.11 and not the imaginary time-evolution matrix of eq. 2.122.

Time-evolution of a single-Gaussian wave function

For the time-evolution of the single Gaussian of the nucleonic wave function, we require the Hesse matrix elements again. Luckily, we already calculated the ingredients for the matrix elements in eq. 3.16 to eq. 3.30. Different to the evolution matrix of the cooling, the evolution matrix of the time-evolution is given by the imaginary contribution of the Hesse matrix elements. Thus, we finally obtain

ĤEoM = 0 0 Ĥ𝑥,𝑘 0 ì 0 ì 0 0 0 0 Ĥ𝜒,𝛾 ì H 𝑐,𝛾 ì 0 Ĥ𝑘,𝑥 0 0 0 ì H 𝑐,𝑥 ì 0 0 Ĥ𝛾,𝜒 0 0 ì H 𝑐,𝜒 ì 0 ì 0 ⊤ ì H ⊤ 𝛾,𝑐 ì H ⊤ 𝑥,𝑐 ì H ⊤ 𝜒,𝑐 0 H 𝜃,𝑐 ì 0 ⊤ ì 0 ⊤ ì 0 ⊤ ì 0 ⊤ H 𝑐,𝜃 0
(3.75)
where ĤEoM = ℑ ì ∇ 𝑞 𝑔 1 ì ∇ 𝑝 𝑔 1 is the equation-of-motion evolution matrix, and its Ĥ𝑞,𝑝 the submatrices of the Hesse matrix where 𝑞, 𝑝 ∈ P 1 are a parameters of the Gaussian. The imaginary part of the sub-matrices have the simple expressions

Ĥ𝑥,𝑘 ℑ = 1 2 0 0 0 1 2 0 0 0 1 2 ⟨𝑔 1 |𝑔 1 ⟩ (3.76) Ĥ𝛾,𝜒 ℑ = 1 2 0 0 0 1 2 0 0 0 1 2 ⟨𝑔 1 |𝑔 1 ⟩ (3.77) ì H 𝑐,𝛾 ℑ = 1 𝑐 1 ì 𝜒 1 ⟨𝑔 1 |𝑔 1 ⟩ (3.78) ì H 𝑐,𝑥 ℑ = - 1 𝑐 1 ì 𝑘 1 ⟨𝑔 1 |𝑔 1 ⟩ (3.79) ì H 𝑐,𝜒 ℑ = - 1 𝑐 1 ì 𝛾 1 ì 𝜒 1 ⟨𝑔 1 |𝑔 1 ⟩ (3.80) H 𝜃,𝑐 ℑ = 1 𝑐 1 ⟨𝑔 1 |𝑔 1 ⟩ (3.81) (3.82) Ĥ𝑝,𝑞 = - Ĥ ⊤ 𝑞,𝑝 (3.83) ì H 𝑝,𝑐 = -ì H ⊤ 𝑐,𝛾 (3.84)
H 𝑐,𝜃 = -H 𝜃,𝑐 (3.85)
where ì 𝛾 1 ì 𝜒 1 corresponds to a point-wise division and the last three lines describe the property of the skew-symmetry of the evolution matrix. The total equation of motion thus reads

2h 0 0 Ĥ𝑥,𝑘 0 ì 0 ì 0 0 0 0 Ĥ𝛾,𝜒 ì H 𝑐,𝛾 ì 0 Ĥ𝑘,𝑥 0 0 0 ì H 𝑐,𝑥 ì 0 0 Ĥ𝛾,𝜒 0 0 ì H 𝑐,𝜒 ì 0 ì 0 ⊤ ì H ⊤ 𝛾,𝑐 ì H ⊤ 𝑥,𝑐 ì H ⊤ 𝜒,𝑐 0 H 𝜃,𝑐 ì 0 ⊤ ì 0 ⊤ ì 0 ⊤ ì 0 ⊤ H 𝑐,𝜃 0 • ì 𝑘 1 ì 𝛾 1 ì 𝑥 1 ì 𝜒 1 𝑐 1 𝜃 1 = h 𝑚 ì 𝑘 ⟨𝑔 1 |𝑔 1 ⟩ 4h 𝑚 ì 𝜒 1 ⊗ ì 𝛾 1 ⟨𝑔 1 |𝑔 1 ⟩ ì ∇ 𝑥 1 E pot ì ∇ 𝜒 1 E kin + ì ∇ 𝜒 1 E pot 𝜕 𝑐 1 E kin + 𝜕 𝑐 1 E pot 0 (3.86)
where ì 𝜒 1 ⊗ ì 𝛾 1 is the point-wise product of those two vectors. This contribution comes from the derivative of the kinetic energy with respect to 𝛾. We recall that the kinetic energy is expressed as

E kin = h 2𝑚 ( ì 𝑘 2 1 + 𝜙 1,𝑥 + 𝜙 1,𝑦 + 𝜙 1,𝑧 )𝑐 2 1 (3.87) with 𝜙 1,𝑖 = 𝜒 1,𝑖 ( 1 4𝜒 2 1,𝑖 + 4𝛾 2 1,𝑖 ) . (3.88)
Let us look at the trivial equations which are already insightful. First, we recover for a single-Gaussian wave function the classical-like evolution in a mean-field

ì 𝑥 1 = ì 𝑘 𝑚 (3.89) ì 𝑘 1 + 2𝑐 1 𝑐 1 h ì 𝛾 1 ì 𝜒 1 = - 1 𝑐 2 1 h ì ∇ 𝑥 1 E pot (3.90)
where the only difference is that the change of momentum is additionally influenced by the change of the coefficient. However, from the last line of the equation of motion we obtain

𝑐 1 = 0 (3.91)
which means that the normalization of the single Gaussian does not change during time. This shows that in the equation of motion, the particle number is indeed inherently implemented. As a result, the second term on the left-hand side of eq. 3.90 reduces to the classical evolution of an independent particle. Furthermore, we can list the evolution equation of the variance and 𝛾 as well

ì 𝜒 1 = 4 h ì 𝜒 1 ⊗ ì 𝛾 1 (3.92) ì 𝛾 1 = h 8𝑚 ì 𝜒 2 1 - 2h 𝑚 ì 𝛾 2 1 - 1 𝑐 2 1 h ì ∇ 𝜒 1 E pot (3.93)
where ì 𝜒 2 1 and ì 𝛾 2 1 are the point-wise products. These are the classical evolution of variances influenced by a mean-field potential. That implies, that a localized wave function described by a single Gaussian, behaves like a free particle in an effective mean-field. This matches the interpretation, that Gaussian wave function are coherent states that follow classical trajectories. In the crude approximation that every nucleonic wave function localizes, we transform the model into a semi-classical molecular dynamics model. Re-introducing an anti-symmetrization of the many body system would recover the fermionic nature, leading us to models such as FMD or AMD. In the precedent model DYWAN, the decoherence assumption decoupled each Gaussian, bringing the model to an evolution of independent test particles within a mean-field similar to Vlasov models.

Time-evolution of a two-Gaussian wave function

In our approach we have also de-localized nucleonic wave function. That means they are decomposed into a superposition of localized gaussian-type functions. If we look at the first case where we extend the nucleonic wave function to have two Gaussian, we can learn how they affect each other. For that reason we take a look at the evolution matrix of a single-Gaussian evolution of a Gaussian 𝑔 1 and 𝑔 2 . That would mean the Gaussians of a wave function would evolve equivalently to an independent particle. However on the off-diagonal positions of the matrix we have submatrices

( Ĥ21 ) 𝑞,𝑝 = ì ∇ 𝑞 𝑔 2 ì ∇ 𝑝 𝑔 1 (3.95) ( Ĥ21 ) 𝑞,𝑝 = -( Ĥ † 12 ) 𝑝,𝑞 (3.96)
which are related to the evolution of the product gaussians 𝑔 21 or 𝑔 12 respectively. Since, the product gaussian depends on each parameter of 𝑔 1 and 𝑔 2 the evolution matrices Ĥ𝑖 𝑗 do not usually have vanishing components. the sub-matrices for the cooling in eq. 3.54 were simplified since we initialized the momentum, 𝛾 and phase parameter with zero. Hence, we were able to decouple some equations.

During the real-time evolution that is not the case. The parameters can take any value. As a result these off-diagonal matrix elements will always couple parameters of one gaussian to parameters of the other. That leads to the effect that the additional term with ì H 𝑐,𝑝 in eq. 3.86 does not drop as it was in the independent case. Rather, these terms are now activated to include the change of the coefficients of the gaussians. In the single-Gaussian case, the coefficient 𝑐 1 was not able to vary since there was no counterpart to absorb the change to keep the normalization of the wave function. In contrast, if there are multiple Gaussians available, the coefficients are able to shift, since the increase of one can be balanced by the decrease of the others. Effectively, an equation such as eq. 3.90 describes the change of the momentum due to the change of the corresponding parameter but also the change due to the flow of weight. As we discussed already for the cooling, the contribution of the off-diagonal Hesse sub-matrices are proportional to the overlap, that means, if the Gaussians are far apart the effects are negligible but when two Gaussians approach the same state, the evolution matrix would become singular in that situation. Therefore the Gaussians are not able to collapse which create an quantum effect that is related to the Pauli-principle. Finally, we understand that a set of Gaussians evolve as the independent evolution in a mean-field potential if they are sufficiently far apart. If the overlap increases significantly the off-diagonal matrix elements add a modification to the evolution equation which takes into account the interferences between each other.

Summary

We have built a mean-field approach which has implemented the quantum effects directly into the evolution matrix by adding off-diagonal elements. These elements correspond to interferences between Gaussians of the same nucleonic wave function. They restrict each wave function to evolve in their respective sub-manifold to ensure orthogonality during the evolution. The model can be cast into a semi-classical evolution by restricting the Gaussians to be sufficiently distant in phase-space. Then the overlap is negligible and each Gaussian would evolve independently in the effective mean-field. This is the case for a decoherence approximation. The approximation fails as soon as two Gaussians start to overlap significantly. Therefore, we do not apply a decoherence but calculate at each timestep the full evolution matrix with all the off-diagonal contributions. By inverting that matrix, we obtain the gradients of the parameters. The gradients are used to simulate the time-evolution of a heavy-ion collision. The integration method follows a predictor-corrector method where the predictor is calculated from a two-step Adams-Bashforth method, that combines the derivatives of the previous two iteration and a corrector which uses the information of the predictor to calculate the parameters of the next time-step. The nuclei are initialized by the kinematic conditions where we put in the incident energy per nucleon, the impact parameter as well as the relative distance between the nuclei. Then an invariant transformation of the wave functions have been performed to ensure the orthogonality during the boost transformation. It was crucial to introduce a phase shift transformation to preserve the nuclear properties that we have set up in the initialization of the nuclei. During that part, we created a ground state by performing first a self-consistency routine based in an harmonic oscillator fit of the effective potential. After the fit converged, we decomposed the harmonic oscillator states into our dynamical gaussian basis functions and used a bottom-up sorting method to construct the three-dimensional nucleus. In the end we performed a cooling process based on an imaginary-time evolution with the same matrix elements as we later need in the real-time evolution. We showed that the cooling process just displaces the configuration of the Gaussians, keeping the momentum, 𝛾 and phase parameter at zero. The other parameters are then following the negative gradients of the total energy to find deeper total energy. However, we had to constrain the coefficients 𝑐 𝑖 in order to maintain particle number conservation and give a lower bound to the cooling evolution.

Chapter 4

Validation of the numerical scheme

In this chapter we will show the results of the dynamical processes, namely the cooling and the time evolution procedure, to validate the numerical stability of the model. In specific, we are interested in the system of 40/48 Ca since there are, as of the time as this thesis is written, a large interest in the isospin diffusion which is investigated at GANIL by the INDRA collaboration [Fable2023]. However, we also want to check if heavier systems are properly represented, hence we are testing as well 124 Sn which was previously also experimentally studied for isospin diffusion with 84 Kr [Piantelli2015]. glance, it might suggest that the density decreases in total, which would mean that the system loses matter. However, since we perform explicitly a renormalization during the cooling procedure, mass conservation is directly conserved. To understand where the missing density escapes we can take a look at fig. 4.3. In these two plots, we show the three-dimensional spatial distribution of 124 Sn where we cut out one half of the distribution to see the cross section. The plots are shown in a color gradient, where red is at over saturation above 𝜌 0 and towards blue it is at low density. On the left, you can see the cross section before the cooling and on the right the distribution after the cooling process. We observe that initially the nucleus has a cross-like internal structure. Further, the central density is indicated to have an over-saturated density, which we have previously confirmed. After the cooling, on the right, the density is now homogeneously distributed in the center at around saturation density. We deduce that the high density has been redistributed evenly in the radial direction to achieve a more favorable ground state. That means, we did not lose matter in the process but rather the internal structure which was enforced due to the harmonic oscillator initialization and lead to the cross-like shape that was over-saturated, has been broken up and reorganized in a way such that there are no internal forces anymore which would make the system excited and unstable. Additionally, to obtain the momentum distribution one may use the Fourier-transform of the gaussians as displayed in eq. 2.17 and sum them up. Effective neutron potential of the approximated harmonic oscillator states before the cooling(dashed blue) and effective neutron potential after the cooling procedure (solid green) for 48 Ca. Right: Effective neutron potential of the approximated harmonic oscillator states before the cooling(dashed blue) and effective neutron potential after the cooling procedure (solid green) for 124 Sn.

Test of the binding energy

The next crucial test to check is the agreement with the binding energy. We note that the potential significantly reduced the repulsion in the center after cooling as well in fig. 4.4. In this figure we display the effective neutron potential before the cooling and after the cooling, again in dashed blue and solid green respectively. We notice, that not only the central part of the potential flattens which reduces the internal forces within the nucleus but also that the potential deepens after the cooling. This is different to the self-consistency routine in subsection 2.2.1 where the artificial flattening resulted to a lift of the potential well. In the case of cooling, we deepen the potential which leads to a higher mean binding energy of the nucleus.

We show the time evolution of 124 Sn in fig. 4.5. As we can see, the cooling process does not affect the mean binding energy significantly. It only adjusts the initial mean binding energy within 0.1 AMeV.

Further, we remark that the evolution is not strictly descending but it occasionally increases slightly.

That is because we added to the cooling process a renormalization of the particle number. If we would follow the pure cooling evolution equation then the mean binding energy would strictly ascend, but due to the restriction to the particle number conservation, we limit the parameter space, which in turn leads to the evolution that we see in fig. 4.5. The evolution can be smoothened with smaller time steps but for reasons of computation time we chose a time step of ∆𝑡 = 0.1 fm/𝑐. Thus, the important role continues afterwards with the real-time evolution equation. In both cases the time evolution after the cooling process as a stable evolution up to 𝑡 = 100 fm/c. Afterwards, the mean binding energy starts to deviate, however the error is still well within 0.1 AMeV. It shows that the model is suited to perform over long times for a dynamical evolution. In addition, we included in the diagrams the experimental data, in blue, for the binding energies of both nuclei. [Audi2003] As we can see, the binding energy for 48 Ca is much closer, within the range of 0.1 AMeV, than it is the case for 124 Sn which is higher than the experimental value by over 0.2 AMeV. This discrepancy can be traced back to the problem, that first the nucleus was initialized from the harmonic oscillator. Since the magic numbers differ at higher levels we observe that for the heavier 124 Sn the binding energy differs much more than for the lighter 48 Ca compared to the experimental value. Second, we also approximated the effective nuclear potential without the spin-orbit effect, which is why we are not able to recreate the actual magic numbers in the imaginary as well as the real time-evolution. Nonetheless, we accept in the current version of the model that the ground states are only approximated, in order to be able to use them in a dynamical of the binding energies for all the nuclei. In this chart the binding energies are calculated directly from the harmonic oscillator with the SkM parameterization of the Skyrme force. As we have discussed before, the cooling does not change the binding energy significantly, but also cannot transition the nuclei to the correct magic numbers. Therefore, we can see that the binding energies do not overlap directly with the experimental data, but they are fluctuating around those data especially at higher mass numbers A. The focus for now, is to initialize a nucleus which is stable over time without strong internal forces, which the cooling achieves. One task in the future might be to include the spin-orbit contribution to tackle exactly this problem of the statics.

Test of the scalar product

One of the crucial tests for this model is to check the scalar product during the time-evolution among nucleonic wave functions, since the first the evolution equations are based on the fact that the de- localized wave functions are orthogonal. This was the advantage compared to AMD models where the fermionic nature is explicitly accounted for by anti-symmetrizing the A-body wave functions which increases the computational time significantly. In our model, we initialized the wave function from orthogonal states and evolve them in such a way, that they keep the orthogonality over time. Therefore, we show the time evolution in fig. 4.9. In this plot, at each time step, each point represents the scalar product between two nucleonic wave functions in the 48 Ca nucleus. As we can see, already in the initialization the wave functions are not completely orthogonal. The reason for that is that we approximated the harmonic wave functions by a sum of finitely many gaussians. Since the fit is not exact, this can lead to some off-set of the orthogonality. The important point is that the majority of the scalar products are less than 0.1 and all of them are below 0.2. Finally, the more important information of this diagram is that even though the initial configuration is not perfectly orthogonal, the equation of motions obey exactly the property to conserve the scalar product over time as we have derived it in eq. 2.124. where shell effects are not accounted for [Urban2012]. Hence, our quantum mechanical model with de-localized wave functions is also able to confirm that the resonance is induced in the neutron-rich 48 Ca as it is the case for other models [Baran2012].

Hence, we validate with the tests in this chapter, that the nuclei are initialized in such a way that they are have the correct density distribution and agree approximately with the experimental binding energies. More crucial is that they keep the initial scalar products constant. That means, the quantum effects are accounted for properly in the evolution equations and therefore the nucleons keep the Fermi distribution contrary to semi-classical models which drive the system to Boltzmann distributions. As a result we see in the dynamics, that the density distribution stays compact and does not smear out over time as in the case of Vlasov models.

Chapter 5

Conclusion

The concept of the model historically began with the idea of wavelets in the precedent model called DYWAN by Francois Sebille. It was a first attempt to define a dynamical function basis to solve the many-body Schrödinger equation. However, the idea was not fully executed but suffered from the problem that the decoherence assumption was performed in the derivation of the equation of motion.

As a result, the original version just complicated a semi-classical evolution equation. In fact, the self-consistent harmonic oscillator routine is inspired from that model. The difference is that in the current method the harmonic oscillator states are accurately constructed and fitted. Previously, the Gaussians were placed according to approximated distributions of the Wigner-transformations of the harmonic oscillator states. Already at that point any quantum effects were washed out by the scheme.

In addition to that, there was no cooling procedure which means that the initial density distribution was approximated by a harmonic oscillator density distribution. We showed that this construction lead to the problem of high central density and thus strong repulsive potentials in the center. Lastly, the main problem of the model was that the equations of motion were heavily simplified such that each Gaussian evolved just as an independent test particle in a mean-field potential. Therefore, we concluded that DYWAN unfortunately did not provide a solid foundation to explore beyond-mean-field effects, because already the mean-field evolution was strongly simplified in contrast to the complicated construction of the initialization. As a matter of fact, the original goal of the project was to explore large-amplitude fluctuations in a quantum mean-field model. However, we took the decision to rework the whole model starting from the equation of motion, to have a consistent quantum-mechanical framework which inherits the original idea of dynamical basis functions.

Thus the main purpose of this thesis shifted to building up a new framework which keeps the detail of a quantum mean-field TDHF description, but establishes additionally a foundation to include beyond mean-field extensions of large-amplitude fluctuations similar to semi-classical models. Therefore, we derived a method which is analogous to the test-particle approach in semi-classical models.

The difference is that in semi-classical models the test-particles describe a phase-space portion and are independent of each other. However, in our approach the test-particles are not fully independent of each other, but rather a set of localized gaussian wave functions keep the information that they belong to the same nucleonic degree of freedom. Based on this idea we derived rigorously the full detail of the equation of motion for a parameterized basis function. It turns out that a Gaussian would behave just as an independent test-particle in a semi-classical model if the overlap with other Gaussians belonging to the same nucleonic wave function goes to zero. This corresponds to the condition of the decoherence approximation which was used in a similar model called DYWAN. The problem of that model was that situations where the overlap is not negligible, the simple evolution equations did not account for the fact that the nucleonic wave function is constrained by Pauli-principle. In our case, we succeeded to introduce the dynamical basis functions properly in a quantum mean-field scheme without going through a Wigner-transformation. We showed that it connects the two leading quantum approaches for nuclear reactions, TDHF-models and Molecular Dynamics models for fermions. In theory, the model can be categorized as a generalization of those Molecular Dynamics models, since they treat the Schrödinger equation in the limit of one localized Gaussian per nucleon, while we increase the function space to superpositions of localized Gaussians. In the limit of infinite Gaussians one would recover the TDHF scheme. Molecular Dynamics models clearly have the advantage that they treat as little parameters as possible, which improves the calculation time. However, it comes with the expense of reducing the description of the mean-field effects which are important at low to Fermi energies. That is the reason why, we increased the number of Gaussians depending on the energy level of the nucleonic wave function. We find for the equation of motion of those Gaussians that, in addition to the semi-classical evolution, there are terms related to the interference between the Gaussians. Due to these terms the cooling manages to reduce the central density and later on keep a stable evolution of the single nucleus.

This was not the case in DYWAN, which started also from an harmonic oscillator decomposition. As we noticed, the harmonic oscillator naturally lead to a very high central density, which consequently brought the strong repulsion of the effective nuclear potential. Thus, since DYWAN is the decoherence approximation of our model, it was not able to escape the problem that was introduced from the harmonic oscillator. Another, difference to many other dynamical models is that the coefficients of the dynamical basis are time-dependent and complex. Originally, the problem arose because in DYWAN the weights of Gaussians were initialized with different values due to the initialization. However, they were time-independent in that model. As a result high weighted Gaussians dominated the evolution of the system. Therefore, in the rework of the approach we included the amplitude parameterization of the dynamical basis. In order to include the evolution of the amplitude in the equation of motion we had to introduce another parameter, so that the skew-symmetric Hesse matrix is invertible. Naturally, since the Gaussians are defined to represent wave functions, the amplitude is not only real valued but complex. Hence we decided to introduce a phase parameter as well to cover the full function space. It turns out that the inclusion of the time-dependence of complex coefficients leads to two additional equations in the equation of motion that help the system keep crucial conservation laws, namely particle number and energy conservation.

All in all, in this thesis we achieved to provide a consistent model with quantum features which is able to be modified in a way that residual terms similar to those in BL-approaches, based on the nucleon-nucleon correlations, can be implemented to create a beyond-mean-field framework which directly propagates the nucleonic wave functions and treats their interferences and resulting quantum fluctuations in the mean-field evolution. At variance from a quantum propagation of coherent states, semi-classical models lose and blur the information of nucleonic degrees of freedom in their decoherent evolution by construction. Therefore, in implementations where nucleon-nucleon collisions and fluctuations are derived from nucleonic degrees of freedom, such as BLOB, a new set of nucleonic semi-classical wave packets has to be redefined in term of distribution functions at each time step.

Energy and momentum conservation should be satisfied, as well as Pauli-blocking factors in phase space, without any propagation from a previous coherent set. The advantage in our framework is that nucleonic degrees of freedom, once implemented in the construction of the initial system, persist and are propagated throughout the evolution. Their de-localization is not the effect of a mere dissolution in a hydrodynamic (Boltzmann) description, but it is the result of tracking quantum interferences among coherent states in the equation of motion along the potential landscape. Furthermore, the collision can be checked on the overlap rather than the occupancy in phase-space. The open question still remains to define a proper treatment of an analogous Boltzmann collision term for this scheme. In the limit of molecular dynamics where we represent one Gaussian per nucleon, the treatment simplifies since the wave function localizes. In our case, a possibility is to represent one Gaussian as a possible localized Dans le chapitre 3, nous nous concentrons sur la mise en oeuvre numérique de l'approche du champ moyen qui a été discutée dans le chapitre 2. Nous implémentons les effets quantiques directement dans la matrice d'évolution en ajoutant des éléments hors-diagonaux. Ces éléments correspondent aux interférences entre les gaussiennes de la même fonction d'onde nucléonique. Ils limitent chaque fonction d'onde à évoluer dans doubtsa sous-variété respective afin de garantir l'orthogonalité au cours de base gaussiennes dynamiques et utilisons une méthode de tri ascendante pour construire le noyau en trois diménsions. Enfin, nous réalisons un processus de refroidissement basé sur une évolution en temps imaginaire avec les mêmes éléments de matrice que dans l'évolution en temps réel. Nous montrons que le processus de refroidissement ne fait que déplacer la configuration des gaussiennes, en maintenant les paramètres de quantité de mouvement, de 𝛾 et de phase à zéro. Les autres paramètres suivent alors les gradients négatifs de l'énergie totale pour trouver une énergie totale plus profonde. Cependant, nous contraignons les coefficients 𝑐 𝑖 afin de maintenir la conservation du nombre de particules et de donner une limite inférieure à l'évolution du refroidissement.

Enfin, dans le chapitre 4, nous validons le modèle par des tests sur différents noyaux, à savoir les isotopes Ca et Sn. Tout d'abord, nous étudions la densité de saturation des noyaux et nous constatons que la configuration initiale de l'oscillateur harmonique sans refroidissement conduit à une sursaturation au centre du noyau et donc à de fortes forces internes provenant du potentiel répulsif. En effectuant la procédure de refroidissement supplémentaire, la densité centrale élevée se disperse dans tout le noyau. La densité de saturation en résulte donc réduite ainsi que le potentiel répulsif, ce qui conduit à une configuration beaucoup plus stable. Ensuite, les tests de l'énergie de liaison montrent que l'énergie de liaison des isotopes ne change pas de manière significative lors du refroidissement.

Ainsi, le rôle important du refroidissement est le réarrangement des gaussiennes pour obtenir une configuration plus proche de la densité de saturation et moins de forces internes pour éviter les instabilités initiales. En outre, nous notons que les énergies de liaison ne sont qu'approximatives puisque nous utilisons un'interaction de Skyrme simplifiée sans les contributions de masse effective et de spin-orbite, et une paramétrisation appelée SV-bas [Klupfel2009] qui n'a pas été réajustée à notre interaction de Skyrme. Néanmoins, les noyaux sont suffisamment bien décrits et stables tout au long de l'évolution.

Nous montrons en effet que le produit scalaire reste constant au cours du temps malgré le fait que les produits scalaires initiaux ne sont pas complètement orthogonaux pour des raisons numériques dans l'initialisation. Finalement, dans une première application du modèle, nous étudions qualitativement une collision périphérique d'isotopes du calcium et observons un effet typique de champ moyen, c'està-dire la résonance pygmée. trouver un schéma pour le terme de collision qui déplace le nombre existant de gaussiennes ou redéfinit complètement l'état final avec un nouveau nombre de gaussiennes. L'avantage est que nous ne sommes théoriquement pas limités au nombre de gaussiennes comme dans la dynamique moléculaire. Il existe généralement une grande liberté pour choisir la manière dont le terme de collision agit. Certaines des conditions qu'il doit remplir sont la conservation de l'énergie et de la quantité de mouvement ainsi que le principe de Pauli. Toutefois, dans l'état actuel du modèle, le principe de Pauli peut être facilement vérifié par le produit scalaire dans l'état final, ce qui signifie que l'état final doit être construit de manière à préserver l'énergie totale et les impulsions. Pendant la période de cette thèse, nous avons réussi à définir. les fondements du modèle dynamique quantique. Nous avons pu essayer des termes de collision simples qui font tourner les fonctions d'onde nucléoniques dans l'espace des impulsions, mais cette méthode simple n'a pas été couronnée de succès. Cependant, la base du modèle semble prometteuse pour ouvrir de nouvelles voies de recherche future, afin d'inclure des termes de collision plus sophistiqués qui permettent d'explorer des phénomènes allant des basses énergies jusqu'aux énergies de Fermi.

87

Titre:

  Introduction d'une base dynamique dans un schéma de champ moyen de fonctions d'ondes orthogonales pour les réactions nucléaires Mots clés: Dynamique nucléaire, Champ moyen quantique, TDHF, Modèles de transport. Résumé: La dynamique des collisions d'ions lourds révèle une variété de mécanismes différents qui sont attribués à la combinaison d'effets collectifs et dissipatifs. Dans les collisions d'ions lourds, un large éventail de phénomènes est exploré en fonction de diverses conditions d'énergie incidente, d'asymétrie d'isospin, de paramètre d'impact et d'autres propriétés. La physique sous-jacente des différents mécanismes peut être comprise qualitativement à partir du régime énergétique auquel le processus se produit. Aux basses énergies jusqu'à environ 15 MeV par nucléon, les collisions nucléon-nucléon à deux corps sont supprimées en raison du blocage de Pauli dans les états finaux. Ainsi, la dynamique est bien décrite par les effets collectifs à longue portée du potentiel de champ moyen. Cependant, dès que le régime énergétique atteint plusieurs dizaines ou centaines de MeV, aux énergies dites de Fermi à intermédiaires, les interactions nucléoniques à deux corps doivent être incluses en plus du comportement collectif. Enfin, aux énergies élevées à relativistes, que l'on peut qualifier de régime participantspectateur, l'interaction à deux corps de courte portée éclipse la contribution du champ moyen et dicte l'évolution principale. Pour chacune de ces plages d'énergie, il existe différents modèles spécialisés qui sont bien adaptés à cette situation spécifique. Par exemple, les mécanismes de faible énergie sont efficacement décrits dans le cadre de la méthode Hartree-Fock dépendant du temps (TDHF), tandis que les fluctuations de grande amplitude aux énergies de Fermi sont mieux décrites dans des approches stochastiques de champ moyen (équation de Boltzmann-Langevin) ou dans des modèles de dynamique moléculaire. C'est un défi de longue date pour les théories nucléaires à N corps de traiter correctement les types d'instabilités concurrentes au seuil entre les énergies de Fermi et les basses énergies dans une image unifiée. À l'état actuel des connaissances, les modèles aux énergies de Fermi ont tendance à simplifier, dans une certaine mesure, la phénoménologie du champ moyen, des courants d'isospin aux modes collectifs, et, à l'inverse, les modèles s'approchant des basses énergies manquent de fluctuations de grande amplitude et de mécanismes où des transformations profondes, voire des bifurcations, se produisent dans le système. Une approche théorique capable de couvrir la transition entre les énergies de Fermi et les basses énergies dans une seule description complète, serait un cadre bien adapté pour étudier l'évolution des fragments et des clusters nucléaires en fonction du temps et de la densité. Nous adaptons une formulation TDHF stochastique en partant de l'équation de Schroedinger à N corps, puisqu'il s'agit de la physique fondamentale qui détermine le comportement collectif de tout système quantique à ces énergies. L'idée principale de ce modèle est la décomposition des fonctions d'onde nucléoniques non locales dans une base de fonctions dynamiques. Dans conséquence, les propriétés du champ moyen sont préservées puisque la fonction d'onde nucléonique n'est pas contrainte d'être localisée et, de plus, le système est préparé dans un ensemble de fonctions de base mobiles afin de suivre un schéma similaire à celui de l'approche semi-classique analogue de Boltzmann-Langevin. Title: Introducing a dynamical basis in a mean-field scheme of orthogonal wave functions for nuclear reactions Keywords: Nuclear dynamics, quantum mean-field, TDHF,transport models Abstract: The dynamics of heavy-ion collisions unveils a diversity of various mechanisms which are related to the combination of collective and dissipative effects. The wide range of phenomena depends on numerous conditions such as incident energy, isospin asymmetry, impact parameter and other properties. The underlying physics for the different mechanisms can be understood qualitatively from the relevant energy regime at which the process occurs. At the low energies up to around 15 MeV per nucleon, two-body nucleon-nucleon collisions are suppressed due to Pauli blocking in the final states. Thus, the dynamics is well described by the long-range collective effects of the mean-field potential. However, as soon as the energy regime reaches up to several tens to hundreds of MeV, at the so-called Fermi to intermediate energies, the two-body nucleonic interactions have to be included in addition to the collective behaviour.

Fig. 1 .

 1 Fig. 1.1 indicates schematically how to distinguish between different incident-energy regimes in heavyion collisions. In this representation, low energies may be identified from the Coulomb barrier to around 15 MeV per nucleon and correspond to projectile wavelengths exceeding the size of the target nucleus.

Figure 1

 1 Figure 1.1: A schematic diagram of the De Broglie wavelength 𝜆 as a function of the incident energy per nucleon in heavy-ion collisions.At low energies, up to 15 AMeV, the probing wavelength is larger than the target nucleus which means that mean-field effects dominate over two-body nucleon-nucleon collisions. In green it is indicated that collisionless models such as the TDHF models describe very well this energy regime. At high energies, around 1 AGeV, the probing wavelength is of the size of the nucleonic interaction which means that collisional effects dominate the dynamical process, which is indicated in red as the participant-spectator regime. Lastly in blue, around Fermi energy, the probing wavelength is of the size of the target nucleus which shows that both mean field as well as two-body correlations have to be taken into account, which means that beyond-mean-field models have to be used to describe this regime accurately. See text.

Figure 1 . 2 :

 12 Figure 1.2: Schematic illustration of the effects of fluctuations in HIC. Left: The dynamical trajectories of identically prepared systems behave deterministically without stochastic fluctuations. Right:The trajectories can spread out due to fluctuations during the evolution and end up in various configurations. This figure is adapted from[Napolitani2013bis]

  .1) The fluctuation term 𝛿𝐼 [𝑓 ] is assumed to be normally distributed with a vanishing mean and a variance to ⟨𝐼 coll [𝑓 ]⟩. As illustrated on the left of fig. 1.2, kinetic equations (Boltzmann or dissipative extensions of the TDHF) predict a single evolution trajectory for HIC due to their deterministic nature (one-body leading contribution) with a variance resulting from dissipative contributions. The inclusion of 𝛿𝐼 [𝑓 ] induces fluctuations in the phase-space density 𝑓 during the evolution and drives the system towards a wide variety of exit channels, as is illustrated on the right side of fig. 1.2. In early approaches the fluctuation seed 𝛿𝐼 [𝑓 ] was introduced as a projection on a suited space [Suraud1992, Belkacem1993, Guarnera1996, Colonna1998] but, in this case, the fluctuation amplitude could result to be underestimated. A more complete approach (BLOB [Napolitani2013, Napolitani2017]) introduces 𝛿𝐼 [𝑓 ] in full phase space intermittently at successive time steps, letting fluctuations develop spontaneously and continuously during the whole dynamical process. In particular, the method ensures that the occupancy variance at equilibrium equals 𝑓 (1 -𝑓 ) in a phase-space cell ℎ 3 and corresponds to the movement of extended portions of phase space which have the size of a nucleon and represent 'semi-classical' wave-packets. In other terms, the residual contribution ⟨𝐼 coll [𝑓 ]⟩ + 𝛿𝐼 [𝑓 ] should carry nucleon-nucleon

Figure 2 . 1 :

 21 Figure 2.1: Gaussian spatial distribution |𝑔| 2 (𝑥) with central position 𝑥 0 and variance 𝜒.

  which have been constrained by the largest possible number of experimental observables to reproduce as close as possible the mean-field results that a fundamental interaction would provide [Chabanat1997,Brown2013, Chappert2015]. The two most used effective interactions are the Gogny [Decharge1980,Hashimoto2012, Chappert2008] and Skyrme[Skyrme1956, Danielewicz2009, Dutra2012, Danielewicz2014]

Figure 2 . 3 :

 23 Figure 2.3: Graphs of a Woods-Saxon type potential for different masses (40,84,124).

Figure 2

 2 Figure 2.4: A schematic overview of self-consistency iteration routine.

Figure 2 . 5 :

 25 Figure 2.5: Decomposition of a second excited harmonic oscillator state (solid black) into a superposition of three Gaussian base function (dashed blue).

Figure 3 . 1 :

 31 Figure3.1: A schematic flow chart which shows the connection between different nuclear transport models. The common feature of all models is the many-body Schrödinger equation which determines the microscopic evolution. The general solution is given by a superposition of Slater determinants. However, practical models start from the single-Slater approximation. Then, there are two paths depending on whether the single-particle states are orthogonal or not. Our model decomposes Hartree-Fock states into dynamical Gaussians. We can draw the connection to F/AMD if we go to the limit of one Gaussian per nucleon. Finally BUU-models are connected through the Wigner transformation.

Figure 3 . 2 :

 32 Figure3.2: Procedure of the harmonic oscillator fit in the self-consistency routine. Left, the initial neutron potential described by the Woods-Saxon potential (solid black) for a 48 Ca nucleus. The smooth approximation of the potential is shown in dash-dotted green and the corresponding harmonic oscillator fit in dashed blue. Right, the same procedure after the first iteration with an effective Skyrme potential.

Figure 3 . 3 :

 33 Figure3.3: Procedure of the harmonic oscillator fit in the self-consistency routine. Left, the initial proton potential described by the Woods-Saxon and the Coulomb potential (solid black) for a 48 Ca nucleus. The smooth approximation of the potential is shown in dash-dotted green and the corresponding harmonic oscillator fit in dashed blue. Right, the same procedure after the first iteration with an effective Skyrme and Coulomb potential.

  3.3 the same plots are shown for the proton potential of the 48 Ca where we have to include additionally the repulsive Coulomb potential. That is the reason why the effective proton potentials converge towards zero at long distances from above whereas the effective neutron potentials in fig.3.2 converge from below. The proton potential is repulsive at long distances because of the Coulomb potential while the neutron potential is only attractive because of the effective nuclear interaction. The addition of the Coulomb potential slightly modifies

Figure 3 . 5 :

 35 Figure3.5: Procedure of the harmonic oscillator fit in the self-consistency routine. Left, the initial proton potential described by the Woods-Saxon and the Coulomb potential (solid black) for a 124 Sn nucleus. The smooth approximation of the potential is shown in dash-dotted green and the corresponding harmonic oscillator fit in dashed blue. Right, the same procedure after the first iteration with an effective Skyrme and Coulomb potential.

Figure 3 . 6 :Figure 3 . 7 :

 3637 Figure 3.6: Evolution of the effective potentials of a 48 Ca nucleus during the self-consistency routine.Left, the progression of the neutron effective potential starting from the initial potential (solid black) through intermediate potentials (dashed black) ending at the final potential (solid blue). Right, the same progression for the proton effective potential.

Figure 3 . 8 :

 38 Figure3.8: Decomposition of the harmonic wave function 𝜑 𝑛 (𝑥) as defined in eq. 2.106. The index 𝑛 denotes the level of the eigenstates which are displayed in solid black in each subplot above. The dashed green lines are the decomposition of the wave function into the sum of gaussian-type functions. The number of Gaussians is determined by the number of peaks of the corresponding wavefunction. For clearity we left out the decomposition for 𝑛 = 4 and 𝑛 = 5. The positions, amplitudes and variances of the Gaussians have been determined using a Gauss-Newton minimization algorithm.

  same structure as in the cooling case in eq. 3.51. On the diagonal positions we find the Hesse matrices ĤEoM 11 and ĤEoM 22

4. 1

 1 fig.4.1. Depicted are the total density distributions after the self-consistency routine of the harmonic oscillator systems for 48 Ca and 124 Sn on the left and right, respectively. The solid black graphs show the density distribution calculated directly from the definition of eq. 2.106, while the dashed blue graphs show the approximate density distributions after the decomposition of the harmonic oscillator wave functions into a superposition of Gaussians. First we notice that the fit is very close to the original density distribution, however, for both nuclei the central density is significantly higher than the saturation density 𝜌 0 = 0.16 fm -3 . That is the reason, why we observed a strong repulsive center in the effective nuclear potentials. It is already the case for lighter nuclei such as 48 Ca which is why

Figure 4 . 1 :

 41 Figure 4.1: Left: Density distribution of the exact harmonic oscillator states (solid black) and the approximation by the decomposition in Gaussians (dashed blue) for 48 Ca. Right: Density distribution of the exact harmonic oscillator states (solid black) and the approximation by the decomposition in Gaussians (dashed blue) for 124 Sn.

Figure 4 . 2 :

 42 Figure4.2: Left: Density distribution of the approximated harmonic oscillator states (dashed blue) and density distribution after the cooling procedure (solid green) for 48 Ca. Right: Density distribution of the approximated harmonic oscillator states (dashed blue) and density distribution after the cooling procedure (solid green) for 124 Sn.

Figure 4

 4 Figure 4.3: 3D cross section plots of the spatial density distribution of 124 Sn with a color graded density scale.

Figure 4

 4 Figure4.4: Left: Effective neutron potential of the approximated harmonic oscillator states before the cooling(dashed blue) and effective neutron potential after the cooling procedure (solid green) for 48 Ca. Right: Effective neutron potential of the approximated harmonic oscillator states before the cooling(dashed blue) and effective neutron potential after the cooling procedure (solid green) for 124 Sn.

Figure 4 . 5 :

 45 Figure 4.5: Time evolution of the mean binding energy of 124 Sn during the cooling procedure. The time step for this calculation was ∆𝑡 = 0.1fm/𝑐

Figure 4 . 6 :

 46 Figure 4.6: Time evolution of the mean binding energy of 48 Ca during the cooling procedure and the real time evolution.

Figure 4 . 7 :Figure 4 . 8 :

 4748 Figure 4.7: Time evolution of the mean binding energy of 124 Sn during the cooling procedure and the real time evolution.

Figure 4 . 9 :

 49 Figure 4.9: Time-evolution of the scalar product of each possible pair of nucleonic wave functions in the 48 Ca nucleus.

Figure 4 . 10 :

 410 Figure 4.10: Snapshots at different times (𝑡 = 0, 40, 80, 110, 150, 200 and 250 fm/c) during a peripheral collision of the isospin asymmetric 40 Ca + 48 Ca collision at 35 AMeV with an impact parameter 𝑏 = 9 fm. The red color indicates the proton distribution and the blue color indicates the neutron distribution.

  state of the wave function and let it collide with another Gaussian of a different nucleon. In this way, we could naturally introduce stochastic fluctuation to the collision term since a nucleonic wave function could take any of the different Gaussian as its representation. The problem to solve is the creation of the final states. The work for the future is to find a scheme for the collision term which either moves the existing number of Gaussians or completely redefine the final state with a new number of Gaussians. The advantage is that we are theoretically not restricted to the number of Gaussians like in molecular dynamics. There is generally a large freedom to choose how the collision term acts. Some of scopique totalement nouvelle de la dynamique nucléaire, basée sur les corrélations particule-particule et les observables liées à la chronologie. Ceci grâce à la granularité plus élevée des nouveaux détecteurs, à une meilleure résolution cinématique et à un pouvoir d'identification amélioré. Par consequent, il sera possible d'extraire l'équation d'état de la matière nucléaire des collisions d'ions lourds en profitant d'observables plus microscopiques, contrairement aux études précédentes, où les observables thermodynamiques devaient être privilégiées. De même, les observables de corrélation à haute résolution permettront de mieux comprendre le processus de formation des fragments et des clusters plutôt que d'étudier les propriétés de l'état final. Dans cette perspective, les théories de la fonctionnelle de la densité et les modèles microscopiques seront essentiels à la fois pour préparer de nouvelles expériences et pour imposer des contraintes encore plus strictes dans la connexion entre la théorie et les observables. Dans le cadre de cette thèse, l'accent est mis sur la description théorique du processus dynamique des collisions d'ions lourds, en particulier aux basses énergies, aux énergies de Fermi et aux énergies intermédiaires, de quelques MeV à des centaines de MeV par nucléon, dans le but de construire un nouveau cadre de modélisation microscopique bien adapté. L'intérêt d'étudier les réactions nucléaires dans cette gamme d'énergie vient du fait qu'elles sont caractérisées par une grande variété de phénomènes, allant des modes collectifs aux processus dissipatifs, en raison de la combinaison de plusieurs contributions comme la structure nucléaire, le champ moyen, les corrélations nucléon-nucléon, les fluctuations et plusieurs types d'instabilités. À cet égard, se baser sur la théorie Hartree-Fock dépendante du temps (TDHF) comme point de départ est une voie privilégiée. L'approche TDHF résout l'équation de Schrödinger dépendante du temps pour la fonction d'onde à N corps dans l'espace des coordonnées. Ce modèle traite les nucléons comme des particules quantiques délocalisées et inclut le potentiel de champ moyen. Il peut fournir une description fiable des effets quantiques et des effets de champ moyen. Cependant, la théorie TDHF présente encore certaines limites, notamment parce que les extensions pour introduire des fluctuations de grande amplitude, au-delà du champ moyen, sont plus compliquées que dans un modèle de Boltzmann semi-classique. C'est pourquoi elle est généralement restreinte à la structure nucléaire et aux réactions à basse énergie. L'objectif de cette thèse est de développer les bases d'un cadre théorique qui couvre les caractéristiques essentielles de la mécanique quantique, mais qui ouvre également une voie pour inclure aussi des fluctuations de grande amplitude. Par conséquent, nous présentons un modèle de transport basé sur la théorie de champ moyen quantique TDHF qui, au même temps, peut être aussi traitée comme un modèle de dynamique moléculaire généralisé qui propage des fonctions d'onde fermioniques délocalisées pour une meilleure description du champ moyen. Il s'agit de l'analogue quantique des modèles semi-classiques dans lesquels des particules tests sont utilisées pour décrire la distribution de la densité. Cependant, contrairement aux modèles semi-classiques, les particules tests sont substitués par des composantes de fonction d'onde qui ne sont plus indépendantes les unes des autres. Puisque certains composants appartiennent à la même fonction d'onde nucléonique, ils gardent l'information sur les degrés de liberté nucléoniques. Par conséquent, les équations d'évolution du champ moyen sont modifiées de telle sorte que les effets quantiques sont directement inclus dans l'équation du mouvement, ce qui nous permet d'explorer également les basses énergies. L'objectif est de donner une nouvelle base qui traite explicitement le comportement collectif des phénomènes de basse énergie mais qui peut être étendue par des effets au-delà du champ moyen grâce à la connaissance des degrés de liberté nucléoniques. Toutefois, l'accent sera mis sur la réalisation d'une telle approche théorique à partir de l'équation du mouvement sous-jacente. L'extension vers l'inclusion des collisions nucléon-nucléon reste toutefois une étape à compléter pour l'avenir. Au chapitre 2, nous expliquons en détail comment nous dérivons les fondements théoriques de notre modèle. doubtNous approximons la fonction d'onde à N corps par un seul déterminant de Slater et ajoutons une condition d'orthogonalité à l'ensemble des fonctions d'onde doubtindividuelles, comme c'est souvent le cas dans les modèles TDHF. Cela nous permet d'abandonner les termes doubtmixtes de l'antisymétrisation dans l'opérateur densité (voir eq. 2.5). Par conséquent, le lagrangien se sépare en une somme de lagrangiens qui ne décrit que le comportement d'une seule fonction d'onde nucléonique. L'idée centrale qui distingue notre approche des autres est que les fonctions d'onde nucléoniques sont décomposées en une superposition de gaussiennes 𝑔 𝑚 caractérisées chacune par une position ì 𝑥 𝑚 , un nombre d'ondes angulaire ì 𝑘 𝑚 , une variance spatiale ì 𝜒 𝑚 , une incertitude de l'espace de phase ì 𝛾 𝑚 et un paramètre de poids complexe 𝑐 𝑚 𝑒 6𝑢𝜃 𝑚 . Pour cette nouvelle représentation, nous dérivons explicitement les éléments de la matrice jacobienne et de la matrice hessienne, qui apparaissent dans les coefficients de la matrice d'évolution. En outre, nous introduisons un potentiel effectif simplifié de type Skyrme et le potentiel de Coulomb pour calculer les valeurs moyennes de l'énergie. Ensuite, nous montrons comment nous prévoyons théoriquement d'initialiser le noyau dans notre système. Nous partons d'une estimation initiale du potentiel nucléaire, le potentiel de Woods-Saxon, et lançons une routine auto-cohérente, qui approxime le potentiel nucléaire avec un oscillateur harmonique afin d'obtenir une famille d'états propres d'énergie orthogonaux. Ces états propres forment une distribution de densité qui définit à son tour les potentiels effectifs de Skyrme et de Coulomb mentionnés précédemment. La routine itère sur cette boucle, où le nouveau potentiel est à nouveau ajusté, jusqu'à ce que l'énergie totale de l'oscillateur harmonique converge. Enfin, nous décomposons les fonctions d'onde harmoniques délocalisées en gaussiennes, où le nombre de gaussiennes dépend du niveau d'excitation de l'oscillateur harmonique, et nous appliquons à la fonction de base gaussienne une procédure de refroidissement, basée sur l'évolution en temps imaginaire. Ainsi, nous devrions être en mesure de trouver une solution statique qui est proche du véritable minimum du potentiel effectif. A la fin du chapitre, nous montrons également que la condition d'orthogonalité, le nombre de masse et l'énergie totale du système sont théoriquement conservés tout au long de l'évolution et à tout moment. Plus précisément, l'introduction de l'amplitude et de la phase dépendant du temps permet de traiter directement l'équation de conservation.

  Dans l'ensemble, dans cette thèse, nous avons réussi à développer un modèle cohérent dans un formalisme quantique qui peut être étendu vers l'inclusion de termes résiduels basées sur les corrélations nucléon-nucléon, dans le même esprit des approches Boltzmann Langevin. De telles extensions pourront être mises en oeuvre pour créer un cadre au-delà du champ moyen qui propage directement les fonctions d'onde nucléoniques et traite leurs interférences et les fluctuations quantiques qui en résul-tent dans l'évolution du champ moyen. À la différence d'une propagation quantique d'états cohérents, les modèles semi-classiques perdent et brouillent l'information des degrés de liberté nucléoniques dans leur évolution décohérente par construction. En particulier, dans les implémentations où les collisions nucléon-nucléon et les fluctuations exploitent les degrés de liberté nucléoniques, comme dans l'approche Boltzmann-Langevin One-Body (BLOB), un nouvel ensemble de paquets d'ondes semi-classiques nucléoniques doit être redéfini en termes de fonctions de distribution à chaque pas de temps. Dans ce cas, la conservation de l'énergie et de la quantité de mouvement doit être satisfaite, ainsi que les facteurs de blocage de Pauli dans l'espace des phases, sans aucune propagation à partir d'un ensemble cohérent précédent. À l'inverse, L'avantage de notre cadre est que les degrés de liberté nucléoniques, une fois mis en oeuvre dans la construction du système initial, persistent et sont propagés tout au long de l'évolution. Leur délocalisation n'est pas l'effet d'une simple dissolution dans une description hydrodynamique (Boltzmann), mais le résultat du suivi des interférences quantiques entre les états cohérents dans l'équation du mouvement le long du potentiel. En outre, dans la prospective de rajouter un terme de collision, le taux de collisions nucléon-nucléon peut être déterminé directement en fonction du recouvrement des fonctions d'onde plutôt que contraint par l'occupation dans l'espace des phases. La question ouverte reste comment définir un traitement approprié d'un terme de collision analogue aux implémentations Boltzmann Langevin dans ce schéma. Dans la limite de la dynamique moléculaire où nous représentons une gaussienne par nucléon, le traitement se simplifie puisque la fonction d'onde est localisée. Dans notre cas, il sera envisageable de représenter une composante gaussienne comme un possible état localisé de la fonction d'onde et de la laisser entrer en collision avec une autre composante gaussienne d'un nucléon différent. De cette manière, nous pourrions naturellement introduire une fluctuation stochastique dans le terme de collision puisqu'une fonction d'onde nucléonique pourrait être approximé par n'importe quelle composante gaussienne qui contribue à la représenter. Le problème à résoudre est la définition des états finaux. Plus précisément, le travail à venir consiste à

  𝜑 𝑖 𝑈 Skyrme 𝜕 𝑞 𝑚 𝜑 𝑖 (2.82) where 𝜌 𝑖 is either the neutron or proton density depending on which one the wave function 𝜑 𝑖 belongs to. The last eq. 2.82 show that we can rewrite the equation either into a derivative of the expectation value of the effective potential or the real part of the matrix element of the effective potential with the wave function 𝜑 𝑖 and its tangent vector 𝜕 𝑞 𝑚 𝜑 𝑖 . Correspondingly, we get the same expression for the derivative of the Coulomb energy 𝜕E Coul 𝜕𝑞 𝑚 = ℜ 𝜑 𝑖 𝑈 Coul 𝜕 𝑞 𝑚 𝜑 𝑖 . (2.83) Altogether, we can combine both potentials to work with a single total effective potential 𝑈 pot = 𝑈 Skyrme + 𝑈 Coul (2.84) to treat the derivatives of the total potential energy E pot = E Skyrme + E Coul . With view to evaluate those expectation values, we prepare a set of necessary integrals. Specifically, we want to know the moments of the potential expectation value of two gaussians 𝑔 𝑛 and 𝑔 𝑚 up to the second order

	𝜕E Skyrme	=	∭	𝛿H Skyrme	𝜕𝜌 𝑖	d 3 𝑥	(2.80)
	𝜕𝑞 𝑚			𝛿𝜌 𝑖		𝜕𝑞 𝑚
		=	∭	𝑈 Skyrme	𝜕𝑞 𝑚 𝜕(𝜑 𝑖 𝜑 𝑖 )	d 3 𝑥	(2.81)
		= ℜ

Using the chain rule we find for the derivative of 𝑞 𝑚 ∈ P 𝑚 of a wave function 𝜑 𝑖

  . If we rewrite the evolution equation eq. 2.10 in terms of the imaginary time 𝜏 = 𝑖𝑡

	2hℑ	𝜕𝜑 𝑖	𝜕𝑔 𝑚	= 2hℑ 𝑖	𝜕𝜑 𝑖	𝜕𝑔 𝑚
		𝜕𝑡	𝜕𝑞 𝜇			𝜕𝜏	𝜕𝑞 𝜇
				= -2hℜ	𝜕𝜑 𝑖	𝜕𝑔 𝑚	=	𝜕E tot	(2.121)
						𝜕𝜏	𝜕𝑞 𝜇	𝜕𝑞 𝜇
	=⇒ 2hℜ	𝜕𝜑 𝑖	𝜕𝑔 𝑚	= -	𝜕E tot	(2.122)
		𝜕𝜏	𝜕𝑞 𝜇		𝜕𝑞 𝜇

  Procedure of the harmonic oscillator fit in the self-consistency routine. Left, the initial neutron potential described by the Woods-Saxon potential (solid black) for a 124 Sn nucleus. The smooth approximation of the potential is shown in dash-dotted green and the corresponding harmonic oscillator fit in dashed blue. Right, the same procedure after the first iteration with an effective Skyrme potential.

			Initial Iteration					Second Iteration
					Woods-Saxon					Skyrme
					Smooth approx.					Smooth approx.
					HO fit					HO fit
						-8	-4	0	4	8
								x [fm]
	-50 -40 -30 -20 -10 0 Figure 3.4: -8 Potential [MeV] 10	-4	0 Initial Iteration	4	8 Woods-Saxon + Coul. Smooth approx. HO fit	-8	-4	0 Second Iteration	4	8 Skyrme + Coul Smooth approx. HO fit
			x [fm]					x [fm]

  𝑁 max 𝑔 𝑁 max Table3.1: Given a maximum level 𝑁 max we list the degeneracy of the level 𝑔 𝑁 max , the sum 𝑁 max

			𝑁 max 𝑁	𝑔 𝑁 2 𝑁 max 𝑁	𝑔 𝑁
	0	1	1	2
	1	3	4	8
	2	6	10	20
	3	10	20	40
	4	15	35	70
	5	21	56	112

l'évolution. Le modèle peut être transformé en une évolution semi-classique en limitant les gaussiennes à une distance suffisante dans l'espace des phases. Le récouvrement est alors négligeable et chaque gaussienne évolue indépendamment dans le champ moyen effectif. C'est le cas de l'approximation de décohérence. L'approximation échoue dès que deux gaussiennes commencent à se récouvrir de manière significative. Par conséquent, nous n'appliquons pas de décohérence, mais nous calculons à chaque pas de temps la matrice d'évolution complète avec toutes les contributions hors diagonale. En inversant cette matrice, nous obtenons les gradients des paramètres. Ces gradients sont utilisés pour simuler l'évolution temporelle d'une collision d'ions lourds. La méthode d'intégration suit une méthode prédicteur-correcteur où le prédicteur est calculé à partir d'une méthode Adams-Bashforth en deux étapes, qui combine les dérivées des deux itérations précédentes et un correcteur qui utilise l'information du prédicteur pour calculer les paramètres du pas de temps suivant. Le système des noyaux cible et projectile est défini dans le centre de masse en assignant aux deux noyaux les conditions cinématiques correspondantes à l'énergie incidente par nucléon, au paramètre d'impact et à la distance relative initiale. Ensuite, une transformation invariante des fonctions d'onde est effectuée pour assurer l'orthogonalité pendant la transformation d'impulsion. Il est crucial d'introduire une transformation de déphasage pour préserver les propriétés nucléaires que nous avons établies lors de l'initialisation des noyaux. Dans l'initialisation, nous définissons un état fondamental en effectuant d'abord une routine d'autoconsistance basée sur un ajustement de l'oscillateur harmonique du potentiel effectif. Une fois que l'ajustement a convergé, nous décomposons les états de l'oscillateur harmonique en fonctions de

Acknowledgements

But as the model is currently set up, Pauli principle can be easily checked by the scalar product in the final state, which means that the final state has to be constructed in a way that we preserve total energy and momenta. During the period of this thesis, we achieved to rework the foundation of the dynamical quantum model. We were able to try out simple collision terms which rotate the nucleonic wave functions in momentum space, but this simple method was not successful. However, the basis of the model looks promising to open up new avenues for future research to include more sophisticated collision terms which allow to explore phenomena ranging from low to Fermi energies.

Résumé étendu en français

Les collisions d'ions lourds représentent un laboratoire terrestre pour sonder l'équation d'état nucléaire [Baldo2012, Baran2005]. L'équation d'état est un objet si fondamental qu'elle ne caractérise pas seulement les noyaux finis, mais aussi une grande variété de scénarios astrophysiques où la matière dense est impliquée [Steiner2005, Oertel2017, Burgio2021, Fantina2022]. Les collisions d'ions lourds Nexplorent l'évolution en densité d'un système nucléaire [Ring1980], des modes collectifs jusqu'aux perturbations désordonnées. En particulier, un large éventail de situations se présentent, où les neutrons et les protons oscillent en phase ou hors phase. Dans les deux cas, cela conduit souvent à une transformation profonde de la structure du système. Étant donnée la propriété de cohésion de l'interaction nucléaire au sein du noyau et de la matière nucléaire, la voie de sortie pourrait également être un processus catastrophique, où l'ensemble du système se désintègre en clusters et fragments [Chomaz2004].

Contrairement aux conditions astrophysiques, cette corrélation entre les collisions d'ions lourds et l'équation d'état est très indirecte pour deux raisons principales. Tout d'abord, surtout au début du processus, les échelles de temps sont comparables au temps de relaxation de l'interaction nucléaire, de sorte qu'un comportement chaotique et des instabilités de volume peuvent apparaître. Deuxièmement, lorsque l'on passe des conditions de la matière nucléaire à des systèmes finis, la surface du système entre également en jeu, imposant des géométries non triviales et des instabilités de surface. Pour traiter une phénoménologie aussi complexe dans un cadre théorique, il faut combiner les caractéristiques à un corps (champ moyen) et au-delà du champ moyenqui, jusqu'à présent, ont été largement approximées dans les modèles de transport ou n'ont pas été corrélées. Il s'agit de la description du caractère non local de la fonction d'onde et des fluctuations quantiques découlant de leur interférence (un corps), ou des fluctuations de grande amplitude déclenchées par des effets dissipatifs (au-delà du champ moyen). Ces aspects du processus nucléaire sont traduits en observables expérimentales qui pourront être étudiés dans les mesures à venir, telles que les corrélations particule-particule et la chronologie d'émission des particules.

Dans cette optique, ces expériences sur les collisions d'ions lourds devraient apporter une vision micro-