TIME-DEPENDENT DESCRIPTION OF HEAVY-ION COLLISIONS

By

Hao Lin

A DISSERTATION

Submitted to

Michigan State University

in partial fulﬁllment of the requirements

for the degree of

Physics — Doctor of Philosophy

2020

ABSTRACT

TIME-DEPENDENT DESCRIPTION OF HEAVY-ION COLLISIONS

By

Hao Lin

In this thesis, we aim to advance the time-dependent transport theories for the description

of heavy-ion collisions, from two perspectives. As an attempt to address multifragmentation

in nuclear collisions, we develop a stochastic transport model based on one-body Langevin

dynamics. The new model is subsequently tested and benchmarked with a series of other

existing models with satisfaction. The model is also applied to address and conﬁrm the

so-called “hierarchy eﬀect” observed in the multifragmentation for certain systems around

Fermi energies. Parallel to the development towards a stochastic theory, we also extend an

approach based on non-equilibrium Green’s function for the description of correlated nuclear

systems in one dimension.

Firstly, we present a new framework to treat the dissipation and ﬂuctuation dynamics

associated with nucleon-nucleon scattering in heavy-ion collisions. The two-body collisions

are eﬀectively described in terms of the diﬀusion of nucleons in viscous nuclear media, gov-

erned by a set of Langevin equations in momentum space. The new framework combined

with the usual mean-ﬁeld dynamics, forming the basis of the new stochastic model, can be

used to simulate heavy-ion collisions at intermediate energies.

Subsequently, as a proof of principle for the new model, we simulate Au + Au reactions

at 100 MeV/nucleon and at 400 MeV/nucleon and look at observables such as rapidity

distribution and ﬂow as a function of rapidity. The results are found to be consistent with

other existing models under the same constrained conditions. To demonstrate the model’s

ability to describe multifragmentation, we also study the formation of fragments in Sn +

Sn reactions at 50 MeV/nucleon, and the fragment distribution and properties are discussed

and compared to two other models commonly employed for collisions.

Next, we move on to tackle the “hierarchy eﬀect” observed experimentally for reactions

around Fermi energies. We simulate Ta + Au at 39.6 MeV/nucleon and compare mainly

the charge and velocity distributions of the fragments from the QP with experimental data.

Our simulation results can reproduce the trends observed in data, and a semi-quantitative

agreement can be reached. This is the ﬁrst time, to our knowledge, that one has succeeded

in addressing the “hierarchy eﬀect” with a dynamical model. The simulation of U + C is

also discussed.

Finally, we present a fully quantum-mechanical model based on non-equilibrium Green’s

function, with short-range two-body correlations incorporated as an extension. We examine

its applications to one-dimensional nuclear systems, such as the preparation and properties

of the ground states, the isovector oscillation of symmetric systems and the boosting of a

“slab” in a periodic box. In particular, the dissipation brought by two-body correlations and

the Galilean covariance of the theory are demonstrated. These studies lay the groundwork

for the future exploration of collisions of correlated nuclear systems in one dimension.

ACKNOWLEDGMENTS

I would ﬁrst like to thank my advisor, Pawel Danielewicz, for not only his academic guidance

but also his support, patience and liberty. He is that kind of liberal and laid-back person, who

would not complain about my not following his words and meandering with little progress. I

can be pessimistic and, at times, depressed, and his kindness and optimism never fail to ease

my nerves and lighten me up a bit. I also appreciate that he has provided me with a great

many support and opportunities to travel and to attend conferences, and that he would go

the extra mile to hunt for job opportunities for me and connect me to people.

I am grateful to my other committee members, Betty Tsang, Carlo Piermarocchi, Filom-

ena Nunes, Gary Westfall, and Morten Hjorth-Jensen, for their guidance, be it critical ques-

tions or heartfelt advice, and for, of course, putting up with my talks every year.

I would also like to thank Scott Pratt and Maria Colonna for fruitful discussions and

suggestion and for their time and eﬀorts in writing references for me. I also want to say

thank you to Hermann Wolter for providing me with a comfortable and hospitable stay in

Munich as well as oﬀering references.

I also want to thank Elizabeth Deliyski and Gillian Olsen for their diligent work and

assistance along the way, especially when I had created mess before and after my trips.

Among my favorite lecturers at MSU are Scott Bogner, Morten Hjorth-Jensen and Vladimir

Zelenvinsky.

The spiritual support and accompany of my peers means a great deal to me.

I shall

take this opportunity to send my gratitude and best wishes to my cohort including, but not

limited to, Maxwell Cao, Mengzhi Chen, Xian-Gai Deng, Linda H’lophe, Jun Hong, Tong

Li, Kuan-Yu Lin, Dan Liu, Jian Liu, Amy Lovell, Hossein Mahzoon, Xingze Mao, Zachary

iv

Matheson, Pierre Nzabahimana, Avik Sarkar, Tommy Tsang, Jonathan Wong, Heda Zhang,

Boyao Zhu, Kuan Zhu ... and, last but not least, whoever is reading this!

The work in this thesis is supported by the U.S. National Science Foundation under

Grants PHY-1403906 and PHY-1520971, the U.S. Department of Energy under Grant DE-

SC0019209, and the Dissertation Completion Fellowship from Michigan State University.

v

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Heavy-Ion Collisions
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Transport Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Fluctuations and Quantum eﬀects . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

1
1
4
7
10

Initialization with the Thomas-Fermi equations

Chapter 2 One-body Langevin transport model . . . . . . . . . . . . . . . .
11
11
2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
2.2 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
2.2.1 The Boltzmann framework . . . . . . . . . . . . . . . . . . . . . . . .
14
2.2.2 The mean-ﬁeld dynamics . . . . . . . . . . . . . . . . . . . . . . . . .
15
2.2.3 The dissipation and ﬂuctuation dynamics . . . . . . . . . . . . . . . .
19
2.2.4
. . . . . . . . . . . .
21
2.2.5 Brownian motions of nucleon wave-packets . . . . . . . . . . . . . . .
2.2.5.1 Partition of test particles into nucleon wave-packets . . . . .
21
2.2.5.2 Evaluations of the Langevin equation’s coeﬃcients on a lattice 23
2.2.5.3 Momentum transfer from the nuclear medium to the nucleon
wave-packet . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.5.4 Recoil for conservation of momentum and energy . . . . . .
2.2.5.5 Pauli-blocking procedure . . . . . . . . . . . . . . . . . . . .
Summary of implementation of Brownian motions . . . . . .
2.2.5.6

24
25
26
26

Chapter 3 Applications to heavy-ion collisions at intermediate energies .
3.1 Benchmarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 One-body observables for Au + Au system . . . . . . . . . . . . . . .
3.1.1.1 Rapidity distribution dN/dy . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
3.1.1.2 Average in-plane ﬂow (cid:104)px/A(cid:105)
3.1.2 Fragmentation dynamics with Sn + Sn at 50 MeV/nucleon . . . . . .
Summary and discussion . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.3

Chapter 4

Investigating “hierarchy eﬀect” in multifragmentation at Fermi
. . . . . . . . . . . .
energies with the Brownian motion model
4.1 Simulation setup: physical inputs and selection of impact parameter b . . . .
4.2 Comparison with experimental data . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Charge distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Velocity distribution . . . . . . . . . . . . . . . . . . . . . . . . . . .

30
30
30
31
33
37
42

45
46
48
54
58

vi

4.3 Absence of the “hierarchy eﬀect” for U + C . . . . . . . . . . . . . . . . . . .
4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 5 Nonequilibrium Green’s function approach for many-body quan-

tum systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Extending the non-equilibrium Green’s function approach for nuclear systems
5.3 Preparing ﬁnite correlated nuclear systems . . . . . . . . . . . . . . . . . . .
5.4 Exploring dissipation with isovector oscillation in a ﬁnite system . . . . . . .
5.5 Galilean Covariance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6 Summary and future prospects . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63
64

66
66
70
76
85
90
92

96

APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

99
APPENDIX A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
APPENDIX B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

vii

LIST OF TABLES

Table 2.1:

Parameters for the mean-ﬁeld interaction. . . . . . . . . . . . . . . .

Table 5.1:

Parameters for the auxiliary ﬁeld in Eq. (5.13).

. . . . . . . . . . .

14

75

viii

LIST OF FIGURES

Figure 1.1:

Figure 1.2:

Schematics of heavy-ion collisions in the center of mass frame. The
left part depicts the two approaching nuclei, and b is the impact
parameter. The dotted lines are boundaries of the overlap regions.
The right part depicts the aftermath of the collisions: the colored
dots in the central region represent the nucleons as well as possible
new particles generated in the participant zone; the two “caps” ﬂying
oﬀ are knowns as the spectators. This ﬁgure is taken from [1].
. . .

Schematic illustration of the eﬀects of ﬂuctuations in HIC. Under
stable conditions without ﬂuctuations on the left panel, the dynamical
trajectories of identically prepared systems tend to be similar, leading
to a unique outcome. The right panel shows that trajectories can
spread out in the prescience of ﬂuctuations during the evolution and
end up in diverse conﬁgurations. This ﬁgure is taken from [2]. . . . .

Figure 2.1: Density proﬁles of stable nuclei 58Ni and 197Au obtained from solu-
. . . . . . . . . . . . . . . . .

tions of the Thomas-Fermi equations.

Figure 2.2: Time evolution of radial density proﬁles for single nuclei 58Ni and
. . . . . . . . . . . . . . . . . . . . . . .

197Au in steps of 40 fm/c.

Figure 3.1: Final rapidity distributions as a function of reduced rapidity for 197Au
+ 197Au at beam energies of 100 MeV/nucleon (upper panel) and
400 MeV/nucleon (lower panel) at an impact parameter b = 7 fm.
Solid curves, dashed curves and the dashed-dot curve correspond to
the Brownian model, BUU-type models, and QMD-type models [3],
respectively.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.2: Final average in-plane ﬂow as a function of reduced rapidity for 197Au
+ 197Au collisions at beam energies of 100 MeV/nucleon (upper
panel) and 400 MeV/nucleon (lower panel) and an impact param-
eter b = 7 fm. Solid curves, dashed curves and the dashed-dot curve
represent the Brownian model, BUU-type models, and QMD-type
models [3], respectively. . . . . . . . . . . . . . . . . . . . . . . . . .

3

8

28

29

32

34

ix

Figure 3.3: Flow slope parameter for diﬀerent transport models [3] for 197Au +
197Au collisions at beam energies of 100 MeV/nucleon (blue squares)
and 400 MeV/nucleon (red triangles) at an impact parameter b = 7
fm. The error bars represent the ﬁtting uncertainties. These become
invisible when they are smaller than the symbols. The colored bands
correspond to roughly 52% conﬁdence intervals from the statistics
of calculations from both the BUU-type and the QMD-type models.
(See text for a more detailed explanation.)
. . . . . . . . . . . . . .

Figure 3.5: Density contours for nucleons projected onto the reaction plane in the
reaction 112Sn + 112Sn at 50 MeV/nucleon at an impact parameter of
b = 0.5 fm, at diﬀerent times during the reaction. Results calculated
with the Brownian motion model are displayed in the ﬁrst row of
the panels. The bottom two rows show results from SMF and AMD
calculations [4]. The densities for the projected contours start at 0.07
. . . . . . . . . . . .
fm−2 and consecutively increase by 0.1 fm−2.

Figure 3.6: Distribution of IMF multiplicity obtained from diﬀerent models for
the reaction 124Sn + 124Sn at 50 MeV/nucleon at an impact param-
eter b = 0.5 fm. The upper panel shows the multiplicity distribution
of IMFs with charge Z > 2 and the lower panel with Z > 6 [4].
. . .

Figure 4.1: Centrality analysis of the events in the simulation. (a) shows the dis-
tribution of the actual impact parameter b sampled in the simulation;
(b) shows the histograms of b for each event group of deﬁnite IMF
multiplicity MIMF; (c) shows the distributions of centrality estima-
tor Et12/Ec.m. for each event group speciﬁed by MIMF; (d) displays
the relationship between Et12/Ec.m. and b within the events that are
color-coded by their corresponding MIMF. . . . . . . . . . . . . . . .

35

38

41

47

Figure 4.2: Charge Z (uppermost panels), parallel velocity Vz (middle panels)

and angular θF/QP distributions for the Ta+Au system at 39.6 MeV/nucl.
The columns correspond to events of diﬀerent IMF multiplicity MIMF.
Fragments are color-coded and sorted by their charge, with Z1 being
the heaviest and ZMIMF being the lightest. Simulation results are
shown as histograms, and data of charge and velocity distributions
extracted from Ref. [5] are shown as solid dots. The theoretical charge
distribution of events with MIMF = 1 is displayed as a semi-opaque
shaded area in the background of the leftmost panel in the ﬁrst row.
(See the text for more discussion.) . . . . . . . . . . . . . . . . . . .

50

x

Figure 4.3: Means and standard deviations of the charge distributions sorted by
the IMF multiplicity (from MIMF = 2 on the leftmost to MIMF = 4 on
the rightmost). Within each panel, the fragments are ranked by their
charge (with the greatest near the top and the lightest near the bot-
tom) and color-coded accordingly. The asterisks represent theoretical
means and the open circles represent experimental means, both ac-
companied by horizontal “error bars” of the size of the corresponding
standard deviation.
. . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.4: Means and standard deviations of the charge distributions sorted by
the multiplicity MIMF in ascending order from left to right and by the
impact parameter b in ascending order from top to bottom. The dis-
tributions within each panel are further ranked by their correspond-
ing charge Z and color-coded. Theoretical means are represented
by asterisks, and experimental ones are by open circles. Standard
. . . . . . . . . . . . . . .
deviations are represented as “error bars”.

Figure 4.5: Means and standard deviations of the distribution of parallel velocity
Vz, organized in a way similar to Fig. 4.3. . . . . . . . . . . . . . . .

Figure 4.6: Density plots of two simulated events at 180 fm/c. Event (a) is cen-
tral, while event (b) is semi-peripheral. The density contours are
ranked by 0.1 ρ0, 0.2 ρ0, 0.4 ρ0, 0.6 ρ0 and 0.8 ρ0 from the outermost
to the innermost. Neither case provides any direct proof of the neck
formation.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.7: Means and standard deviations of the distribution of parallel velocity
Vz, sorted by the impact parameter b in the simulation and organized
in a way similar to Fig. 4.4.
. . . . . . . . . . . . . . . . . . . . . .

55

57

59

60

62

Figure 5.1:

Second-order direct Born diagrams for neutron and proton self-energies
Σn/p. Note that the particle species in the loop can be the same as or
the opposite of the species of the propagator in the line to the left of
the diagram. In the present work, the interaction transfers no isospin
within the Born diagram, making the scattering independent of isospin. 72

Figure 5.2: One-dimensional nuclear equation of state in the NGF approach.
Crosses represent the results obtained within pure mean-ﬁeld approx-
imation. Dashed lines represent a cubic spline ﬁt through the crosses.
Open boxes show results obtained within complete calculations with
correlations, where the pure mean-ﬁeld results were matched by ad-
justing the parameters in the auxiliary ﬁeld in Eq. (5.13).
. . . . . .

75

xi

Figure 5.3: Occupation of the neutron natural orbitals. The open circles repre-
sent those of the correlated ﬁnal state after the adiabatic evolution,
and the crosses correspond to the uncorrelated initial state. The lines
serve to guide the eye. . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5.4: Neutron spectral functions vs energy for an N = Z = 6 self-bound
system in an L = 15 fm box at discretized momenta. The red dashed-
dot line indicates the kinetic energy of a free neutron with momentum
p. See the text for discussion. . . . . . . . . . . . . . . . . . . . . . .

Figure 5.5: Momentum state occupation for neutrons in the same N = Z = 6
ground-state system as in the preceding two ﬁgures and for a Fermi
gas. The dashed, dash-dotted and solid lines represent the Fermi
gas at temperature T = 3.6 MeV, the mean ﬁeld system and the
correlated system, respectively. . . . . . . . . . . . . . . . . . . . . .

Figure 5.6:

Isovector-like oscillation of the neutron slab and the proton slab with
mean-ﬁeld approximation. The red line represents the neutron den-
.
sity proﬁle and the blue line represents the proton density proﬁle.

Figure 5.7: Evolution of the lowest three moments of the neutron slab. The red
lines and the black lines correspond to calculations with mean ﬁeld
. . . . . . . . . . . . . .
only and the full calculations, respectively.

Figure 5.8: Evolution of density for a correlated N = Z = 2 slab boosted with
velocity V = 0.176 c. The initial density is illustrated with a dash-
dotted line in every panel. The panels (b)-(d) show evolution ad-
vanced subsequently over 30 fm/c, with the latest proﬁle illustrated
with a solid line. Intermediate proﬁles at every 15 fm/c are illustrated
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
with dashed lines.

78

83

84

87

88

92

xii

Chapter 1

Introduction

The work presented in this thesis aims to bring new developments to transport theories on

nuclear collisions. To address the phenomenon of multifragmentation, I developed a new

stochastic transport model that treats the dissipation and ﬂuctuations tied to two-body cor-

relations on an equal footing. Towards a fully quantum-mechanical transport theory, I stud-

ied the non-equilibrium Green’s function theory and applied it to describe one-dimensional

nuclear systems.

In the following sections, I will brieﬂy introduce heavy-ion collisions and the transport

models for describing them. The introduction will be followed by a discussion of incorpo-

rating ﬂuctuation and quantum eﬀects into transport models, where both motivation and

some state-of-the-art attempts will be presented. An outline of the upcoming chapters will

be given in the end.

1.1 Heavy-Ion Collisions

In this thesis, we will focus on heavy-ion collisions (HIC) at intermediate energies. These are

very involved nuclear processes, in which a heavy projectile, typically with a mass number

A > 4, impinges on another heavy target at an incident energy ranging from tens of MeVs

up to a few GeVs per nucleon.

The two nuclei come into contact quickly afterwards, and a ﬁnite system composed of

1

hot and dense nuclear matter is formed as the two nuclei compress against one another. The

temperature and density achieved in the compression phase depend on the incident energies

and the alignment of the centers of the two nuclei (also known as centrality). In general, a

higher incident energy and more central collisions will result in hotter and denser matter. By

varying incident energies and projectile-target pairs, HIC allows one to probe the properties

of nuclear Equation of State (EoS), especially those with strong dependence on density, and

is therefore a popular playground for both experimentalists and theorists. The study of

nuclear EoS through HIC is of particular interest to astrophysicists, as the behaviors of the

nuclear EoS around twice saturation density (ρ0 = 0.16 fm−3) aﬀect properties of neutron

stars, such as the mass-radius relation [6] and the tidal deformability [7]. For now, HIC

experiments are the closest thing that one can perform to achieve the high density nuclear

environment of neutron stars on Earth. At the other end of the spectrum concerning the low-

density behaviors of the nuclear EoS, the possible liquid-gas phase transition [8, 9] and the

possible consequences of multifragmentation [10, 11] and formation of exotic structures [12]

also received a fair amount of attention in recent decades.

2

Figure 1.1: Schematics of heavy-ion collisions in the center of mass frame. The left part
depicts the two approaching nuclei, and b is the impact parameter. The dotted lines are
boundaries of the overlap regions. The right part depicts the aftermath of the collisions: the
colored dots in the central region represent the nucleons as well as possible new particles
generated in the participant zone; the two “caps” ﬂying oﬀ are knowns as the spectators.
This ﬁgure is taken from [1].

Fig. 1.1 depicts schematically a simpliﬁed scenario of the collisions of two heavy-nuclei in

the center-of-mass frame, with the initial stage on the left and the nearly-ﬁnal stage on the

right. Such a phenomenological description of HIC is usually known as the ﬁreball model [13].

On the left-hand-side of Fig. 1.1, b is known as the impact parameter, which is a geometrical

measure of the overlap of the two nuclei. The smaller b is, the greater the overlap is. If

one ignores the Fermi motion in the transverse directions, only the nucleons lying within

the overlap region of the two nuclei can interact and slow down and they are known as the

participants. The other nucleons, referred to as the spectators, outside the overlap region

ﬂy away as there are no opposing nucleons to stop them. It is in the participant regions

that nucleons are actively and frequently interacting. Nuclear matter in the central region

will ultimately decompress and fall apart; emitted clusters can be caught by detectors. By

3

measuring the distributions of the ﬁnal particles, clusters and fragments and comparing the

results to theoretical simulations, one may hopefully obtain some constraints on the expected

behaviors of the nuclear EoS.

1.2 Transport Models

HIC diﬀers from few-body direct reactions in terms of the large degrees of freedom involved

as well as the more complex reaction processes with many more steps. The ﬁreball model [13]

mentioned previously is a simple phenomenological model that captures the essence of HIC.

As more details and sophistication were presented by HIC experiments over time, fully

many-body and microscopic theoretical descriptions of HIC were called for.

One of the very ﬁrst microscopic models used for comparison with experimental data is

the cascade model [14]. The cascade model describes both the positions and momenta of all

the nucleons in the system in a time-dependent manner. The interaction among nucleons

comes solely from two-body nucleon-nucleon scattering, which are usually simulated by a

geometric means. Within the same time step, one checks the closest approaches of all

possible nucleon-nucleon pairs. For pairs with a closest approach d <(cid:112)σ(√s)/π, where √s

is the relative energy of the two nucleons and σ is the energy-dependent total cross section,

one samples randomly the ﬁnal momenta of the pair from the underlying angular diﬀerential

cross section. This type of geometric treatment is still used as a paradigm of two-body

scatterings in many of the current transport models.

A serious problem with the cascade model is that nuclei tend to fall apart far too quickly

due to the lack of a potential that binds the nucleons. The cascade plus mean ﬁeld model

came later as an improvement. A popular choice of mean ﬁeld parameterization comes from

4

the Skyrme interactions. The simplest form of a mean ﬁeld potential is U (ρ) = A(ρ/ρ0) +

B(ρ/ρ0)σ, where A < 0, B > 0 and σ > 1. The parameters are adjusted to reproduce

nuclear matter properties. More sophisticated forms of the mean-ﬁeld parametrization as

well as discussion will be presented in Chapter 2. With the inclusion of mean ﬁeld, positions

and momenta of individual nucleons are updated according to Hamilton’s equation (see more

details in Chapter 2), in addition to the geometric treatment of the hard nucleon-nucleon

scatterings.

Another issue with the naive cascade model is the violation of Pauli exclusion principle.

As a remedy, one can impose a Pauli-blocking probability for each hard two-body scattering.

[15] Let the initial momenta of a pair of colliding nucleons be (r1, p1) and (r2, p2), and
the ﬁnal momenta be (r1, p1(cid:48)) and (r2, p2(cid:48)). Note that initial and ﬁnal positions are not
altered in the cascade approach. One builds two spherical phase space regions of radius R in
coordinate space and of radius P in momentum space around (r1, p1(cid:48)) and (r2, p2(cid:48)) and ﬁnds
out the respective phase space occupancy f = N/Nmax, where Nmax = (4π/3)2(RP/h)3 and

N is the number of nucleons in the region under consideration. The scattering is allowed
with a probability of max{(1 − f1(cid:48))(1 − f2(cid:48)), 0}. There is no universal prescription for the
radii in coordinate or momentum space. For a smoother and more accurate description of

the phase space occupancy, one can also consider representing a nucleon by a group of test

particles. Some relevant discussions can be found in Chapter 2.

With the incorporation of mean-ﬁeld dynamics and the Pauli-blocking procedure in ad-

dition to the raw cascade model, we obtain a simulation method based on the Boltzmann-

5

Uehling-Uhlenbeck (BUU) equation.

∂f1
∂t

+

p1
m

∂f1
∂r −

∂U
∂r

∂f1
∂p1

=

g

dσ
dΩ

h3(cid:90) dp2 dΩ v12
×(cid:16) ˜f1 ˜f2f1(cid:48)f2(cid:48) − f1f2 ˜f1(cid:48) ˜f2(cid:48)(cid:17) ,

(1.1)

where f1 = f (r, p1, t) is the phase space density, U is the mean-ﬁeld potential, g is the
dΩ is the NN-diﬀerential cross section, and ˜f = 1− f. The left-hand-side
of the equation governs mean-ﬁeld dynamics, and the right is the two-body collision term,

degeneracy factor, dσ

which describes how binary collisions of the form (1, 2) ←→ (1(cid:48), 2(cid:48)) induce depletion (→)
and accretion (←) of particles at the phase space site 1.

The BUU equation [Eq. (1.1)] serves as the foundation of the two major types of trans-

port models for HIC frequently used currently: the BUU type and the Quantum Molecular

Dynamics (QMD) type. The BUU-type models tend to be more truthful in solving the BUU

equation with the test particle techniques, combining elements of the mean-ﬁeld evolution

and the cascade model. The QMD-type models represent nucleons as rigid wavepackets of

ﬁnite extension and simulate the mean-ﬁeld evolution and the binary collision of the cen-

troids of those wavepackets. Numerical ﬂuctuations in the QMD-type models are generally

not washed out as they are in the case of the BUU-type models, but QMD-type models tend

to suﬀer from issues related to violation of Pauli principle. More discussion regarding the

Boltzmann equation as well as the traditional transport models will be given in Chapter 2.

6

1.3 Fluctuations and Quantum eﬀects

The Boltzmann equation [Eq. (1.1)] represents a truncation up to two-body correlations

in the Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY) hierarchy, and thus is not the

whole story. The deterministic nature of the Boltzmann equation predicts a single averaged

evolution trajectory for HIC, as is illustrated on the left side of Fig. 1.2, which is clearly at

odds with the diverse fragmentation in HIC observed in experiments [16, 17].

The formation of light clusters relies on correlations beyond two body (e.g. 3-body, etc.),

and the production of IMFs is thought to be the result of the spinodal decomposition seeded

by density ﬂuctuations. Incorporating many-body correlations, even on a semi-classical level,

is prohibitively complex. Alternatively, one may introduce the eﬀects of higher-order corre-

lations as ﬂuctuations in phase space. The eﬀects of higher-order correlations scale like ρn,

where n is the order of the correlations, and they are perturbatively small when the density

is low. The erratic and complex behavior of higher-order correlations, although not random

in a strict sense in the semi-classical theory, render them highly sensitive to initial conditions

and hence unpredictable, which justiﬁes a stochastic treatment of the perturbations, leading

to the introduction of ﬂuctuations into transport models.

7

Figure 1.2: Schematic illustration of the eﬀects of ﬂuctuations in HIC. Under stable condi-
tions without ﬂuctuations on the left panel, the dynamical trajectories of identically prepared
systems tend to be similar, leading to a unique outcome. The right panel shows that tra-
jectories can spread out in the prescience of ﬂuctuations during the evolution and end up in
diverse conﬁgurations. This ﬁgure is taken from [2].

The so-called Boltzmann-Langevin (BL) equation [2, 18] has been introduced to include

the ﬂuctuation δI[f ] of the averaged collison integral (cid:104)Icoll[f ](cid:105),

∂f
∂t

+

p
m

∂f
∂r −

∂U
∂r

∂f
∂p

= (cid:104)Icoll[f ](cid:105) + δI[f ].

(1.2)

The ﬂucutation term δI[f ] is typically assumed to be normally distributed with a vanishing
mean and a variance equal to (cid:104)Icoll[f ](cid:105). [2] The inclusion of δI[f ] induces ﬂuctuations in the
phase space density f 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.

The BL equation provides a formalism to incorporate ﬂuctuations, the most truthful

and exact numerical implementation of which is only available in two dimensions [19–21] to

my best knowledge. Some attempts of treating ﬂuctuations [22] in transport theory indeed

predate the birth of BL equation, while many of the other modern stochastic transport

8

Figure1:Illustrationoftheeﬀectofﬂuctuationsonadynamicaltrajectory.Instableconditionsamoderatespreadoftrajectories,aroundtheaverage(associatedwiththeaveragecollisionterm¯Icoll[f])isobserved(leftpanel).TherightpanelshowsthatinpresenceofinstabilitiestheﬂuctuatingBLtermleadstobifurcationoftrajectories.Takenfrom[63].2.3FluctuationsinfullphasespaceandtheBLOBmodelThelatterprocedurecanbeimplementedbyreplacingtheresidualterms(¯Icoll+δI)byasimilarUehlingUhlenbeck(UU)-liketerm,involvingnucleonpackets,whichrespectstheFermistatisticsbothfortheoccupancymeanvalueandfortheoccupancyvariance.ThusonemayconsiderarescaledUUcollisiontermwhereasinglebinarycollisioninvolvesextendedphase-spaceportionsofnucleondistributionofthesametype(neutronsorprotons),A,B,tosimulatenucleonwavepackets,andPauli-blockingfactorsactonthecorrespondingﬁnalstatesC,D,alsotreatedasextendedphase-spaceportions.Thechoiceofdeﬁningeachphase-spaceportionA,B,CandDsothattheisospinnumberiseither1or−1isnecessarytopreservetheFermistatisticsforbothneutronsandprotons,anditimposesthatblockingfactorsaredeﬁnedaccordinglyinphase-spacecellsforthegivenisospinspecies.TheaboveprescriptionsleadtotheBoltzmannLangevinOneBody(BLOB)equations[64]:∂f∂t+{f,H}=¯Icoll+δI==g!dpbh3!W(AB↔CD)F(AB→CD)dΩ.(12)Intheaboveequation,Wisthetransitionrate,intermsofrelativevelocitybetweenthetwocollidingphase-spaceportionsanddiﬀerentialnucleon-nucleoncrosssection:W(AB↔CD)=|vA−vB|dσdΩ.(13)ThetermFcontainstheproductsofoccupanciesandvacanciesofinitialandﬁnalstatesovertheirfullphase-spaceextensions:F(AB→CD)="(1−fA)(1−fB)fCfD−fAfB(1−fC)(1−fD)#.(14)Inpractice,ifthetest-particlemethodisemployed,thephase-spaceportionsA,B,CandDshouldbeagglomeratesofNtesttest-particleseach,andthenucleon-nucleoncrosssectionusedinEq.(13)8models [23, 24] amount to approximate solutions of the BL equation in one way or another.

The stochasticity of most of the QMD-type models is largely statistical in nature, with the

use of wavepackets of ﬁnite extension, whose tie to the BL equation is perhaps even less

clean-cut. As such, the BL equation needs not be treated as the gold standard of stochastic

treatment in semi-classical transport theories. In Chapter 2, I will discuss my treatment of

the dissipation and ﬂuctuation dynamics in HIC, from a rather diﬀerent starting point and

perspective.

So far, we have discussed only semi-classical transport theories, which come with limita-

tions, in spite of their simplicity and success in interpreting HIC physics: 1. New approaches

and techniques in other ﬁelds of nuclear physics are mostly on a quantum mechanical basis,

which makes it diﬃcult to transfer and apply them to the semi-classical transport models;

2. An exact antisymmetrization is also diﬃcult to achieve in semi-classical models, as they

often suﬀer from imperfections in enforcing the Pauli-blocking in binary collisions [25]; 3. It

remains largely unexplored how large a role quantum correlations play in HIC physics.

The Time-Dependent Hatree Fock (TDHF) theory has been developed for low-energy nu-

clear collisions [26,27]. The dissipation in TDHF, however, originates solely from mean-ﬁeld

dynamics, which is insuﬃcient for collisions at higher energies. The Time-Dependent Density

Matrix (TDDM) theory [28], for instance, has been developed to take into account of two-

body correlations, but its application has mostly been restricted to collective vibrations [29]

and fusion [30], due to its complexity. Alternative attempts, such as the Stochastic Mean-

Field (SMF) approach by Ayik et al. [28,31] and the Stochastic-TDHF (STDHF) [32] theory

have also been formulated to incorporate beyond-mean-ﬁeld correlations by introducing ei-

ther initial ﬂuctuations or continuous random noise during the evolution, which simplify the

calculations in principle. Yet, applications of these stochastic quantum approaches target

9

predominantly low-energy nuclear collisions [28]. In Chapter 5, we will present a quantum-

mechanical transport model based on the theory of non-equilibrium Green’s functions. The

model includes short-range two-body correlations to mimic binary collisions in semi-classical

models, which should provide suﬃcient stopping for intermediate-energy HIC.

1.4 Outline

Chapter 2 will give a formal introduction and formulation of the stochastic transport model

based on one-body Langevin dynamics, followed by some applications to heavy-ion collisions

and benchmarks with other models in Chapter 3. In Chapter 4, I will attempt to address

the “hierarchy eﬀect” in multifragmentation observed in experiment with this new model. In

Chapter 5, I will turn to the quantum transport approach based on non-equilibrium Green’s

functions and discuss its application to nuclear systems in one dimension. Finally, a summary

and conclusion will be given in Chapter 6.

10

Chapter 2

One-body Langevin transport model

2.1 Introduction

In central heavy-ion collisions at modest energies, two nuclei approach and collide to form a

composite nuclear system. At very low incident energies and moderate charges, the system

tends to remain fused and de-excites by emission of a few nucleons and light clusters. The

picture gains in complexity as the incident energy increases, since more energy is available

for the system to populate a greater volume of phase space leading up to a plethora of exit

channels. In violent collisions, the composite nuclear system formed is highly excited, and

its evolution is very sensitive to the instabilities present in the system. These instabilities

may deform the shape of the system in phase space, resulting in a breakup into multiple

fragments.

The fragmentation phenomenon was experimentally observed [33] as early as in the late

1970s. Since the early 1990s, more experimental eﬀorts have been devoted to the study of

intermediate-mass-fragment multiplicities [34].

On the other hand, while diﬀerent transport models have been successfully applied to

describe many one-body observables, our understanding and treatment of the fragmentation

mechanism has yet to be reconciled. The inclusion of ﬂuctuations into transport theories

is expected to be of particular importance. Heavy-ion collision experiments make measure-

11

ments over an ensemble of nearly identically prepared colliding systems, and that these

experimental observables reﬂect the distribution of all possible outcomes of that ensemble.

Angular cross sections, for example, are directly obtained from the angular distribution of

the deﬂected particles, but not from any individual, isolated event. On the contrary, a vast

number of the semi-classical transport models are deterministic in nature. These models

predict a single exit channel in principle. This poses little trouble when the underlying dis-

tribution of outcomes is sharp and narrow. Ensemble averages from transport calculations

also converge very well to one-body observables with only weak dependence on channels, such

as those concerning collective ﬂows [35,36]. However, the conﬁgurations in multifragmenta-

tion are obviously heavily dependent on the exit channels, and on the intermediate channels

as well. It requires the transport models to be able to explore a wide range of dynamical

trajectories. The inclusion of ﬂuctuations creates branching points in the evolution of the

system allowing for jumps among diﬀerent states including those susceptible to instabilities.

In general, there are two major types of transport models for simulating heavy-ion col-

lisions. One type of approaches are essentially molecular dynamics of nucleons represented

by single-particle wave-packets, augmented by a phenomenological two-body collision term

of the wave-packets [37–58]. The propagation and scattering of localized wave-packets help

preserve the many-body correlations, and the stochastic treatment of the two-body scat-

tering introduces ﬂuctuations. The other type of approaches aim at directly solving the

Boltzmann-Uehling-Uhlenbeck (BUU) equation, with the system characterized by a one-

body phase space distribution function [23,24,59–75]. Solving the BUU equation yields the

deterministic time evolution of the one-body distribution function leading to a single exit

channel. In recognition of the importance of ﬂuctuations, many eﬀorts have been made to

extend the Boltzmann framework, such as the derivation of the Boltzmann-Langevin equa-

12

tion by Ayik and Gregoire [76, 77], the Stochastic Mean Field (SMF) model by Colonna et

al. [23], and recently the Boltzmann-Langevin-One-Body (BLOB) dynamics by Napolitani

and Colonna [24]. Meanwhile, the inclusion of few-body scatterings in the pBUU model by

Danielewicz [78] also describes the production of light clusters with mass A < 4.

2.2 Formulation

In this section, we will explain the formulation of the Brownian motion model and discuss

the details of the implementation of the simulation code.

2.2.1 The Boltzmann framework

In semi-classical transport theories [15], the nuclear system is often characterized by the

one-body phase space distribution function f (r, p, t). The time evolution of the distribution

function f is approximated with the Boltzmann equation,

∂f
∂t

+ {f,H} = Icoll.

(2.1)

The self-consistent Hamiltonian H encompasses all information about the nuclear mean
ﬁeld interaction as well as Coulomb interaction, while the residual two-body interaction,

mainly nucleon-nucleon scattering, enters through the collision integral Icoll. The Boltz-

mann equation provides us with a simple deterministic model to study heavy-ion collisions

theoretically. Numerical simulations under the Boltzmann framework can be carried out by

means of the test-particle method [15, 79].

13

2.2.2 The mean-ﬁeld dynamics

Neglecting the collision integral in the Boltzmann equation (2.1), we recover the so-called

Vlasov equation,

∂f
∂t

+ {f,H} =

∂f
∂t

+

p
m · ∇rf − ∇rU · ∇pf = 0,

(2.2)

where H = T + V and the mean ﬁeld U = δV/δρ.

The Vlasov equation retains only one-body information. The interaction between any

individual particle and the rest of the system is approximated by a mean-ﬁeld interaction. In

practical calculations, the phenomenological mean-ﬁeld interactions are usually employed,

Un/p(cid:0)ρ(r), δ(r)(cid:1) = A(cid:18) ρ(r)

ρ0 (cid:19) + B(cid:18) ρ(r)
ρ0 (cid:19)D
ρ0 (cid:19) ± 2Siso(cid:18) ρ(r)
ρ0 (cid:19)δ(r),

∇2(cid:18) ρ(r)

+

C
2/3
0

ρ

(2.3)

where δ = (ρn − ρp)/ρ is the isospin asymmetry and parameters A, B, C, D, and Siso,
summarized in Table I, are ﬁtted to reproduce nuclear matter properties at normal density

ρ0 = 0.16 fm−3: the binding energy of 16 MeV/nucleon, the incompressibility of 240 MeV

and the symmetry energy of 30.3 MeV at normal density ρ0 [3, 79]. Spin dependence and

momentum dependence is ignored for simplicity in this parametrization.

A [MeV]

-209.2

B [MeV]

156.4

C [MeV]

-6

D

1.35

Siso [MeV]

18

Table 2.1: Parameters for the mean-ﬁeld interaction.

The Coulomb potential UCoul(ρch(r)) can be determined from the Poisson’s equation for

14

electrostatics,

∇2UCoul = −

1
0

ρch(r).

(2.4)

In the current model, we consider two species of particles only: neutrons and protons.

The numerical scheme of solving the Vlasov equation is adapted from the lattice Hamiltonian

method with test particles proposed by Lenk and Pandharipande [79]. The coordinate space

is discretized into a cubic lattice with the lattice spacing l = 1 fm. Each test particle has a

triangular-shaped form factor and contributes to the nearest eight lattice sites.

2.2.3 The dissipation and ﬂuctuation dynamics

In heavy-ion collisions, nucleon-nucleon scattering acts like a dissipative force, driving the

system towards thermal equilibrium. In accordance with the ﬂuctuation-dissipation theorem,

the dissipation of the beam energy heats up the system and is thus inevitably accompanied

by thermal ﬂuctuations. The thermal ﬂuctuations may manifest themselves in terms of

ﬂuctuations in phase space density, which are expected to be linked with multifragmentation

observed in intermediate-energy heavy-ion collisions. In this subsection, we aim to develop

a framework that oﬀers a consistent and simultaneous description of the dissipation and

ﬂuctuation dynamics.

The collision integral tied to two-body scattering reads

Icoll =

g

h3(cid:90) d3pb(cid:90) dΩ

dσab
dΩ

vab(cid:2)(1 − fa)(1 − fb)fa(cid:48)fb(cid:48)

−fafb(1 − fa(cid:48))(1 − fb(cid:48))(cid:3),

(2.5)

where the degeneracy factor is g = 4 for nucleons (when no distinction is made between

15

neutrons and protons; see the discussion later), and vab is the relative velocity between

nucleons a and b, and dσ/dΩ the diﬀerential nucleon-nucleon cross section. The indices a(cid:48)

and b(cid:48) denote the ﬁnal states of the colliding pair.

As the scattering energy increases, nucleon-nucleon scattering peaks more sharply for-

ward. One can reduce the collision integral Icoll (2.5) into a Fokker-Planck form by making

an expansion over the scattering angle θ (see Appendix 6 for more details) [80],

Icoll → −(cid:88)i

where

∂
∂pi

a(cid:110)fa ·

1

a(cid:3)(cid:111)
2(cid:2) ˜Ri
a + (1 − fa)Ri
a(cid:17),
a(cid:16)faD
+(cid:88)i,j

∂2
a∂p

∂pi

ij

j

Ri
a = −
˜Ri
a = −
g
ij
a =

D

g

g

ab,

h3(cid:90) d3 pb fb Fab qi
h3(cid:90) d3pb fb (1 − fb) Fab qi
4h3(cid:90) d3pb fb(1 − fb)Fab(q2

ab,
abδij− qi
abq

(2.6)

(2.7)

(2.8)

(2.9)

j
ab),

with qab = pa − pb and

Fab = (πvab/2)(cid:90) 1

0

θ2 (dσab/dΩ) d cos θ.

(2.10)

The beyond-mean-ﬁeld dynamics, depicted traditionally as two-body scattering processes,

is transformed, by the Fokker-Planck Eq. (2.6), into diﬀusion processes of nucleons in the

viscous nuclear system. The vector coeﬃcients R and ˜R, usually known as the drag coeﬃ-

16

cients, are connected to the viscosity of the system. The tensor coeﬃcient D is referred to

as the diﬀusion tensor, which describes the anisotropic diﬀusion of particles.

Note that the diﬀerential cross section dσab/dΩ enters the equation through the func-

tion Fab, and hence has an eﬀect on all coeﬃcients in the equation. Admittedly, only the

integrated eﬀects of the cross section can be accounted for, so the overall magnitudes of the

coeﬃcients do not depend sensitively on the detailed angular dependence of the cross section.

The anisotropic diﬀusion is predominantly governed by the anisotropy of the momentum dis-

tribution of particles in the medium, as can be seen in the deﬁnition of the diﬀusion tensor in

Eq. (2.9). One can carry out the expansion to higher orders. The resulted diﬀerential equa-

tion beyond second order, however, is not tractable with the aid of diﬀusion or Langevin

dynamics. Furthermore, similar expansions such as the Kramers-Moyal expansion of the

Boltzmann integral operator in terms of diﬀerential operators beyond second order do not

behave well [81]. Therefore, we stick to the second-order Fokker-Planck equation and choose

to test its validity.

Consider an ensemble of systems with identical initial conditions.

In the presence of

ﬂuctuations, the evolution of the ensemble will diverge. The Fokker-Planck equation provides

a mathematical description of the distribution and ensemble-averaged behavior of these

identically prepared systems.

Indeed, the Fokker-Planck approach has been employed to

study the motion of an ensemble of Brownian particles, classical or quantal, in a medium at

constant temperature, as the stationary solution of the Fokker-Planck equation yields the

equilibrium distribution representing the correct statistics [82].

Inspired by the ideas of Brownian motion in a heat bath, we intend to encapsulate the

beyond-mean-ﬁeld dynamics altogether into the Brownian motion of nucleons in the typically

non-equilibrium nuclear medium, through the Fokker-Planck Eq. (2.6). While the Fokker-

17

Planck equation is deterministic, one may simulate the diﬀerent dynamical trajectories of

the system by use of the corresponding Langevin equation. The diﬀerential form of the

nonlinear Langevin equation for nucleons undergoing Brownian motion reads,

dpa =

1

2(cid:2) ˜R + (1 − fa)R(cid:3)dt + σdBt,

(2.11)

where R and ˜R carry the same deﬁnitions and meanings as in equations (6) and (7), σ is a
3×3 positive deﬁnite matrix such that

Dij =

1

2(cid:88)k

σikσjk.

Bt denotes a Guassian random process with properties

(cid:104)dBt(cid:105) = 0,
tdB

j

(cid:104)dBi

t (cid:105) = dt δij.

(2.12)

(2.13)

(2.14)

This equation describes the momentum transfer, or the “kick”, experienced by a nucleon

due to its interaction with the medium within a time interval ∆t. The ﬁrst term is dissipative

and connected to the viscosity of the nuclear medium, while the second term is stochastic and

gives rise to the ﬂuctuations in the dynamics. In the limit of thermodynamical equilibrium,

coeﬃcients in the Langevin equation are related by the equilibrium temperature, in a manner

akin to the classical Einstein relation, as is shown in Appendix 6.

For the time being, we do not distinguish between nn, pp, and np scatterings by employ-

ing a spin-isospin averaged nucleon-nucleon cross section. We further restrict our attention

to elastic scatterings by targeting collisions at energies near or below pion production thresh-

18

old. Extensions to incorporate elastic and inelastic collisions between diﬀerent species can

be made in the future at little cost by adjusting the degeneracy factor g and adopting the

suitable diﬀerential cross sections dσab/dΩ. More care needs to be taken for the change of

species though in the case of inelastic scatterings.

The Langevin equation meets our goal of treating the dissipation and ﬂuctuations in the

dynamics both consistently and simultaneously. We evolve the system by application of the

Langevin equation to every nucleon in the system in addition to the mean-ﬁeld dynamics.

Details of the numerical implementations will be discussed in a subsequent subsection.

2.2.4

Initialization with the Thomas-Fermi equations

For nucleons inside a stable nucleus, two-body scattering is strongly suppressed by the Pauli

blocking. Hence, for any given mean-ﬁeld potential, the initial conﬁguration of nucleons

in phase space should coincide with the stationary solution to the Vlasov equation. This

solution amounts to that to the coupled Thomas-Fermi equations [79],

Un(cid:0)ρ(r), δ(r)(cid:1) +
Up(cid:0)ρ(r), δ(r)(cid:1) + UCoul +

F(cid:0)ρn(r)(cid:1) = µn,
F(cid:0)ρp(r)(cid:1) = µp,

k2

k2

2
2mn
2
2mp

(2.15)

(2.16)

where UX (ρ, δ) is the self-consistent mean-ﬁeld potential as in Eq. (2.3), µX is the chemical

potential and kF is the Fermi momentum. The subscript X denotes the particle species.

Because of the surface term in the mean ﬁeld, the Thomas-Fermi equations are second-

order ordinary diﬀerential equations. They are to be solved with the boundary conditions

ρn/p(r → ∞) = 0 and (dρn/p/dr)|r=0 = 0. One can solve them numerically by employing
Ansatzes for ρn/p(r) and adjusting µn/p iteratively [66].

19

In this work, we propose a diﬀerent method to solve the coupled Thomas-Fermi equations.

We rewrite the equations by multiplying both sides by density ρX,

hX(cid:0)ρn(r), ρp(r)(cid:1) ρX (r) = µX ρX (r)

(2.17)

with the single-particle hamiltonian hX = UX + 2k2
F /2mX + UCoul δX,p. Eq. (2.17) has
the same structure as the Hartree-Fock equation, allowing us to tackle it as an eigenvalue

problem. Using a discretized position basis, we can obtain a matrix representation for hX

and a vector representation for ρX,

(X)
ij

h

(X)
i

ρ

= (cid:104)ri|ˆhX|rj(cid:105) = hX(cid:0)ρn/p(ri), ρn/p(rj)(cid:1),

= (cid:104)ri|ρX(cid:105) = ρX (ri).

(2.18)

(2.19)

Note that hX is not diagonal, because the derivatives involved are computed in terms of ﬁnite

diﬀerences in the basis. In Eq. (2.17), ρX plays the role of the eigenvector om hX and µX the

eigenvalue. We use a self-consistent iterative method [83] to ﬁnd eigenvectors and eigenvalues

for hn and hp. The pair of eigenvectors {ρn, ρp} in the position basis corresponding to the
smallest eigenvalues, i.e., lowest chemical potentials {µn, µp}, is chosen to generate the ﬁelds
hn/p(ρn, ρp) in the next iteration and picked as the actual density proﬁles in the end.

We illustrate this method by computing the density proﬁles for a medium-sized nucleus

58Ni and a large-sized nucleus 197Au using the mean-ﬁeld potential with a K = 240 MeV

mentioned above. The computed radial density proﬁles are shown in Fig. 2.1. The density

proﬁles with the Thomas-Fermi approximation lacks the ripples associated with shell eﬀects

in typical Hartree-Fock calculations. The tails of the density proﬁles exhibit rapid fall-oﬀs,

20

which is also typical of Thomas-Fermi calculations [83]. In practice, the unphysical fall-oﬀ of

the tails is mitigated by numerical sampling of test particles and coarse graining, as is seen

in the initial distributions at t = 0 fm/c in Fig. 2.2.

2.2.5 Brownian motions of nucleon wave-packets

In this model, the beyond-mean-ﬁeld residual-interactions, between individual nucleons and

the nuclear medium they are locally immersed in, are presumed to be dissipative and random,

and governed by the proposed Langevin Eq. (2.11). We refer to these momentum and energy

exchanges between particles and medium as Brownian motions. Below we will describe

at length the perspective on these Brownian motions and their numerical implementation.

Given the mesoscopic nature of nuclear systems, physics and practical details are entangled.

2.2.5.1 Partition of test particles into nucleon wave-packets

While ensuring a suﬃciently smooth coverage of the phase space in the simulation of mean-

ﬁeld dynamics, the large quantity of test particles used, typically Ntest = 102 - 103 test

particles per nucleon, have adverse eﬀects on the ﬂuctuation dynamics. The scatterings

of test particles supposedly representing the same nucleon are uncorrelated, which would

inevitably wash out most of the ﬂuctuations in the dynamics. This, to a large extent,

explains why BUU-type approaches typically have vanishingly small ﬂuctuations compared

with QMD-type approaches, whose degrees of freedom are nucleons. Diﬀerent attempts have

been made to restore the nucleonic degrees of freedom in two-body scatterings in the BUU

framework [24, 35]. The main idea is to associate test particles adjacent in phase space into

so-called nucleon wave-packets and to move them collectively as a whole.

We adopt a similar approach to enhance the eﬀects of ﬂuctuations. In each time step,

21

we partition the test particles into nucleon wave-packets and execute Brownian motion with

these nucleon wave-packets as the degrees of freedom. The pre-partition is accomplished

through the k-means clustering algorithm [84] with a metric in phase space parametrized in

the following form,

d2 =

(ri − rj)2

d2
r

+

(pi − pj)2

d2
p

(2.20)

where subscripts i and j denote two points in phase space. The parameters dr and dp address

the compactness in the coordinate and the momentum spaces respectively. We run the k-

means clustering algorithm to partition both the neutron test particles and the proton ones

separately. The algorithm is set to terminate after several iterations, and the values dr =

1.2 fm and dp = 130.5 MeV/c are used. It is found in practice that the ﬁnal results are

not sensitive to either the early termination of the clustering algorithm or the values of the

metric parameters.

After the pre-partition, the system is divided into N neutron subspaces and Z proton

subspaces. The pre-partition is simple but somewhat arbitrary, and thus only the centroids

are used. We identify these centroids as the scattering centers for the nucleons.

For each centroid (ri, pi) for which the local nucleon density is above 0.08 fm−3, we
consider a spherical region centered at ri of radius R ∼ 2 fm. This value corresponds
roughly to the sum(cid:112)σNN/π +(cid:113)(cid:104)r2
ch(cid:105), where the nucleon-nucleon cross section σNN (cid:39) 40
mb and the root-mean-square proton charge radius(cid:113)(cid:104)r2
ch(cid:105) (cid:39) 0.86 fm. R can also be made
density-dependent, with the empirical choice of R(ρ) = [1.99 − 0.18ρ− 4
3 (ρ − ρ0)]fm. The
latter empirical formula is what actually used in the code. Inside the spherical region, we

search for test particles close to the centroid (ri, pi) using the following phase-space metric:

22

d2 =

(ri − r(cid:48))2

σ2
r

+

(pi − p(cid:48))2

σ2
p

(2.21)

with σr = R and σp = /2R. This metric emphasizes the compactness in momentum space,

while connected to the Heisenberg uncertainty principle. The Langevin Eq. (2.11) is origi-

nally intended for point-like particles, and hence wave-packets well-localized in momentum

space are preferred. The Ntest test particles of the same species closest to the centroid

form the wave-packet to undergo Brownian motion. The rest of the particles constitute the

medium, with which the wave-packet interacts.

2.2.5.2 Evaluations of the Langevin equation’s coeﬃcients on a lattice

The coeﬃcients R, ˜R and D involve integrals folded over the momentum space, which

requires the knowledge of occupation at diﬀerent momenta. To this end, we construct a three-

dimensional cubic lattice over the entire momentum space inside each spherical scattering

region and we evaluate the occupation at diﬀerent sites.

The lattice spacing Lp needs to be chosen with care to faithfully reﬂect the actual spread

and spacing of the underlying test particles. We use the standard deviation σ

wp
p

of the

momentum of test particles belonging to the nucleon wave-packet as a measure to constrain
p , /(2R)}, where R is
lattice spacing. We normally choose the spacing Lp = max{1.2σ
the radius of the spherical scattering region under consideration. With such constraints, the

wp

values of the spacing Lp typically fall between 100 MeV/c and 140 MeV/c, which ensures a

sensible coarse graining in the momentum space.

Occupation f (p) at each momentum lattice site is evaluated in the same fashion as spatial

densities are on a spatial lattice in mean-ﬁeld dynamics simulations. Test particles have a

23

triangular-shaped form factor, contributing to the eight nearest lattice sites only. Integrals of

R, ˜R and D are computed as summations over all sites of the three-dimensional momentum

lattice.

2.2.5.3 Momentum transfer from the nuclear medium to the nucleon wave-

packet

With a reasonable time step size ∆t ∼ 0.25 - 0.5 fm/c, the ﬁrst term of the Langevin Eq.
(2.11) can be readily calculated. The stochastic term involves a matrix σ, which needs to

be extracted as the square root of the diﬀusion matrix D. Note that, by the deﬁnition in

equation (8), D is a real symmetric positive semi-deﬁnite matrix. It follows that D can be

diagonalized as D = OΛO(cid:62), where Λ is a diagonal matrix and O an orthogonal matrix.

We then represent σ, cf. Eq. (A.12), as

σ = OΛ

1
2 O(cid:62),

(2.22)

where Λ

1
2 is diagonal whose diagonal elements are the unique square roots of the correspond-

ing diagonal elements in Λ. It can be easily veriﬁed that σ constructed according to (2.22)

satisﬁes (8).

With a time step of size ∆t, the diﬀerential notation dBt is interpreted as a 3-dimensional

random vector, whose components are independent Gaussian random numbers. The under-

lying Gaussian distribution has a mean equal to zero and a variance equal to ∆t.

In summary, the momentum transfer ∆p within a time step ∆t in a spherical scattering

24

region is simulated as

∆pa =

1

2(cid:2) ˜R + (1 − fa)R(cid:3)∆t + σ g(0, ∆t)

(2.23)

with g(0, ∆t) being a random vector comprised of 3 independent Gaussian random numbers

sampled with mean = 0 and variance = ∆t.

2.2.5.4 Recoil for conservation of momentum and energy

After a nucleon wave-packet is shifted in the nuclear medium inside a spherical scattering

region, the recoil of the nuclear medium needs to be accounted for in order to preserve the

conservation of total momentum and total energy. The interaction between the nucleon

wave-packet and the nuclear medium is reciprocal.

Indeed, by exchanging the subscripts

a and b in the Langevin Eq. (2.11), one obtains an expression for how the nucleon wave-

packet induces recoil of particles in the medium. Thus, the recoil can be treated in principle

precisely.

On the other hand, owing to the facts that the number of particles involved in the medium

is large and that the recoils are coupled in a nontrivial manner, we instead adopt a collective

and approximate treatment of the recoil eﬀects: the center of momentum of the nuclear

medium is shifted to conserve total momentum, and all particles in the medium are scaled

with respect to the new center of momentum to conserve total energy.

Additionally, it is worth noting that the nuclear medium almost always contains more

than one nucleon. The collective shift-and-scale adjustment, in eﬀect, introduces many-body

correlations in the nuclear medium.

25

2.2.5.5 Pauli-blocking procedure

Within the scattering region, after the kick of the nucleon wave-packet and the adjustment

for recoil in the nuclear medium, we compute again the occupation over the entire lattice

in momentum space. The Brownian motion is ﬁnalized only if none of the occupation at

any lattice site exceeds 1. Otherwise, we deem the Brownian motion unphysical and revert

all changes. This Pauli blocking procedure is eﬀective. For single ground state nuclei, over

97% of the attempted Brownian motions are blocked in the current model prescription, and

stability of the nuclei are demonstrated in Fig. 2.2.

Fig. 2.2 shows the time evolution of the radial density proﬁles for a single 58Ni and a

single 197Au up to 200 fm/c in steps of 40 fm/c, simulated by our Brownian code under

the same controlled conditions speciﬁed in Ref. [3]. Ideally, the density distribution should

remain unchanged over time for single nuclei. In our case, the density proﬁles show only

small scale ﬂuctuations over the entire course of simulation. It indicates, in particular, that

numerical solutions of the Thomas-Fermi Eq. (2.15) approximate the true solutions of the

Vlasov Eq. (2.2) reasonably well. Further, the majority of the spurious large momentum

transfers are eﬀectively blocked by the Pauli-blocking procedure.

2.2.5.6 Summary of implementation of Brownian motions

In each time step, the occupied phase space is partitioned into N + Z subspaces of roughly

equal volume. Scattering regions are constructed spherically around the spatial centroids of

each subspace. These regions are to be examined successively in a random order. Within

each scattering region, a separation of the nucleon wave-packet from the nuclear medium is

made. The Langevin equation is evaluated on a 3-dimensional lattice in momentum space,

and the resulted momentum transfer is applied to the nucleon wave-packet. In observance of

26

the conservation laws, the recoil eﬀects are taken into account through an adjustment of the

momentum distribution of particles in the medium. A Pauli-blocking procedure is applied

in the end to preserve the Pauli exclusion principle.

27

Figure 2.1: Density proﬁles of stable nuclei 58Ni and 197Au obtained from solutions of the
Thomas-Fermi equations.

28

00.050.10.150.20246810r (fm)00.050.10.15 (fm-3)np58Ni197AuFigure 2.2: Time evolution of radial density proﬁles for single nuclei 58Ni and 197Au in steps
of 40 fm/c.

29

00.050.10.150.2     0   40   80 120 160 200Time (fm/c)0510r (fm)00.050.10.150.2 (fm-3)197 Au58 NiChapter 3

Applications to heavy-ion collisions at

intermediate energies

3.1 Benchmarks

In this section, we ﬁrst demonstrate practical applicability of our model to heavy-ion colli-

sions by simulating the Au-Au collisions at both 100 MeV/nucleon and 400 MeV/nucleon,

at an impact parameter of b = 7 fm. We study the nucleon rapidity distribution and the
average in-plane ﬂow (cid:104)px/A(cid:105) and compare our simulation results with those from other trans-
port codes in the code evaluation project of Ref. [3]. After conﬁrming that our model yields

reasonable results for one-body observables, we proceed to investigate its ability to describe

multi-fragmentation processes. To this end, we study the time evolution of the systems:

112Sn + 112Sn and 124Sn + 124Sn at 50 MeV/nucleon at an impact parameter b = 0.5

fm. A preliminary comparison of results with the Stochastic Mean Field (SMF) and the

Antisymmetrized Molecular Dynamics (AMD) models from Ref. [4] is also made.

3.1.1 One-body observables for Au + Au system

Fluctuations associated with Brownian motions enable our model to probe a broader range of

intermediate and ﬁnal channels. It is of interest to study whether the diversity of intermediate

30

and ﬁnal channels aﬀects the description of one-body observables. Our model is applied

to simulate the Au + Au reactions at two incident energies, 100 MeV/nucleon and 400

MeV/nucleon. These speciﬁc reactions were also studied and compared in a transport code

evaluation project under controlled conditions [3]. The same impact parameter b = 7 fm as

in Ref. [3] is chosen. Identical mean-ﬁeld interactions and nucleon-nucleon cross sections as

in Ref. [3] are employed.

3.1.1.1 Rapidity distribution dN/dy

The rapidity distribution dN

dy in the ﬁnal state gives a Lorentz invariant measure of the degree
of stopping of nucleons attained in heavy-ion collisions [3]. Rapidity y is deﬁned as follows

y =

1
2

ln

E + pzc
E − pzc

.

(3.1)

In the non-relativistic limit, where pzc << E, rapidity is roughly proportional to the velocity
component parallel to the beam axis: y ≈ vz/c. The more particles populate the mid-rapidity
(y ≈ 0) region in the center of mass frame, the stronger the stopping eﬀects are.

31

Figure 3.1: Final rapidity distributions as a function of reduced rapidity for 197Au + 197Au at
beam energies of 100 MeV/nucleon (upper panel) and 400 MeV/nucleon (lower panel) at an
impact parameter b = 7 fm. Solid curves, dashed curves and the dashed-dot curve correspond
to the Brownian model, BUU-type models, and QMD-type models [3], respectively.

In Fig. 3.1, we display the ﬁnal rapidity distributions from our calculations accompanied

by results of selected BUU and QMD calculations from Ref. [3]. At 100 MeV/nucleon, there

is a large amount of ﬁlling of the mid-rapidity region, indicating a relatively strong stopping.

While all codes except for AMD exhibit a shallow double-humped structure, diﬀerences in

details of the rapidity distributions are not negligible. This is probably tied to diﬀerences

32

05001000-2-1012y/ybeam0200400dN/dy Brownian pBUU SMFAMD IQMD100 MeV/u400 MeV/uin treating Pauli principles in diﬀerent codes and to the delicate competition of mean-ﬁeld

interaction and many-body correlations at this incident energy [3]. At higher incident en-

ergy of 400 MeV/nucleon, fewer particles populate the mid-rapidity region compared to the

outer regions and the double-peaked feature is more pronounced, due to a weaker stopping

and shrinking Fermi momentum compared to incident momentum per nucleon. General

consistency is found among most calculations, although nonrelativistic models such as the

Brownian motion model and AMD predict mildly stronger stopping.

3.1.1.2 Average in-plane ﬂow (cid:104)px/A(cid:105)
Use of a ﬁnite impact parameter, b = 7 fm, in this study, breaks the macroscopic rotational

symmetry around the beam axis in the system, and therefore anisotropy appears in the

transverse collective momentum distribution. We focus on the average in-plane ﬂow (cid:104)px/A(cid:105),
known as the transverse ﬂow, as a function of the reduced rapidity y/ybeam in the center

of mass frame. When quantiﬁed, the transverse ﬂow is commonly described in terms of an

“S-shaped” curve. The slope at the origin, commonly known as the “slope parameter”, is of

importance. The attractive nature of the mean ﬁeld interaction above the Coulomb barrier,

when combined with a ﬁnite impact parameter, tends to draw the two nuclei closer in the

x direction and to yield a negative slope. Two-body scattering, on the other hand, acts

to deﬂect nucleons, bending the slope towards positive. The slope parameter reﬂects the

competition between the mean ﬁeld and the two-body scattering. [3] Particles in the mid-

rapidity region are expected to come from the compressed region during the collision, and

thus this study can also shed light on the behavior of the nuclear equation of state beyond

normal density.

33

Figure 3.2: Final average in-plane ﬂow as a function of reduced rapidity for 197Au + 197Au
collisions at beam energies of 100 MeV/nucleon (upper panel) and 400 MeV/nucleon (lower
panel) and an impact parameter b = 7 fm. Solid curves, dashed curves and the dashed-
dot curve represent the Brownian model, BUU-type models, and QMD-type models [3],
respectively.

In Fig. 3.2 we show the average in-plane ﬂow for our calculations together with results

from selected transport models from Ref. [3].

In all calculations, the expected S-shaped

curves are produced at both energies. The positive slopes at the origin indicate that the

eﬀects of nucleon-nucleon scattering dominate over those of the mean ﬁeld. At 100 MeV/nu-

34

-40-2002040-1-0.500.51y/ybeam-80-40040 px/A  (MeV/c) Brownian pBUU SMFAMDIQMD100 MeV/u400 MeV/ucleon, the BUU-type models clearly produce greater inﬂections than the QMD-type models.

The prediction from the Brownian motion model lies between them. At 400 MeV/nucleon,

it appears that all ﬁve transport models yield very consistent results in the mid-rapidity

region.

Figure 3.3: Flow slope parameter for diﬀerent transport models [3] for 197Au + 197Au
collisions at beam energies of 100 MeV/nucleon (blue squares) and 400 MeV/nucleon (red
triangles) at an impact parameter b = 7 fm. The error bars represent the ﬁtting uncertain-
ties. These become invisible when they are smaller than the symbols. The colored bands
correspond to roughly 52% conﬁdence intervals from the statistics of calculations from both
the BUU-type and the QMD-type models. (See text for a more detailed explanation.)

The slope parameters at mid-rapidity can be extracted through a linear ﬁt in a small

interval centered at the origin. The values of the slope parameters at two energies for diﬀerent

35

GIBUU-RMDGIBUU-SkyrmeIBLIBUUpBUURBUURVUUSMFBrownianAMDIQMD-BNUIQMDCoMDImQMD-CIAEIQMD-IMPIQMD-SINAPTuQMDUrQMD153045607590105120135150165180dpx/A/dyred (MeV/c)100 MeV/u400 MeV/utransport simulations [3] are summarized in Fig. 3.3. The error bars take into account the

ﬁtting uncertainties only. On top of the simulation results, we also added shaded bands

to indicate regions in which calculations are considered to be statistically consistent with

majority of the BUU-type and QMD-type models. To obtain the statistically consistent
regions, we ﬁrst computed the √2σ-intervals centered at the mean values for the BUU sample

and the QMD sample independently. σ stands for the standard deviation of the sample. The
statistically consistent regions were taken to be the overlap of the √2σ-intervals from the

two samples. If we further assume that the BUU sample and the QMD sample follow an

identical Gaussian distribution, the consistent regions can also be interpreted as roughly 52%

conﬁdence intervals. Note that throughout the statistical analysis, results from the Brownian

motion model were deliberately excluded to avoid any possible bias. Nevertheless, the slope

parameters extracted from the Brownian simulations are found to be statistically consistent

with majority of the other calculations.

The reassuring consistency of results between the Brownian motion model and other

current transport models provides evidence that the one-body Brownian motion picture can

successfully capture the eﬀects of two-body scattering and mean ﬁeld in heavy-ion collisions.

Last but not least, in analyzing results from the Brownian motion model, we averaged the

rapidity and the in-plane ﬂow over 32 independent events, and the averaged results exhibit

good parities as functions of the reduced rapidity. The ﬂuctuation dynamics described by the

Langevin Eq. (2.11) does not automatically preserve forward-backward reﬂection symmetry,

and parity symmetry breaking or other types of symmetry breaking can be observed in

individual events. Since the “directions” of the symmetry breakings are essentially random,

symmetry can be restored by ensemble-averaging over independent runs.

In fact, these

symmetry breakings, resulted from the broadening of dynamical trajectories, can be crucial

36

for the formation of fragments to be discussed in the next subsection. For instance, an

uneven break-up of two nuclei in symmetric central collisions will be likely associated with

an asymmetry in the rapidity distribution.

3.1.2 Fragmentation dynamics with Sn + Sn at 50 MeV/nucleon

In this subsection, we will demonstrate how the ability of our model to probe a plethora of

intermediate and ﬁnal states is connected with the formation of intermediate mass fragments

(IMF) in central heavy-ion collisions. The systems we study are 112Sn + 112Sn and 124Sn +

124Sn at 50 MeV/nucleon at an impact parameter b = 0.5 fm. These systems have already

been studied by Colonna et al. using both SMF and AMD [4]. We will make a preliminary

comparison between our simulation and those in SMF and AMD. A density-dependent and

energy-dependent nucleon-nucleon cross section is used,

σN N (Elab, ρ) = σ

free

N N (Elab) exp(cid:20) − α

ρ/ρ0

1 + Elab/E0(cid:21)

(3.2)

with α = 0.3, ρ0 = 0.16 fm−3, E0 = 150 MeV and a maximum cut-oﬀ at 150 mb. σfree
is taken as the cross section parametrization by Li and Machleidt at zero density [85]. The

N N (Elab)

evolution of the systems is followed up to 280 fm/c after initial contact.

Fig. 3.5 shows the density contour plots from projecting nucleons on the reaction plane for

the 112Sn + 112Sn system at diﬀerent stages of the reaction calculated with diﬀerent trans-

port models. All three models give a qualitatively consistent description of the compression-

expansion dynamics, and cluster structures are formed in the expansion phase. During the

approach and compression up to 40 fm/c, calculations from the three models with ﬂuctu-

ations and many-body correlation eﬀects do not appear to be distinguishable from what

37

Figure 3.5: Density contours for nucleons projected onto the reaction plane in the reaction
112Sn + 112Sn at 50 MeV/nucleon at an impact parameter of b = 0.5 fm, at diﬀerent times
during the reaction. Results calculated with the Brownian motion model are displayed
in the ﬁrst row of the panels. The bottom two rows show results from SMF and AMD
calculations [4]. The densities for the projected contours start at 0.07 fm−2 and consecutively
increase by 0.1 fm−2.

would be expected from conventional transport models without ﬂuctuations. The suppres-

sion of the role of ﬂuctuations is linked to the limited volume of phase space available for

the system to populate, since it is far from thermalized in the early stage. At 120 fm/c,

we can already observe, in all simulations, that the expanding systems turn inhomogeneous.

These inhomogeneities, which would not have existed without ﬂuctuations, provide seeds for

fragmentation, and there are cluster structures forming in the core. As the system continues

to expand, the lumps of matter move away from each other and escape from the central

region. For all three models, sizable fragments can be identiﬁed after 200 fm/c, and changes

between 200 fm/c and 240 fm/c are sporadic and moderate. As a result, we assume that the

38

-20-20-10010-20-10010-20-10010-20-10010-20-10010-20-10010-20-1001020SMFAMDz (fm)x (fm)0 fm/c-100102040 fm/c120 fm/c200 fm/c240 fm/cBrownianconﬁgurations have frozen out by 240 fm/c, and we terminate the simulations at 280 fm/c.

The three models diﬀer quite substantially in details of the predictions for the expansion

phase. In the Brownian motion model calculations, the degree of stopping is comparatively

low and the system tends to expand more along the beam direction. On the other hand,

in AMD, the system expands quickly with a focus around the x-direction, which indicates

a very strong stopping that is also seen in the Au-Au simulations. The relatively isotropic

and slow expansion in SMF can probably be explained by the spinodal decomposition of a

nearly homogeneous source at low density [86]. We also count the number of nucleons in the

“gas” phase (ρ < 1

6ρ0) predicted by our model and compare the number to results of SMF
and AMD [87]. It is found that our model yields more gas-phase nucleons than AMD, but

only slightly fewer than SMF.

Regarding the fragmentation mechanism in our model, in-medium Brownian motions in-

troduce branching points in the dynamical trajectories and allow for “jumps” among a greater

range of intermediate and ﬁnal conﬁgurations. These jumps, stemming from nucleon-nucleon

scatterings, are abrupt and discontinuous in time. The Langevin Eq. (2.11), with its stochas-

tic term introducing discontinuities in time, serves the purpose. By solving the Langevin

equations, we attempt to simulate the abrupt jumps from one n-body conﬁguration to an-

other. As the deviations from the ensemble-averaged trajectory predicted by the Boltzmann

Eq. (2.1) accrue from the jumps, exotic conﬁgurations including those with fragmentation

eventually become accessible. In the quantum-mechanical picture, conﬁgurations are repre-

sented by superpositions of Slater determinants. While mean-ﬁeld evolution is coherent in

principle, the residual incoherent many-body correlations, such as two-body scattering, result

in decoherence and transitions between diﬀerent Slater determinants. Among the stochas-

tic approximations of the quantum many-body problems are the AMD model [37, 38, 88]

39

and the Stochastic Time-Dependent Hartree-Fock (STDHF) theory [32].

In treating the

jumps between diﬀerent conﬁgurations as a stochastic process, the Brownian motion model

is conceptually consistent with these quantal approaches.

40

Figure 3.6: Distribution of IMF multiplicity obtained from diﬀerent models for the reaction
124Sn + 124Sn at 50 MeV/nucleon at an impact parameter b = 0.5 fm. The upper panel
shows the multiplicity distribution of IMFs with charge Z > 2 and the lower panel with Z >
6 [4].

41

IMF multiplicity M0.10.20.30.4BrownianSMFAMDZ > 2024681012IMF multiplicity M00.10.20.30.40P(M)BrownianSMFAMDZ > 6A comparison of the IMF multiplicity between the Brownian motion and the other two

transport models for 124Sn + 124Sn is shown in Fig. 3.6. In our case, fragments are iden-

tiﬁed with a simple coalescence algorithm with a cut-oﬀ density ρc = 0.02 fm−3 [35]. The

distribution is obtained from 100 independent simulations.

For IMFs with charge Z > 2, our multiplicity distribution looks compatible with that

from AMD. Both distributions maximize around multiplicity = 7 or 8 and share a similar

spread. For larger IMFs with charge Z > 6, our calculations yield a slightly lower multi-

plicity, while still predicting essentially the same spread as the others. Since our model, as

well as SMF, predicts more free nucleon emission than AMD, fewer nucleons are available

in the “liquid” phase (ρ > 1

6 ρ0) for the assembly of IMFs. Moreover, due to the large num-
ber of free nucleons and light fragments in our simulations, the eﬀective surface-to-volume

ratio is also higher, which might lead to spurious evaporation of test particles and hinder

the formation of fragments. Nonetheless, the Brownian motion model breaks through the

limitations of traditional Boltzmann transport framework and proves to have great potential

for the description of multifragmentation.

3.1.3 Summary and discussion

In the last two sections, we have reformulated the beyond-mean-ﬁeld dynamics in heavy-ion

collisions in terms of Brownian motions of nucleons in the viscous, out-of-equilibrium nuclear

medium, as opposed to the typical two-body scatterings. The Brownian motions represent,

in eﬀect, the momentum and energy exchange between a nucleon and the nuclear medium

in which it is immersed. They are governed by a set of Langevin equations consisting of a

friction-like term and a stochastic term. This approach describes the dissipation and ﬂuc-

tuation dynamics consistently and simultaneously. Furthermore, each simulation generates

42

a unique dynamical trajectory, enabling us to probe diﬀerent exit channels and obtain the

distribution of possible outcomes of the ensemble.

The details of the numerical simulations, including a new method to initialize stable

nuclei from the Thomas-Fermi approximation, have been presented. We have applied our

model to the time evolution of isolated stable nuclei. The stabilities of the simulations and

the nuclei are well established.

To demonstrate that our model’s ability to describe one-body observables is on par with

that of other transport models, we have studied the ﬁnal rapidity distribution and average

in-plane ﬂow in the reaction 197Au + 197Au at two incident energies and showed that our

results are comparable with those obtained from either QMD-type models or BUU-type

models [3].

We have also investigated formation of fragments in heavy-ion collisions with our model

and conﬁrmed the crucial role ﬂuctuations play in seeding multifragmentation. We have

repeated the calculations of Sn + Sn at E = 50 MeV/nucleon, previously done with the

SMF and the AMD models [4]. As seen from the time evolution of density contours for

nucleons projected on the reaction plane, and all three models depict a fragmented system

with similar general features. Regarding the distribution of IMF multiplicity, We ﬁnd that

the yield of light IMFs with Z > 2 in the Brownian motion model is comparable to that in

AMD, but our yield of large IMFs with Z > 6 is slightly lower than that in SMF or AMD.

So far, we have successfully demonstrated the abilities and potential of the Brownian

motion model to describe various scenarios in heavy-ion collisions at intermediate energies.

Its ability to traverse diﬀerent dynamical trajectories makes it particularly suitable for the

study of multifragmentation, which is beyond the reach of many traditional transport mod-

els. It is also superior to many models with stochastic extensions in that it treats dissipation

43

and ﬂuctuation on an equal footing. More explicit introductions of many-body correlations

possibly with connections to quantal expectation values are under consideration. In the fu-

ture, it is of great interest to confront the optimized Brownian motion model to experimental

data. We have also in mind the goal of studying the fragmentation mechanism. For exam-

ple, we have looked at whether and when the central region of a fragmented system enters

the mechanically unstable region and the ﬁndings seem to favor the spinodal decomposition

mechanism. A more careful inspection is planned.

As a ﬁnal note, while the Fokker-Planck/Langevin approach has been motivated here by

the Boltzmann equation, one can think about circumventing the latter in the future and

deriving the coeﬃcients for the equation directly from microscopic theory, after separating

slow coarse-grained nucleon motion from fast inter-nucleon motion [89].

44

Chapter 4

Investigating “hierarchy eﬀect” in

multifragmentation at Fermi energies

with the Brownian motion model

The INDRA collaboration performed a series of experiments for diﬀerent colliding systems

at Fermi energies and studied the fragmentation of the quasiprojectiles [5]. For many of the

systems studied, the experimental observations seemed to suggest the presence of a “hierarchy

eﬀect” between the charge and the parallel velocity of the intermediate mass fragments (IMF)

emitted by the quasiprojectile (QP): “the ranking in charge induces in average a ranking in

the average parallel velocity, the heaviest fragment is the fastest and is focused in the forward

direction relative to the QP recoil velocity” [5]. The hierarchy eﬀect suggests that memories

of the entrance channel play an important role in collisions at intermediate energies. This

essentially rules out the statistical decay approach under the equilibrium assumption and

calls for a dynamical treatment of the collision and the multifragmentation processes.

We will revisit the hierarchy eﬀect by simulating the collision of Ta+Au at 39.6 MeV/nu-

cleon with the Brownian motion model and compare our results with the experimental ﬁnd-

ings.

45

4.1 Simulation setup: physical inputs and selection of

impact parameter b

The transport model used is the Brownian motion model, as described in Chapter 2, with

essentially identical physical inputs. The mean-ﬁeld potential used is of the same form as

Eq. (2.3) with parameters ﬁtted to reproduce the binding energy E/A = -16 MeV/nucleon,

the incompressibility K = 200 MeV and the symmetry energy Esym = 30.3 MeV at nor-

mal density ρ0. The nucleon-nucleon cross section employed is energy-dependent and with

empirical in-medium modiﬁcations,

σNN(Elab, ρ) = σ

free

NN (Elab) exp(cid:20) − α

ρ/ρ0

(1 + Elab/E0)β(cid:21)

(4.1)

with α = 0.3, ρ0 = 0.16 fm−3, E0 = 137 MeV, β = 1.4 and a maximum cuﬀ-oﬀ at 150
NN (Elab) is a parametrization of the free energy-dependent

mb, similar to Eq.

(3.2). σfree

nucleon-nucleon cross section [85].

Another important input for running the simulations is the selection of the impact pa-

rameter b. The impact parameter is not an experimental observable and therefore can only

be inferred from the relevant centrality observable. In the experiment, the level of violence

of the collisions is measured by Et12/Ec.m. (the subscripts 1 and 2 refer to the charge of
particles), the ratio of the total transverse energies of light charge particles with Z ≤ 2 to the
center of mass energy [5]. The experimental data for the Ta+Au system at 39.6 MeV/nu-

cleon shows that for majority of events, with IMF multiplicity MIMF ranging from 2 to 4

and with the charge of the heaviest fragment Z1 between 10 and 55, the distributions of

Et12/Ec.m. are densest in the range between 10% and 20% [5]. A sensible goal of choosing

46

the impact parameters for the simulations is to mimic the distributions of Et12/Ec.m. to our

best ability.

Figure 4.1: Centrality analysis of the events in the simulation. (a) shows the distribution of
the actual impact parameter b sampled in the simulation; (b) shows the histograms of b for
each event group of deﬁnite IMF multiplicity MIMF; (c) shows the distributions of centrality
estimator Et12/Ec.m. for each event group speciﬁed by MIMF; (d) displays the relationship
between Et12/Ec.m. and b within the events that are color-coded by their corresponding
MIMF.

The distribution of impact parameter we used is a uniform distribution from 0 to bmax ≈
12.7 fm. The actual distribution sampled in the simulation is shown in panel (a) of Fig. 4.1.

The total number of events simulated is 1498. In panel (b) of Fig. 4.1, we sorted the events by

their corresponding IMF multiplicity MIMF. Note that events with MIMF > 4 are extremely

rare and are thus omitted in this and the following panels. Most of the fragmentation events

of the quasiprojectile in the simulation took place with an impact parameter less than 10 fm.

Note, in particular, there is also a surge of fragmentation events, with MIMF ≥ 2, and

47

024681012b (fm)020406080counts024681012b (fm)020406080countsM = 1M = 2M = 3M = 40.10.150.2Et12/ECM020406080100countsM = 1M = 2M = 3M = 4024681012b (fm)0.10.150.2Et12/ECMM = 1M = 2M = 3M = 4a depletion of non-fragmentation events, MIMF = 1, between 4 and 9 fm. Exit channels

corresponding to fragmentation seem to be associated with semi-central or semi-peripheral

centralities. The centrality estimator Et12/Ec.m. sorted by MIMF is summarized in panel

(c) of Fig. 4.1. For MIMF between 2 and 4, the bulk part of the events have their estimators

Et12/Ec.m. to be in the range between 12% and 18%, which is in a reasonable agreement

with what was reported in the experimental data [5]. Hence, this choice of uniform sampling

of the impact parameter b is justiﬁed. In panel (d) of Fig. 4.1, we tested the sensitivity of

the centrality estimator Et12/Ec.m. to the impact parameter b in the simulation. From the

perspective of the simulation, we did not ﬁnd a strong correlation between Et12/Ec.m. and b

when b is under 10 fm. This echoes with the observation in (b) that most of the fragmentation

are limited to events with b under 10 fm. Normally, one would expect a larger Et12/Ec.m.

with a smaller b, as there are more participant nucleons in the reaction zone getting deﬂected

into transverse directions in more central events. However, when fragmentation is present,

the initial center of mass energy can transfer not only into the transverse kinetic energy of

free nucleons, but also into energy necessary to break up the colliding systems and to form

clusters. It is the fragmentation involved that make the expected trend between Et12/Ec.m.

and b less clear. When b is greater than 10 fm and practically no fragmentation occurs, the

expected trend is recovered.

4.2 Comparison with experimental data

In this section, we will ﬁrst compare the charge Z, the parallel velocity Vz and certain angular

θF/QP (deﬁned in details in the subsequent text) distributions of the IMFs between the

experimental data and the simulation results, in the form of histograms. The two subsections

48

will present more details of the charge and parallel velocity distributions by comparisons of

the extracted means and standard deviations and discuss the observation and existence of

the “hierarchy eﬀect”.

Before confronting the experimental data, we point out some of the diﬀerences of the

fragments predicted by the model and to justify the comparison. The fragments produced

in the simulations should be treated as primary fragments, which may undergo secondary

statistical decay before reaching detectors. The simulation of statistical decay of primary

fragments is out of the scope of the underlying transport code and is thus not taken into

account. Still, we would like to discuss some qualitative and semi-quantitative comparisons

between the primary fragments from transport theory and the actual fragments analyzed

in the experimental data. We aim to demonstrate the so-called “hierarchy eﬀect”, due to

the entrance channel eﬀects in the distribution and properties of the IMFs. These entrance

channel eﬀects are not expected to emerge from modiﬁcations from the statistical decay of

primary fragments, but rather as a result of the collision dynamics. Hence, for the purpose

of studying entrance channel eﬀects in multifragmentation, it is justiﬁable to look only at

primary framgments in most cases.

49

Figure 4.2: Charge Z (uppermost panels), parallel velocity Vz (middle panels) and angular
θF/QP distributions for the Ta+Au system at 39.6 MeV/nucl. The columns correspond to
events of diﬀerent IMF multiplicity MIMF. Fragments are color-coded and sorted by their
charge, with Z1 being the heaviest and ZMIMF being the lightest. Simulation results are
shown as histograms, and data of charge and velocity distributions extracted from Ref. [5]
are shown as solid dots. The theoretical charge distribution of events with MIMF = 1 is
displayed as a semi-opaque shaded area in the background of the leftmost panel in the ﬁrst
row. (See the text for more discussion.)

50

MIMF = 2020406080Z020406080countsMIMF = 3020406080Z0204060MIMF = 4020406080Z05101520Z1Z2Z3Z4ExpExpExpExp02468Vz (cm/ns)050100150counts02468Vz (cm/ns)020406002468Vz (cm/ns)05101520-101cos F/QP020406080100counts-101cos F/QP020406080100-101cos F/QP05101520MIMF = 1For the comparison, we will focus on the charge Z distribution, the distribution of the

velocity component Vz in the direction of the beam, and the angular distribution of the IMFs

for events with MIMF between 2 and 4. The angle θFi/QP chosen for the angular distribution
is the angle between the velocity of the ith fragment in the frame of the quasiprojectile (QP)

and the quasiprojectile velocity VQP in the center of mass frame. To be more speciﬁc,
cos (θFi/QP ) is deﬁned as

cos (θFi/QP ) =

(Vi − VQP ) · VQP
|Vi − VQP||VQP|

,

(4.2)

where Vi is the velocity of the ith largest fragment in the center of mass frame and VQP is

the velocity of the constructed quasiprojectile,

VQP = (cid:80)MIMF
(cid:80)MIMF

i=1 ZiVi
i=1 Zi

.

(4.3)

In both the experimental and theoretical analysis, the IMFs were selected as fragments with

charge Z > 2 and a positive velocity Vz > 0 in the center of mass frame. The events were

sorted in terms of their IMF multiplicity. Within each type of events, the IMFs are sorted

by their charge with Z1 being the fragment with the greatest charge and ZMIMF being the
one with smallest charge.

In the ﬁrst two rows of Fig. 4.2, we show both the experimental data and the simula-

tion results for the charge Z distribution and the velocity Vz distribution. The theoretical

distributions are normalized to have the same total counts of events as the experimental

counterparts.

The data and the simulation agree coarsely on a semi-quantitative level. The positions

51

of peaks and the spread of the experimental distributions are roughly reproduced by the

distributions of the simulated primary fragments. An exact match is not expected, as the

comparison is done between the primary fragments predicted in theory and the actual frag-

ments picked up in experiments. The distributions in experimental data tend to be wider

than those predicted by the theory. Statistical decay is likely to help broaden the narrower

theoretical distributions in both charge Z and velocity Vz. For the charge distribution cor-

responding to events with MIMF = 2, two peaks appear in the data for both of the IMFs,

Z1 and Z2, while the model predicts a single peak only. The two peaks in the data may

originate from a wide plateau-like distribution though. The theoretical peak of Z1-fragments

roughly corresponds to the lower charge peak in the data, while the higher charge peak in

the data is missing. This is again believed to be because of the lack of treatment for the

primary fragments, especially those produced in events with MIMF = 1 in the simulations.

We put the charge distribution of the projectile-like remnants for MIMF = 1 events in the

simulations in the background of the charge distribution for MIMF = 2. It is observed that

the MIMF = 1 remnants in the simulations have largely ﬁlled up the missing peaks and areas

around the higher charge peak in data for Z1. Thus, it can be conjectured that taking into
account the statistical decay of the remnants in low multiplicity events with MIMF ≤ 2 will
be likely to help recover the missing peak in the case of MIMF = 2 as well as to contribute

to fragments with charge Z > 50 in high multiplicity events with MIMF > 2.

The experimental claims regarding the velocity and the angular distributions of the Z1-

fragments, in relation to other fragments, are conﬁrmed. For all multiplicities, it is observed

in both experiment and simulations that the largest fragment has a larger velocity component

52

along the beam axis z compared to the other fragments, i.e.,

Z1
z

(cid:104)V

(cid:105) > (cid:104)V

Zi
z (cid:105)

(4.4)

for 2 ≤ i ≤ MIMF. Similarly for the angular distributions, regardless of the multiplicity, the
largest fragment tends to peak forward when the quasi-projectiles break up, i.e.,

(cid:104)cos (θF1/QP )(cid:105) > 1,
(cid:104)cos (θFi/QP )(cid:105) < 1,

(4.5)

(4.6)

for 2 ≤ i ≤ MIMF. Both of the observations suggest that the largest fragment, being
fastest and forward-peaked, retains most of the momentum of the projectile, which serves as

evidence for the presence of entrance channel eﬀects for the fragmentation of the system Ta

+ Au at 39.6 MeV/nucleon, in agreement with the experimental observation of the “hierarchy

eﬀects” [5].

It remains to be investigated whether a similar positive correlation between charge and

velocity exists for the fragments Z2 − ZMIMF, which is not immediately clear and evident
from the histograms. In what follows, we will compare the means and the standard deviations

extracted from the distributions in Fig. 4.2, which provide us with succinct abstractions of

the underlying distributions. Extra care should be taken though when the distributions have

more peaks than one, as the mean is no longer indicative of the positions of the peaks.

53

4.2.1 Charge distribution

In Fig. 4.3, we show the means and standard deviations extracted from the theory and the

data for events with MIMF between 2 and 4. As is mentioned before, the charge distribution

of the Z1-fragment in the data for MIMF = 2 has two peaks, so the mean extracted in that

case does not correspond to any of the peak in the data but rather lies somewhere between

two peaks. The mean extracted from the simulation for Z1 for MIMF = 2 mostly corresponds

to the lower charge peak in experiment and is hence to the left of the experimental mean.

For MIMF = 2 and 3, the experimental means are reproduced by the theoretical simulations

reasonably well, and the deviations do not exceed 2 protons in most cases. The theoretical

standard deviations tend to be narrower than their experimental counterparts, which suggest

that the theory still fails to probe the entire variety of exit channels accessible in the actual

collisions. We reiterate here that parts of those missing exit channels might be recoverable

by feeding the simulated primary fragments to secondary statistical decay.

54

Figure 4.3: Means and standard deviations of the charge distributions sorted by the IMF
multiplicity (from MIMF = 2 on the leftmost to MIMF = 4 on the rightmost). Within each
panel, the fragments are ranked by their charge (with the greatest near the top and the
lightest near the bottom) and color-coded accordingly. The asterisks represent theoretical
means and the open circles represent experimental means, both accompanied by horizontal
“error bars” of the size of the corresponding standard deviation.

In Fig. 4.4, we categorize the simulated events further by their impact parameter b

in the range from 0 to 10 fm. Each row of Fig. 4.4 corresponds to a 2-fm range of impact

parameter. However, note that, since it is impossible to diﬀerentiate the experimental events

by b or any other centrality measure with what is available in Ref. [5], each column has the

identical experimental means and standard deviations in all panels. The standard deviations

for the simulations are comparatively narrow, as they take into account the variations among

events within the narrow 2-fm impact parameter range, and thus they should not be directly

compared with the experimental standard deviations, which cover the entire spectrum of

centrality. As the impact parameter b increases, i.e., as the collisions get more peripheral,

the average charge of Z1 shifts towards greater charge value with greater resemblance to the

original projectile. It is thus expected that the entrance channel eﬀects for Z1 grow as b

grows. One should also notice that, even in the most central range of impact parameters

55

0204060ZZ4Z3Z2Z1low Z <----------------------> high ZMIMF = 20204060ZMIMF = 3020406080ZMIMF = 4Z1ExpZ2ExpZ3ExpZ4Exp[0, 2] fm, it is still possible to distinguish the largest fragment with Z1 from all the other

fragments with smaller charge. This shows that the system does not form a compound

nucleus at thermal equilibrium even in the most central collisions. The Z1-fragment can

largely be regarded as a remnant of the projectile. Other fragments with smaller charge are

much less distinguishable among themselves, which suggests that they possibly have a similar

origin that is diﬀerent from the largest fragment, especially for the more peripheral events.

A relatively good agreement for the mean charges between the theory and the experiment is

obtained for the range of impact parameter between 4 and 8 fm for MIMF = 3 and 4. That

interval of impact parameter also overlaps largely with the interval of b, [4, 9] fm, with a

increase of fragmentation events in the simulation.

56

Figure 4.4: Means and standard deviations of the charge distributions sorted by the multi-
plicity MIMF in ascending order from left to right and by the impact parameter b in ascending
order from top to bottom. The distributions within each panel are further ranked by their
corresponding charge Z and color-coded. Theoretical means are represented by asterisks,
and experimental ones are by open circles. Standard deviations are represented as “error
bars”.

57

Z4Z3Z2Z1MIMF = 2MIMF = 3MIMF = 4Z4Z3Z2Z1Z4Z3Z2Z1Z4Z3Z2Z10204060ZZ4Z3Z2Z10204060Z020406080ZZ1ExpZ2ExpZ3ExpZ4Exp[0, 2] fm[2, 4] fm[4, 6] fm[6, 8] fm[8, 10] fm4.2.2 Velocity distribution

We display both the theoretically and the experimentally extracted means and standard

deviations of the velocity Vz distributions, sorted by MIMF as before, in Fig. 4.5. Just like

in the case of charge Z, the theoretical standard deviations are narrower than the experi-

mental ones, which again suﬀers from the fact that only primary fragments are produced

in the simulations. The agreement of the extracted means between theory and experiment

is acceptable, when one considers that we did not adjust any parameters in the model to

force a ﬁt to the data. For the largest fragment with Z1, the theory tend to predict a faster

velocity Vz along the beam axis. The over-prediction gets even more serious when the mul-

tiplicity MIMF gets larger. Vz of the other fragments with smaller charge are accordingly

underpredicted in the simulations, following from the conservation laws. This indicates that

the collisions simulated tend to have a weaker stopping than the actual collisions in the

experiment. Possible reasons for the weaker stopping include our choice of range of impact

parameters and our parametrization of the nucleon-nucleon cross sections in Eq. (4.1).

Putting the issue of somewhat weaker stopping aside, we observe similar trends and

properties for the velocity distributions across the simulations and the data. In particular,

the Z1-fragment has the greatest Vz, which distinctively diﬀerentiates it from other fragments

with smaller charge as well as smaller Vz. Such a divide between the Z1 fragment and the

other fragments has already been seen in the previous discussion about charge distribution.

58

Figure 4.5: Means and standard deviations of the distribution of parallel velocity Vz, orga-
nized in a way similar to Fig. 4.3.

We attribute these divides in the form of charge and velocity to the fact that these frag-

ments have two distinct sources: the projectile-like remnant and the mid-rapidity remnant.

The projectile-like remnant retains most of the entrance channel memories. The mid-rapidity

source is conjectured to be linked to the formation of a neck in the experimental paper [5].

In the simulations we have not been able to gather enough information and evidence to

support the transient existence of a neck. In Fig. 4.6, we show two events in the simulation,

with one corresponding to a central event and the other to a semi-peripheral event. In the

more central event (a), the projectile and the target were strongly deformed, with a lot of

mass accumulating in the central zone. The nucleons in the central zone were substantially

decelerated and formed clusters as they expanded and escaped from the central zone. In

the semi-peripheral event (b), the projectile and the target passed through each other much

more swiftly, leaving behind a thinning tail. The depletion of matter in the central region

also occurred much more rapidly. Both the thinning tails and the matter in the central

zone belong to the mid-rapidity region. Neither scenario in (a) or (b) would ﬁt exactly the

59

024Vz (cm/ns)Z4Z3Z2Z1low Z <------------------------> high ZMIMF = 2024Vz (cm/ns)MIMF = 3024Vz (cm/ns)MIMF = 4Z1ExpZ2ExpZ3ExpZ4Exp“neck” formation conjecture. It could be that the neck is just far too transient and fragile to

be picked out or that many competing fragmentation mechanisms were involved. Note also

that the two pictures presented here by no means exhaust all scenarios that occurred in the

simulation.

A rather weak hierarchy eﬀect, i.e., a positive correlation between charge Z and velocity

Vz, can be seen in both the simulations and the experimental data for all the fragments for

all multiplicities in Fig. 4.5. This hierarchy eﬀect, albeit weak, serves as a manifestation of

the entrance channel eﬀects in play nonetheless.

(a) Simulation with b ≈ 3 fm.

(b) Simulation with b ≈ 9 fm.

Figure 4.6: Density plots of two simulated events at 180 fm/c. Event (a) is central, while
event (b) is semi-peripheral. The density contours are ranked by 0.1 ρ0, 0.2 ρ0, 0.4 ρ0, 0.6
ρ0 and 0.8 ρ0 from the outermost to the innermost. Neither case provides any direct proof
of the neck formation.

In Fig. 4.7, we further divide the events based on their impact parameters. The theo-

retical standard deviations are narrow because of the small range of impact parameters in

each group of events. Other than events corresponding to MIMF = 4 and b ∈ [0, 2] fm,

60

which suﬀers from a small sample size, all other event groups exhibit the “hierarchy eﬀect"

similar to those in the experimental data. Such ranking eﬀects seem to be strongest for

semi-central events with impact parameter b ranging from 2 to 6 fm, particularly for highly

fragmented events. The strong ranking eﬀects observed in the more central events are mani-

festations of the entrance channel eﬀects and can potentially be attributed to the decay of a

strongly deformed and out-of-equilibrium source, such as the fracture of a “neck” connecting

the projectile and the target as proposed in the experimental paper [5]. For greater impact

parameters with b > 8 fm, the divide between the Z1 fragment and the others is more promi-

nent, and the velocity gradient of the smaller fragments also appear much weaker. As the

collisions get more peripheral, the smaller fragments become more homogeneous in terms of

both charge and velocity. One may infer that fragments in peripheral collisions are more

likely to come from a nearly equilibrated source possibly originating from the small and

constrained reaction zone at mid-rapidity.

61

Figure 4.7: Means and standard deviations of the distribution of parallel velocity Vz, sorted
by the impact parameter b in the simulation and organized in a way similar to Fig. 4.4.

62

Z4Z3Z2Z1MIMF = 2MIMF = 3MIMF = 4Z4Z3Z2Z1Z4Z3Z2Z1Z4Z3Z2Z1024Vz (cm/ns)Z4Z3Z2Z1024Vz (cm/ns)024Vz (cm/ns)Z1ExpZ2ExpZ3ExpZ4Exp[0, 2] fm[2, 4] fm[4, 6] fm[6, 8] fm[8, 10] fm4.3 Absence of the “hierarchy eﬀect” for U + C

As is described in Ref. [5], no “hierarchy eﬀect” was reported for the U + C system at the

lab energy E = 24 MeV/nucleon, especially when the multiplicity MIMF is low.

We investigated the same collision system within our model. We noticed in the simu-

lations is that, regardless of the impact parameter, the two incoming nuclei tend to fuse

together to form a compound nucleus with minor deformations and non-negligible rotational

motion. The initial kinetic energy appears to be dissipated by emission of free nucleons and

in the form of internal and rotational energy of the compound nucleus formed. Within the

time scale we considered, up to 300 fm/c, we observe no fragmentation of the compound

nucleus and thus, we are unable to produce the charge and velocity distribution of the IMFs,

as we did for the case of Ta + Au at 39.6 MeV/nucleon. Nonetheless, we speculate that

the compound nuclei predicted by our simulations are susceptible to statistical decay, given

enough time and accuracy of the description of the energy and structure of the compound

nucleus, both of which are lacking in our semi-classical transport approach.

The absence of the correlation between the ranking of charge Z and velocity component Vz

for the U + C system in experiments is more characteristics of the ﬁssion of an equilibrated

source. This is, in fact, in agreement with the theoretical prediction of a single primary

source.

The strikingly diﬀerent behaviors between the Ta + Au system and the U + C in the

simulations can be understood as results of their diﬀerent dominant decay channels: for Ta

+ Au, where the projectile and target are comparable, fragmentation dominates; for U + C,

the ﬁssion channels dominate, due to the high ﬁssility of U. Fission or other sorts of decay of

an equilibrated system typically require a much longer time scale and a rather sophisticated

63

description of the system and is thus beyond the scope of our model.

4.4 Summary

In this chapter, we investigated the behaviors of heavy-ion collisions at Fermi energies and,

in particular, the entrance channel eﬀect reported experimentally [5] for the IMFs. We sim-

ulated the primary fragments in the Ta + Au system at 39.6 MeV/nucleon. Despite the fact

that we did not reproduce exactly the ﬁnal fragments that were picked up in the experiments

for a lack of secondary statistical decay treatment, we ﬁnd that the overall properties and

trends of the experimental fragment distributions also survive in the distributions of the

simulated primary fragments. In particular, we are able to conﬁrm the so-called “hierarchy

eﬀect” in the simulation, in sharp contrast to what might be expected from decay from an

equilibrated source. The fact that the IMF with the greatest charge also travels fastest can

be attributed to the retention of the projectile memory. For certain semi-central collisions,

we also observe a strong correlation between the charge ranking and the velocity ranking

even for the smaller fragments, which might be seen as circumstantial supporting evidence

for the “neck formation” conjectured by the experimental eﬀorts [5].

We also brieﬂy looked at the U + C system at 24 MeV/nucleon, for which the experiment

collaboration reported no “hierarchy eﬀect”. Our simulations of U + C does not predict any

sort of fragmentation or ﬁssion, but rather a rotating compound source. This, however, is

not entirely at odds with the experimental observation, as ﬁssion or statistical decay of a

nearly equilibrated source is indeed consistent with the lack of “hierarchy eﬀect” observed.

Finally, this project, to our best knowledge, for the ﬁrst time, attempts to address the

“hierarchy eﬀect” from a dynamical standpoint with the use of a stochastic transport theory.

64

Our simulations conﬁrm the “hierarchy eﬀect” as an important entrance channel eﬀect for

collisions at Fermi energies, when the multifragmentation mechanism dominates.

In this

study, we did not adjust either the mean-ﬁeld potentials or the N N-cross sections. It will

be interesting to test whether there is any sensitivity of this eﬀect to these inputs of nuclear

interactions in the future. For a more thorough and detailed comparison with data, it will

also be important to combine the use of transport simulations and that of statistical decay

models.

65

Chapter 5

Nonequilibrium Green’s function

approach for many-body quantum

systems

The work in this chapter is now published in Ref. [90].

5.1 Introduction

Quantum mechanics has taken root in contemporary nuclear theory, from the Hartree-Fock

method [91–93] to the chiral eﬀective ﬁeld theory [94–96] in nuclear structure, from the

distorted-wave Born approximation [97, 98] to the coupled-channel techniques [99–101] in

low-energy nuclear reactions, just to name a few. The development of time-dependent quan-

tum approaches is not on par with their static counterparts, with the most notable one

being the time-dependent Hartree-Fock (TDHF) theory [26,27,102–105]. However, the bare

form of the TDHF theory is a mean ﬁeld theory neglecting many-body correlations, which

does not provide suﬃcient stopping and thermalization to describe intermediate or high en-

ergy heavy-ion collisions [106]. The community of intermediate-energy heavy-ion collisions

continues to rely predominantly on semi-classical transport models based on the Boltzmann-

Uehling-Uhlenbeck (BUU) equation, where quantum eﬀects such as Pauli-blocking [107,108],

66

the uncertainty principle, through the use of wavepackets [109, 110], and antisymmetriza-

tion [111, 112], can be partially incorporated. While a vast majority of the semi-classical

transport models have demonstrated their applicability and strengths through comparison

with experimental data [66,113–115], questions concerning the validity of the quasi-particle

picture and the eﬀects of quantum correlations linger. Moreover, as already pointed out in

Ref. [116], the semi-classical methods, detached from the rapidly developing quantum tech-

niques employed in other ﬁelds of nuclear physics, are diﬃcult to improve systematically.

Even attempts to take into account the missing classical many-body correlations through

the incorporation of ﬂuctuations [18, 23, 117, 118] are not able to recover all eﬀects of quan-

tum correlations, as they are of fundamentally diﬀerent nature. A fully quantum-mechanical

transport model for heavy-ion collisions is needed.

Solving the many-body Schrödinger equation exactly is almost always a forbiddingly

daunting task. Fortunately, more often than not, one is interested in one-body observables,

for which a complete many-body wavefunction is not necessary. One can obtain enough

information from one-body reductions, such as one-body density matrix or one-body Green’s

function, by integrating out the irrelevant degrees of freedom. In obtaining a closed equation

of motion, the TDHF theory approximates the two-body density matrix by products of

one-body density matrix and consequently neglects all two-body correlations [26, 83].

In

contrast, in the non-equilibrium Green’s function (NGF) theory, two-body correlations are

encapsulated in the self-energy and a systematic order-by-order approximation of the self-

energy is possible [116, 119]. As a result, the NGF theory is more versatile and promising

than the TDHF theory in the description of intermediate-energy heavy-ion collisions, where

eﬀects of two-body correlations are usually non-negligible.

An NGF code for the description of nuclear collisions in one dimension was developed

67

by Rios et al., and the mean ﬁeld dynamics within NGF was discussed in Ref. [116]. The

current work serves as an extension to the previous work. Before we present the extensions

and discuss the new results, let us brieﬂy go over some relevant concepts and deﬁnitions in

the NGF theory. (See Refs. [116, 119] for more details.) Within the NGF theory, fermionic

many-body systems can be characterized in terms of the single-particle Green’s functions

G<(x1, t1; x2, t2) = i(cid:104)φ†(x2, t2) φ(x1, t1)(cid:105) ,
G>(x1, t1; x2, t2) = −i(cid:104)φ(x1, t1) φ†(x2, t2)(cid:105) .

(5.1)

(5.2)

Here φ and φ† are annihilation and creation operators in the Heisenberg picture, and the

expectation value is taken with respect to the total wavefunction also in the Heisenberg

picture. When t1 = t2, the same-time lesser Green’s function G< is equivalent to the one-

body density matrix up to a constant imaginary factor i. Thus, the lesser and greater Green’s

functions are often associated with the particle density and the hole density, respectively.

One may also deﬁne Green’s functions as functions of momentum and time arguments, by

replacing the annihilation and creation operators in coordinate space by those in momentum

space. Green’s functions in coordinate space and Green’s functions in momentum space can

be converted into one another through Fourier transformation.

The equations of motion for the lesser and greater Green’s functions are the so-called

68

Kadanoﬀ-Baym equations [120],

(cid:20)i ∂

∂t1

2
2m

+

(cid:20) − i ∂

∂t2

2
2m

+

≶

∂2
∂x2

1(cid:21) G

≶

∂2
∂x2

2(cid:21) G

(1, 2) =(cid:90) dx3 ΣHF(1, 3) G

≶

(3, 2)

≶
d3 Σ+(1, 3) G

(3, 2)

t0

t0

+(cid:90) t1
+(cid:90) t2
(1, 2) =(cid:90) dx3 G
+(cid:90) t1
+(cid:90) t2

t0

t0

≶

(1, 3) G−(3, 2) ,

(5.3)

d3 Σ

≶

(1, 3) ΣHF(3, 2)

≶
d3 G+(1, 3) Σ

(3, 2)

≶

d3 G

(1, 3) Σ−(3, 2) ,

(5.4)

where we have used shorthand notations such as 1 to represent (x1, t1). The diﬀerent Σ(1, 2)

are the so-called self-energy, and provide an approximate description of correlations in the

system. They can be expanded diagrammatically, and will be discussed further in the next
≶). The retarded (+) and

(5.10) and (5.11) for Σ

section (see Eq.

(5.9) for ΣHF, Eq.

advanced (−) functions are deﬁned through

F±(1, 2) = F δ(1, 2) ± θ[±(t1 − t2)] [F >(1, 2) − F <(1, 2)] ,

(5.5)

where F δ represents a singular part at t1 = t2. It has been shown that one can derive the

BUU equation from the Kadanoﬀ-Baym equation by approximating the interacting Green’s

functions with the free Green’s functions in a uniform system [119], which demonstrates a

strong link between the NGF approach and the semi-classical models based on the BUU

equation. On the other hand, in contrast to classical transport theories, the correlation

integrals, integrating over the entire history, lead to memory eﬀects.

69

To advance the application of these techniques to nuclear physics, we draw on rigorous

NGF results and combine them with phenomenology taking into account the context of

numerical limitations and the expectations regarding physical characteristics of the nuclear

systems. The aim of this work is to show the applicability of NGF in a variety of settings that

are relevant for nuclear theory, in preparation for the application to simulations of nuclear

collisions. The organization is as follows.

In Sec. 5.2, we will discuss the introduction

of isospin degrees of freedom and the incorporation of short-range two-body interactions.

In Sec. 5.3, we will describe the preparation of ﬁnite correlated nuclear systems within the

NGF approach. Sec. 5.4 examines the isovector oscillation of slabs in one dimension and

Sec. 5.5 demonstrates how to boost a ﬁnite system to move at a constant velocity. In the

end, we will give a summary and remark on the prospects in Sec. 5.6.

5.2 Extending the non-equilibrium Green’s function ap-

proach for nuclear systems

To describe realistic nuclear systems at moderate energies, an obvious but nonetheless crucial

extension to the previous NGF model [116] is to introduce isospin dependence in the Green’s

functions. This can be readily achieved by introducing two versions of the Green’s functions,

of the self-energies, and of other relevant quantities, one for the neutron subsystem and the

other for the proton subsystem. The separation is rather straightforward in practice, but

comes with the cost of doubling the memory requirement for the computation and doubling

the number of Kadanoﬀ-Baym equations to be solved simultaneously. The mapping from

1D neutron or proton densities to 3D densities continues [116] to follow the relation from the

Hugenholtz-van Hove theorem: ρ3D

n/p = ξ ρ1D

n/p with ξ ≈ 0.1217 fm−2. In what follows, the
70

subscripts n and p, indicating the particle species, will be dropped for the sake of brevity

whenever no confusion arises.

The bulk part of mean-ﬁeld interaction U independent of isospin follows closely from that

in Ref. [116],

U (ρ) =

3
4

t0ρ +

2 + σ

16

t3ρ1+d,

(5.6)

with the parametres t0 = −2150.1MeV fm3, t3 = 14562MeV fm3(1+d) and d = 0.257 ﬁtted to
reproduce properties of symmetric nuclear matter at saturation density. With the mean-ﬁeld

U dependent on the total baryon density, the evolutions of the neutron Green’s functions and

the proton Green’s functions are coupled right from the very beginning. Upon diﬀerentiating

the isospin degrees of freedom, we also introduce a dependence on isospin imbalance, δ =
(ρn− ρp)/ρ, into the nuclear equation of state where the form valid up to the second order in
the imbalance is E
A (ρ, 0) + S(ρ) δ2, with S(ρ) = Skin(ρ) + Sint(ρ). The ﬁrst term
Skin(ρ) stems from the diﬀerence in kinetic energies of the neutrons and protons for ﬁnite

A (ρ, δ) = E

imbalance, while the true isospin dependence of the nuclear forces is reﬂected in the second

term taken here in the form

ρ0(cid:19)σ
Sint(ρ) = S0(cid:18) ρ

,

(5.7)

where the values S0 = 20.1 MeV, ρ0 = 0.16 fm−3 and σ = 0.35 are currently used in the
sym
n/p (ρ, δ) is given by

model. The isospin-dependent component of the mean ﬁeld potential U

U

sym

n/p (ρ, δ) = ∂[ρ Sint(ρ) δ2](cid:46)∂ρn/p.

(5.8)

The Hartree-Fock self-energy ΣHF is approximated as a sum of an isospin-independent com-

71

ponent and an isospin dependent component,

ΣHF(1, 2) = δ(1 − 2){U [ρ(1)] + U

sym
n/p [ρ(1), δ(1)]}.

(5.9)

The introduction of explicit isospin degrees of freedom to the system and of the isospin

dependence to the mean ﬁeld paves the way for future study on the nuclear symmetry

energy as well as other isospin-dependent eﬀects [121].

Figure 5.1: Second-order direct Born diagrams for neutron and proton self-energies Σn/p.
Note that the particle species in the loop can be the same as or the opposite of the species
of the propagator in the line to the left of the diagram. In the present work, the interaction
transfers no isospin within the Born diagram, making the scattering independent of isospin.

One of the advantages of the non-equilibrium Green’s function approach, over the more

conventional TDHF theory, lies in its ability to incorporate two-body interactions into the
≶ in the correlation integrals. The incorporation of two-body inter-

residual self-energies Σ

actions induces quantum correlations in the system, which have been lacking in both mean

ﬁeld approaches and semi-classical approaches at the explicit level. In the current approach,

we have included in the self-energies the correlation contribution represented with the dia-

grams in Fig. 5.1 by invoking the second-order Born approximation. Eﬀectively, this direct

diagram describes the scattering within nn, pp and np pairs with diﬀerent lines representing

72

nn/pn/ppp/np/nGreen’s functions of the suitable species. The corresponding self-energies Σ

≶ read

(p, t; p(cid:48), t(cid:48)) =(cid:90) dp1

2π
≶
×G

dp2
2π V (p − p1) V (p(cid:48) − p2)
≶
(p1, t; p2, t(cid:48)) Π

(p − p1, t; p(cid:48) − p2, t(cid:48)) ,

(5.10)

≶

Σ

where

≶

Π

(p, t; p(cid:48), t(cid:48)) = 2(cid:90) dp1

2π

≶
dp2
2π G

≷
(p1, t; p2, t(cid:48)) G

(p2 − p(cid:48), t(cid:48); p1 − p, t) ,

(5.11)

where V (q) is the two-body potential with q being the relative momentum and the factor of

2 accounts for spin degeneracy.

While some previous attempts were made to include residual two-body interactions [122–

124], the choice of the two-body potential V (q) in those attempts might not have been

optimal for the one-dimensional nuclear systems under consideration. Instead, to reﬂect the

short-range nature of the residual nucleon-nucleon interactions and to mimic the eﬀects of

nucleon-nucleon scattering in the kinetic limit, we adopted a new form for the two-body

potential:

V (q) = V0 |q| exp(cid:26)(cid:104) −(cid:16) ηq

2(cid:17)2(cid:105)(cid:27) ,

(5.12)

where the values V0 = 205.491 MeV and η = 0.57 fm have been chosen to establish a cor-

respondence between the one-dimensional collision rates and the three-dimensional collision

rates in symmetric nuclear matter in the semi-classical limit. Note that this potential yields

no diﬀerence in the scattering between np and nn or pp pairs.

The residual interactions introduce an extra contribution termed Ecorr to the total energy

and thus eﬀectively alter the equation of state. This issue arises, as we have ﬁtted the nuclear

73

matter in the mean-ﬁeld approximation to the nuclear matter equation of states in the ﬁrst

place, before introducing two-body interactions. Furthermore, in principle, the interactions
≶ should be the

entering the Hartree-Fock self-energy ΣHF and the residual self-energy Σ

same, which is beyond the naive separation of mean ﬁeld dynamics and many-body dynamics

employed in the current picture. Getting around the ﬁrst issue, to get back to the previous

satisfactory equation of state, we add an auxiliary parameterized contribution to the mean-

ﬁeld energy

Eaux

A

(ρ, δ) = A1(cid:18) ρ
ρ0(cid:19)τ1
exp(cid:26)(cid:18) − b1
+A2(cid:18) ρ
ρ0(cid:19)τ2

ρ0(cid:19)(cid:27)
exp(cid:26)(cid:18) − b2

ρ

ρ

ρ0(cid:19)(cid:27) δ2 ,

(5.13)

and adjust parameters to reproduce the nuclear matter equation of state from the mean-ﬁeld

approximation in Eq. (5.9) as closely as possible. The parameter values are summarized in

Table 5.1 and the results for the equation of state are illustrated in Fig. 5.2. The calculations

are carried out for a uniform system enclosed in a box of size L, with ∆x = 0.57 fm and

periodic boundary condition.

74

A1 [MeV]

39.33

τ1 [MeV]

1.157

b1

1.192

A2 [MeV]

-14.66

τ2 [MeV]

1.533

b2
0.70

Table 5.1: Parameters for the auxiliary ﬁeld in Eq. (5.13).

Figure 5.2: One-dimensional nuclear equation of state in the NGF approach. Crosses rep-
resent the results obtained within pure mean-ﬁeld approximation. Dashed lines represent
a cubic spline ﬁt through the crosses. Open boxes show results obtained within complete
calculations with correlations, where the pure mean-ﬁeld results were matched by adjusting
the parameters in the auxiliary ﬁeld in Eq. (5.13).

To give context to the above results, it should be pointed out that the discretization

interval ∆x introduces a cut-oﬀ for particle momenta in the system, pmax = h/(2∆x). That

cut-oﬀ limits the integration over the momenta in the correlation self-energies, represented

75

00.050.10.150.20.250.3 (fm-3)-20020406080100E/A (MeV)Pure Neutron Matter with MFSymmetric Nuclear Matter with MFMF + Corr (PNM)MF + Corr (SNM)in Fig. 5.1, in addition to what the range parameter η does. Also the size L of the periodic

box determines the spacing for the mesh in momentum, ∆p = h/L, and this impacts both

kinetic and correlation contributions to the energy. With this the auxiliary energy should

have a dependence on ∆x and L, Eaux

A (ρ, δ, ∆x, L). We treat ∆x as a parameter
of our dynamical model, resigned to the fact that to get to analogous dynamic results we

A = Eaux

might need to adjust some other parameters, such as V0 or η in (5.12), when changing ∆x.

Here we choose ∆x = η. As to the dependence on L, we have carried out a number of

calculations at diﬀerent values of L ranging from about 5 fm to 25 fm and diﬀerent nucleon

numbers corresponding to a wide range of density values of interest to us up to 2ρ0. In these

calculations we see variation less than 1 MeV in E/A for uniform matter at ﬁxed ρ, δ and

∆x.

5.3 Preparing ﬁnite correlated nuclear systems

Following the extensions to accommodate the isospin degrees of freedom and two-body resid-

ual interactions and the study of the equation of state for inﬁnite nuclear matter, we turn to

the preparation of ﬁnite nuclear systems in one dimension. Obtaining a good approximation

of the ground state is the starting point for performing simulations of nuclear collisions.

Within the full-ﬂedged NGF approach, we continue to use adiabatic switching tech-

niques [116, 122–125] to evolve the system towards an approximate ground state. At t = 0,

the system is initialized in an uncorrelated harmonic oscillator ground state occupying the

lowest N/2 neutron shells and Z/2 proton shells. The remaining spin degeneracy g = 2

limits us to even-even nuclear systems. With the help of a switching function f (t) varying

smoothly from 0 to 1 over the switching period [0, Ts], the total Hamiltonian H(t) switches

76

smoothly and slowly from the harmonic oscillator ﬁeld HHO to the nuclear ﬁeld Hnucl ac-
cording to H(t) = [1− f (t)] HHO + f (t) Hnucl. The nuclear Hamiltonian Hnucl contains both
the mean-ﬁeld component and the two-body interaction component. Since the second order

diagrams in Fig. 5.1 contain eﬀectively interactions at two times, some extra care needs to
≶ which gets transformed

be taken when adiabatically turning on the residual self-energy Σ

according to:

≶

Σ

≶
(t, t(cid:48)) → f (t) Σ

(t, t(cid:48)) f (t(cid:48)) .

(5.14)

During the adiabatic switching process, the system evolves towards the nuclear ground state.

The quality of the approximation to the ground state depends crucially on the switching

time Ts, with a longer switching time yielding a better approximation. However, the entire

preparation history enters the correlation integrals in Eq. (5.3) in every single time step in

the subsequent evolution. Too long a switching time can stall the numerical calculations

quickly as the burden of history to be remembered grows.

The eﬀects of a cooling friction during switching have been studied previously [123, 124,

126]. Imposing an additional friction-like potential Ufric

n/p in the switching period can to a
certain extent relax the need of an excessively long switching time. The operational form

that we employ for the potential is

fric
n/p(x, t) = f0

U


∆t

∂ρn/p(x, t)

∂t

t
Ts

exp(cid:26)(cid:20) −

t

Ts(cid:21)(cid:27) ,

(5.15)

where f0 is the strength of the friction and ∆t is the time step size in the numerical simula-

tions. While major desired dependencies have been factored out in Eq. (5.15), the strength

parameter f0 remains fragile, with optimal values varying from system to system and from

conditions to conditions. In addition, even for optimal value any improvement in the switch-

77

ing with the inclusion of the cooling friction can be marginal. In practice, extending the

switching time still easily outperforms any delicately designed switching functions or cooling

frictions, when a high-quality ground state is in need. For the time being, we simply im-

plement the cooling friction parametrized in Eq. (5.15) in the code with a relatively small

strength f0 = 4 − 6 fm4/c, without claiming the optimality of the choice.

Figure 5.3: Occupation of the neutron natural orbitals. The open circles represent those
of the correlated ﬁnal state after the adiabatic evolution, and the crosses correspond to the
uncorrelated initial state. The lines serve to guide the eye.

We now study the adiabatic evolution of a symmetric system with N = Z = 6. Initially,

the system is comprised of 3 neutron shells and 3 proton shells in an external harmonic

oscillator potential. The initial wave function can be written as an antisymmetrized product

of the harmonic-oscillator single-particle states, and it is therefore an uncorrelated state. In

78

12345678910Neutron Natural Orbitals00.20.40.60.811.2OccupationCorrelatedUncorrelatedN = Z = 6 (1D)Fig. 5.3, we denote the occupation of the neutron natural orbitals, found by diagonalizing the

one-body neutron density matrix, with red crosses. In the uncorrelated system, the natural

orbitals coincide with the eigenstates of the single-particle harmonic oscillator hamiltonian

hHO, and the lowest three are fully occupied as expected. The case for protons is exactly

the same due to the isospin symmetry present in the system and thus is not shown here.

If the mean ﬁeld were the only interactions switched-on in the system during the adiabatic

evolution, the dynamics would be identical to that in the TDHF theory. Indeed, the NGF

approach with mean-ﬁeld approximation has been studied previously [116, 122]. The ﬁnal

state would simply be uncorrelated, albeit as a product of complicated single-particle states.

In the momentum basis {φk ∼ eikx}, however, even such an uncorrelated state becomes a
complicated superposition of uncorrelated momentum states with all sorts of particle-hole

excitations. While in terms of the complicated states the density matrix can be diagonal,

the density matrix acquires oﬀ-diagonal elements in the momentum representation.

The second-order Born diagram in Fig. 5.1 allows for swapping of orbital pairs in an

uncorrelated state due to the interaction and for mixing of that swapped state with the

original, thus introducing a correlation impacting reinteraction. Depletion in the probability

of staying in the same uncorrelated state is compensated with a probability of populating

other uncorrelated states. In a system with weak nonuniformity, the density matrix becomes

nearly diagonal in the momentum representation. If the system is sparse or mixing is strongly

suppressed by antisymmetrization, impact of a mixing persists over a long time. Eventually,

a kinetic limit is reached where integrations in the residual self-energies can be carried

out to arrive at energy-momentum conservation in a scattering process. The peculiarity of

the conservation in the two-particle scattering in one dimension is that the orbitals only

either retain or interchange momenta. For the species of the same type this leads to the

79

persistence of the uncorrelated states, when composed out of near momentum eigenstates,

and thus no apparent change compared to the mean-ﬁeld dynamics. However, because of the

self-consistency in the diagram in Fig. 5.1, the admixtures of states with swapped orbitals

pile up over time, leading to correlations between more than two orbitals even when residual

interactions act between two particles only. The scattering processes in the limit of energy-

momentum conservation involve then more than two particles and more complicated changes

occur in the uncorrelated states, potentially more similar to analogous processes in three

dimensions than for the isolated two-particle scattering. Also, for a more nonuniform system

the diﬀerent orbitals may share the same momentum, so even if the momenta are retained

or swapped the occupations of uncorrelated states may change. In Fig. 5.3 we observe a

signiﬁcant depletion of the occupation for the low-lying states and increased occupation for

the higher, in spite of the one-dimensional peculiarity of the on-shell two-body scattering

and of the two-body scattering rate being coarsely set to be the same as for the three

dimensions. On the same note, even in three dimensions, when insisting on a quasiparticle

picture for a ﬁnite system with energy conservation for single-particle energies, the departure

of occupations from 1 and 0 would never occur.

If a particle is removed from or placed in a natural orbital from a stationary system with

weak correlations, the system largely remains in a stationary state. This can be quantiﬁed

in terms of a spectral function with a unique single-particle energy, equal to the diﬀerence in

the energies for the two stationary states. Inspired by Fermi liquid theory, this independent

quasiparticle picture is pushed to the extreme in the semiclassical transport theory for nuclear

reactions where speciﬁc energies are attributed to quasiparticles with deﬁnite momenta [66].

The spectral function naturally relates to NGF. One characteristic of increasing correlation

strength is that identiﬁcation of natural orbitals becomes a challenge. They may be always

80

obtained by diagonalizing the density matrix but, as far as evolution of NGF in the relative

time is concerned, the orbitals can be discussed only in a perturbative context. The spectral

function is not diagonal in any orbitals anymore, but can be near-diagonal and/or diagonal

elements of the spectral function can be considered.

Given the above, in the correlated system we look at diagonal elements of the spectral

function in momentum representation. For particular species with particular spin projection,

denoted with ν (also speciﬁcally used for neutrons in the shell-model notation), the spectral

function Sν(p, E) is related to the retarded Green’s function G+

ν (p, E) with

Sν(p, E) = −

1
π

Im G+

ν (p, E) .

(5.16)

The retarded Green’s function G+

ν (p, E) admits a Lehmann representation,

G+

ν (p, E) =(cid:88)p

|(cid:104)Ψ
E − (E

(N +1)
p

(N +1)
p

|φ†ν(p)|Ψ
− E
(N−1)
h

(N )
0
(N )
0

(cid:105)|2
) + iη
|φν(p)|Ψ
− E

(N−1)
h

(N )
0
(N )
0

(cid:105)|2
) + iη

+(cid:88)h

|(cid:104)Ψ
E − (E

,

(5.17)

where η → 0+, φ and φ† are annihilation and creation operators in momentum representation
(N−1)
(N +1)
respectively, Ψ
are eigenstates of the total Hamiltonian with N+1 neutrons and Ψ
p
h
are eigenstates with N−1 neutrons. The spectral function Sν(p, E) picks out the poles
of retarded Green’s function, which are the excitation energies of adding or removing a

neutron of momentum p with respect to the N-neutron ground state Ψ0. The strength of

the excitation is regulated by the overlap function in the numerators in the retarded Green’s

function. On account of the completeness of the sets of states and commutation relation

for the operators, the spectral function integrates over energy to 1. This allows one to

81

interpret Sν(p, E) as the probability density for a particle ν with momentum p to have an

energy E [127].

In a noninteracting Fermi gas the spectral function peaks at the kinetic energy of a

particle with momentum p, i.e. S(p, E) = δ[E − p2(cid:14)(2m)]. For a self-bound slab in a box, in

the semiclassical limit behind the transport for reactions, two such peaks are expected, one

at lower energy for slab interior and one at higher for the exterior, with relative intensity

of the peaks representing the relative share of the size of the box for the slab interior and

exterior. In the Fermi liquid theory the spectral function of a correlated uniform system
(E−Ep)2+[1/(2τp)]2 + S(cid:48)(p, E), where Ep
is often sought in the Lorentzian form: S(p, E) ≈
is the energy of the quasi-particle with momentum p, τp is the lifetime, Zp is quasiparticle

Zp/(2τp)

strength and S(cid:48) is contribution of incoherent background that ensures that the sum rule is

satisﬁed [127].

In Fig. 5.4, we show the spectral functions for neutrons in a symmetric system with

N = Z = 6 in an L = 15 fm box, at six of the discrete momenta, in the mean-ﬁeld calculation

and in the full calculation with two-body residual interaction. The system were evolved

up to 300 fm/c, and we Fourier transformed the retarded Green’s functions of the form
G+(p, t1; p, t2) with t1 ≥ t2 to get the retarded Green’s functions G+(p, E) in momentum
and energy. The spectral functions Sν were computed according to Eq. 5.16. For reference,

we provide in the ﬁgure the location where the kinetic energy for a given momentum is.

Multiple spikes are present in both calculations. For three of the lower momenta, strength

is present at about 55 MeV down from the free-space kinetic energy, about the depth of

the potential well in nuclear systems. The spikes in the mean-ﬁeld case are sharper, and,

for the three lowest momenta, a few of the peaks in the two diﬀerent cases overlap almost

perfectly, suggesting that the mean ﬁeld alone is already responsible for a signiﬁcant part

82

of the fragmentation of natural orbitals in momentum representation. Notably in the limit

of perfectly static inﬁnite-time evolution, the mean-ﬁeld peaks would turn into δ-functions,

meaning that any width observed for these peaks in Fig. 5.4 is due to our time domain

methodology.

Inclusion of residual interactions smoothes out and spreads the peaks and

moves some well beyond any resolution impact of our methodology.

Figure 5.4: Neutron spectral functions vs energy for an N = Z = 6 self-bound system in an
L = 15 fm box at discretized momenta. The red dashed-dot line indicates the kinetic energy
of a free neutron with momentum p. See the text for discussion.

The structures observed in Fig. 5.4 are outside of the reach of semiclassical transport

theory. Telling is the fact that to resolve the structures in the spectral function reasonably

well, we had to evolve the system for over 300 fm/c. Key stages of reactions to which

transport is applied commonly last a twentieth of that time [66, 115], meaning that energy

structure will generally be poorly resolved and conservation in binary collisions will not

hold [128].

83

-100-5005010000.020.040.060.08Neutron S(p,E) (MeV -1)p = 0 MeV/cEpkin-100-5005010000.020.040.060.08p = 82.7 MeV/cEpkin-100-5005010000.020.040.060.08p = 165.3 MeV/cEpkinCorrelatedUncorrelated (MF)-5005010015000.040.080.120.16p = 330.6 MeV/cEpkin50100150200250E (MeV)00.040.080.120.16p = 578.6 MeV/cEpkin35040045050055000.040.080.120.16p = 909.2 MeV/cEpkinFigure 5.5: Momentum state occupation for neutrons in the same N = Z = 6 ground-state
system as in the preceding two ﬁgures and for a Fermi gas. The dashed, dash-dotted and
solid lines represent the Fermi gas at temperature T = 3.6 MeV, the mean ﬁeld system and
the correlated system, respectively.

In Fig. 5.5, we show the momentum distribution for neutrons in the same N = Z = 6 sys-

tem we have been discussing and for a noninteracting Fermi gas at temperature T = 3.6 MeV.

The high-momentum tails seen in the uncorrelated mean-ﬁeld case and the case with short-

range two-body interaction lack the features of a free Fermi gas at ﬁnite temperature, sug-

gesting again that the underlying picture cannot be entirely captured in the Fermi liquid

theory. In particular, the high-momentum (> p1D

F ≈ 195 MeV/c) tail accounts for about 25%
of the neutrons in the correlated system, comparable to what is reported to be around 20%

in nuclei in nature [129–131], while the corresponding proportion in the uncorrelated system

is just about 10%. A signiﬁcant number of the high-momentum nucleons are present due to

84

0100200300400500600700800900p (MeV/c)10-610-510-410-310-210-1100f(p)Fermi Gas (T = 3.6 MeV)Uncorrelated (MF)Correlatedthe short-range two-body interaction.

5.4 Exploring dissipation with isovector oscillation in a

ﬁnite system

To describe heavy-ion collisions at the intermediate energies ranging from tens of to hun-

dreds of MeV/u, it is important to include dissipation of the relative energy in the center

of mass frame in the dynamics, which is lacking in the conventional TDHF theory. How-

ever, as discussed in the preceding section, eﬀects of the short-range two-body interactions

included in the correlation integrals in the NGF approach can be interpreted as scattering

and rearranging population of orbitals. In fact, these integrals turn into collision integrals

in the kinetic limit for the theory.

To explore the dissipative eﬀects introduced by the short-range two-body interactions

in the NGF approach, we simulate the isovector-like oscillation of neutrons and protons

in a ﬁnite system. The simulations are performed in two diﬀerent settings, one with only

mean ﬁeld dynamics, which is equivalent to the case in TDHF theory, and the other with

both mean ﬁeld dynamics and residual two-body interactions. For each setting, we prepare

a symmetric nuclear system with N = Z = 4 in its ground state at the center of a box,

employing the adiabatic switching technique discussed before. The time is reset to t = 0

fm/c after completion of the adiabatic evolution. At t = 0 fm/c, we then perturb the

neutron slab and the proton slab, which have been overlapping with each other perfectly.

Speciﬁcally, we apply opposite jolts to the neutron and proton slabs, that take the form of

85

Galilean boosts to the same-time lesser Green’s functions G<

n/p at t = 0:

≶, after
n/p

G

(x, 0; x(cid:48), 0) = exp(cid:26)(cid:16) ± i

≶, before
n/p

×G

P x

 (cid:17)(cid:27)
(x, 0; x(cid:48), 0) exp(cid:26)(cid:16) ∓ i

P x(cid:48)

 (cid:17)(cid:27),

(5.18)

where P is the small momentum boost of the order of 55 MeV/c, equivalent to 1.6 MeV in

kinetic energy per nucleon. After the jolts, the neutron and proton center of mass start

to move away from one another, but as particles reach the edge of their optical potentials,

they retreat and the centers of mass move back. Fig. 5.6 shows the evolution of the density

proﬁles for neutrons and protons for the setting with mean ﬁeld only. The position of the

center of mass for the neutron slab ¯X can be seen, for both scenarios, in the top panel of

Fig. 5.7. When the residual interaction is present, the isovector oscillation is much more

strongly damped.

86

Figure 5.6: Isovector-like oscillation of the neutron slab and the proton slab with mean-
ﬁeld approximation. The red line represents the neutron density proﬁle and the blue line
represents the proton density proﬁle.

On top of the isovector oscillation in Fig. 5.6, one may observe there deformations of the

neutron slab and the proton slab with respect to their own center of mass. Since the slabs

are not rigid, it is conceivable that one induces other forms of oscillation. The additional

changes may progress at the same or multiple of the isovector frequency, or at a completely

diﬀerent frequency. For more insight we evaluate two more moments for a neutron slab and

plot them as a function of time in the lower two panels of Fig. 5.7. Speciﬁcally, the center
1
2 ,

panel displays the second moment about the center of mass of the neutron slab, (cid:104)(X− ¯X)2(cid:105)
describing the breathing mode for the neutron slab and also for the whole system after the

isovector mode dies down. Finally, the bottom panel displays the third moment about the

87

00.020.040.060.080.1 (fm-3)-10-505x (fm)00.020.040.060.080.1-10-50510ProtonNeutront = 59.8 fm/ct = 46.2 fm/ct = 0 fm/ct = 14.4 fm/ccenter of mass of the neutron slab, (cid:104)(X − ¯X)3(cid:105)
the center of mass.

1
3 , representing skeweness of the slab around

Figure 5.7: Evolution of the lowest three moments of the neutron slab. The red lines and
the black lines correspond to calculations with mean ﬁeld only and the full calculations,
respectively.

The higher moments, quantifying behavior of the farther-out features of the density

distribution, reach higher values in Fig. 5.7 than does the center of mass position. No

matter what moment is considered the damping of its oscillations is much stronger for the

88

-0.500.51.81.92050100150200250-202dynamics with correlations than without. In fact within the considered time interval hardly

any damping is observed for the pure mean-ﬁeld dynamics. While the oscillations in the

center of mass are fairly sinusoidal, indicating the dominance of a single frequency, this is

less so for the higher moments. The oscillations in the third moment are in phase with

those in the ﬁrst, but are deﬁnitely anharmonic, suggesting that initialization of the Green’s

functions, with the same boost parameter at low densities as at high, might be too naive.

The second moments exhibit diﬀerent periodicity than the odd moments, suggesting that an

isoscalar breathing mode got excited and progresses at a diﬀerent frequency than isovector.

Both the third and second moment damp away at a slower rate than the ﬁrst moment, when

the correlations are present.

In the case of the third moment, the diﬀerence in damping

compared to the ﬁrst moment is subtle and can be due to the fact that the correlations

weaken at the low densities that get tested by the high moment. In the case of the second

moment the damping is weak at longer time and we suspect that this is due to the peculiar

nature of two-body scattering in the semiclassical limit in one dimension, where the orbitals

either persist or interchange. The interchange will shuﬄe the momentum between neutron

and proton slabs, damping the isovector mode, but do nothing to the isoscalar mode once

the isovector mode dies out. The glitch in the third moment occurring at 25 − 30 fm/c is
likely due to ripples in the density tails, produced by the boost, coming from both sides and

colliding at the boundary of the periodic box. Some early time transient behavior in the

correlated dynamics is also due to the fact that the jolts are taken to be of a short duration

compared to the relaxation of correlations.

89

5.5 Galilean Covariance

The Kadanoﬀ-Baym equations transform covariantly under Galilean transformation of the

reference frame. That means that Galilean transformed solutions from one frame should

remain solutions of the equations in another frame. Obviously there is no guarantee that

this remains true in a numerical solution. Discretization in space introduces a cut-oﬀ in

momentum space that does not transform covariantly. Practical violations of the covariance

may impair studies of moving matter whether in collective oscillations or in collisions. In the

context of the latter, we need to prepare slabs in their ground state, boost them and direct

against each other. Due to violation of covariance, the boost can potentially excite the slab

and lead to its break-up ahead of any collision.

In this section, we explore impact of a boost on a slab in practice. In the calculation we

use a box of size L = 15 fm with periodic boundary conditions. We prepare a symmetric

system with N = Z = 2 at the center of the box as usual. We then boost the slab to provide

it with a momentum per nucleon P . When working with a periodic box, the momenta in

representing the Green’s functions are discretized with the spacing of ∆p = h/L. We simplify
the boost by taking P = 2∆p, or wavevector k = 2(2π/L) (cid:39) 0.838 fm−1, corresponding to
boost velocity V = P/m = 0.176 c. With this the Green’s functions can be shifted over

momentum mesh without interpolations.

The correlation integrals in the Kadanoﬀ-Baym equation involve the integration of the

entire history of the system, rendering the dynamics non-Markovian, in sharp contrast to

both semi-classical transport theories and the TDHF theory. If the boost were applied to the

end of the adiabatic evolution, we would be boosting the single-particle characteristics, but

not two-particle correlations. We can aﬀord it to a degree when exciting isovector oscillations,

90

risking some heating of the system in addition to that occurring anyway due to dissipation

of the collective oscillation. The local Galilean boost operator T (x, t; V ), that transforms a
Green’s function from a frame of a still system with positions in terms of coordinate x to

one in terms of x(cid:48) = x + V t, according to

G

≶
≶
moving(x(cid:48)1, t1; x(cid:48)2, t2) = T (x1, t1; V ) G
still(x1, t1; x2, t2)

×T ∗(x2, t2; V ) ,

(5.19)

is

T (x, t; V ) = exp(cid:26)(cid:104) iP

 (x + V t/2)(cid:105)(cid:27) .

(5.20)

In momentum representation this yields

≶
moving(p1, t1; p2, t2) = e−
G

iV t1 (p1−P/2)

≶
still(p1 − P, t1; p2 − P, t2)
G
iV t2 (p2−P/2)

×e

.

(5.21)

When starting investigation of a boosted slab at t0 = 0, the transformation is applied to

functions where one of the arguments coincides with t0 and the other either coincides with

or precedes t0.

After applying the boost transformation in momentum, we evolve the system according

to the Kadanoﬀ-Baym equations in the standard manner. The evolution of the density proﬁle

for a boosted slab is illustrated in Fig. 5.8. Speciﬁcally in panel (a) we show the initial proﬁle

with a dash-dot line. In the subsequent panels (b)-(d), we subsequently advance the evolution

by 30 fm/c in each, displaying the latest proﬁle with a solid line and intermediate proﬁles, at

15 fm/c intervals, with dashed lines. It is apparent that the slab moves at a uniform pace.

91

Some small variations in the proﬁle are seen, potentially related to spatial discretization, but

these are not signiﬁcant enough to be of concern in exploring slab collisions. With this, we

demonstrate covariance of the NGF approach is adequate enough for future studies of the

collisions.

Figure 5.8: Evolution of density for a correlated N = Z = 2 slab boosted with velocity
V = 0.176 c. The initial density is illustrated with a dash-dotted line in every panel. The
panels (b)-(d) show evolution advanced subsequently over 30 fm/c, with the latest proﬁle
illustrated with a solid line.
Intermediate proﬁles at every 15 fm/c are illustrated with
dashed lines.

5.6 Summary and future prospects

In this chapter, we have discussed the application of the NGF theory to the description of the

dynamics of one-dimensional correlated nuclear systems in diﬀerent scenarios. The prepa-

ration of ﬁnite nuclear systems is relevant for the initialization of projectile and target in

92

00.050.10.150.20.250 fm/c0 fm/c15 fm/c30 fm/c-6-4-20246X (fm)00.050.10.150.2 (fm-3)0 fm/c15 fm/c30 fm/c45 fm/c60 fm/c-6-4-202460 fm/c15 fm/c30 fm/c45 fm/c60 fm/c75 fm/c90 fm/c(d)(c)(a)(b)simulating nuclear collisions. The dissipation arising from the short-range two-body interac-

tions, demonstrated in the isovector-like slab oscillation, can give rise to the stopping eﬀects

in collisions. One should bear in mind that the dissipation in one dimension simulated in

the current numerical model may not have captured the entirety of the actual dissipation in

three dimensions. In a three-dimensional description, relative longitudinal momentum can

be eﬀectively redirected into the transverse directions, presumably enhancing the dissipation.

In the current one-dimensional description, dissipation eﬀects are still largely conﬁned to the

internal excitation of the systems in terms of exchanging momentum. The translation of a

stable slab partially mimics the situation of two nuclei approaching each other before mak-

ing contact, but we do acknowledge that boosting two stable slabs in the opposite directions

requires more careful design and is thus out of the scope of this work.

We shall also brieﬂy touch on the limitations and possible future extensions of the cur-

rent NGF model. The tremendous amount of CPU hours and memories needed to solve the

Kadanoﬀ-Baym equation numerically renders direct simulations in higher dimension com-

putationally infeasible. In fact, even in the current model, we needed to look for techniques

to reduce the computational cost and the calculations have only been possible on high-

performance computer clusters. We have used OpenMP [132] to parallelize the code so that

independent loop iterations can be carried out simultaneously on multiple cores. We also

made the observation that, at any given moment t, in propagating the Green’s functions
G<(·, t(cid:48),·, t) and G>(·, t,·, t(cid:48)) with t ≥ t(cid:48), only the history correlated with the current is
relevant, i.e., quantities with time arguments (t, t(cid:48)) or (t(cid:48), t). This allows for the use of arrays

with just one dimension in time for t(cid:48), which can be overwritten during the propagation of

t, in sharp contrast to the naive way of storing quantities in both of the temporal dimen-

sions. The optimization made so far makes no compromise on the exactness of the numerical

93

solutions in the stated framework. Further reduction in time or memory is possible under

the assumption that elements far oﬀ the diagonal are negligible. In fact, in Ref. [116] it was

demonstrated that far oﬀ-diagonal spatial elements of Green’s functions may be dropped

with impunity and one may hope to be able to introduce a ﬁnite memory time limiting

extent in time. These could reduce the time and memory cost of computations.

Regarding extension of the scope, the current separate parametrizations of mean ﬁeld

and residual interactions represents some limitation. Progress towards consistency between

the interactions could be reached by solving the T -matrix equation [119] with an interaction

more representative of a bare interaction. With a separable approximation to the latter, the

computational eﬀort would not be much more serious than the current. At low energies,

excitations of chaotic collective modes can compete with nucleon-nucleon collisions in equili-

brating larger colliding systems. Accounting for such excitations can be achieved with a ring

approximation to the self energy, again requiring quite similar computational eﬀort to the

current. Beyond the short-range correlations, extensions can be conjectured to the current

NGF model to treat pairing correlations upon diﬀerentiation of the spin degrees of freedom.

Pairing eﬀects have been incorporated into collisions in a three-dimensional TDHF code

through the initialization of Bardeen-Cooper-Schrieﬀer (BCS) ground states [105], with the

caveat that the pairing amplitudes are kept ﬁxed during the evolution. A more proper way

to treat pairing correlations in time-dependent processes is to solve the TDHF-Bogoliubov

(TDHFB) equation [27, 103, 104]. It is unclear how pairing can enter the current picture of

the NGF theory. Still, following the footprints of the TDHFB theory, one may try including

a time-diagonal pairing-Green’s functions of the form (cid:104)φ(x1, t) φ(x2, t)(cid:105), analogous to the
abnormal density or the pairing tensor κ(t) in the TDHFB theory. A pairing term ∆ shall

enter the Hamiltonian when the Kadanoﬀ-Baym equations for the normal Green’s functions

94

are solved. An additional equation for the time-diagonal pairing-Green’s functions void of

the short-range interactions, mimicking that for the pairing tensor [104], also needs to be

solved simultaneously. Within such a naive conjecture, violation of total particle number

and total energy might arise though.

We conclude that this work extends the application of the NFG theory in nuclear physics

beyond the mean ﬁeld approximation and lays the ground for the future study of nuclear

collisions, taking us one step closer to a fully quantum-mechanical and realistic transport

model. A multitude of extensions and optimization can also be considered and implemented

for the improvement of this model.

95

Chapter 6

Conclusion

In this thesis, we present two time-dependent transport models for HIC, which represent the

two major directions for future developments of transport theory. The one-body Langevin

transport model incorporates ﬂuctuations to take into account the eﬀects of many-body

correlations and attempts to address fragmentation in HIC at intermediate energies; the NGF

transport model treats nuclear collisions in a fully quantum-mechanical and self-consistent

manner, capable of describing both the mean-ﬁeld dynamics and the short-range two-body

correlations.

In Chapter 2, we described how we had arrived at a purely one-body transport model

based on Langevin dynamics starting from the conventional BUU equation. The formu-

lation of the new model was discussed at length and the technicalities in the numerical

implementation, including a self-consistent method to solve the Thomas-Fermi equation for

the initialization of stable nuclei, were carefully documented.

In Chapter 3, we presented a series of benchmark results to demonstrate the reliability

and compatibility of the one-body Langevin transport model. We ﬁrst simulated the Au

+ Au collisions at two incident energies, 100 MeV/nucl and 400 MeV/nucl, at an impact

parameter of 7 fm. We compared the rapidity distributions dN/dy, the average in-plane
ﬂow (cid:104)px/A(cid:105) and the ﬂow parameter d(cid:104)px/A(cid:105)/dyred and found good agreement with other
existing transport codes at both incident energies under the same controlled conditions. We

96

also investigated the nearly central collision Sn + Sn at 50 MeV/nucl, where our model was

again shown to be robust in describing the fragmentation dynamics involved. Results were

compared with two other mature models, SMF and AMD. Despite the drastic diﬀerences in

their treatment of stochastic dynamics, all models yield similar distributions for the primary

IMF multiplicity.

In Chapter 4, we investigated the so-called “hierarchy eﬀect” in the correlation between

charge and velocity of IMFs reported in experiment. The collision of focus was Ta+Au at

39.6 MeV/nucl. By analyzing the properties of the primary IMFs generated in the simulation

with the Brownian motion model, we conﬁrmed that the IMF with the greatest charge

tends to have the greatest velocity component and be forward-peaked, and the existence

of a weak positive correlation between charge and velocity. These observation contrast the

idea of IMF being emitted from an equilibrated source, and can only be interpreted with

a dynamical transport model, such as the one presented in this thesis. We also simulated

U+C at 24 MeV/nucl and could not predict any relevant fragmentation phenomenon, due

to the model’s inability of describing statistical decay. This, however, does not contradict

the lack of “hierarch eﬀect” reported for the same system in experiment.

In Chapter 5, we ﬁrst presented a quantum-mechanical transport model based on NGF,

and, in particular, showed how short-range two-body correlations were naturally included in

the full-ﬂedged Kadanoﬀ-Baym equation. The EoS E/A(ρ) of correlated nuclear matter in

one dimension was calculated. We demonstrated the correlation eﬀects by comparing natural

orbital occupancy n, spectral functions S(p, E) as well as the momentum distribution f (p) of

nuclear systems in 1D, with or without two-body correlation. The dissipative eﬀects brought

by two-body correlation was evidenced by simulation of the damping of spatial multipole

moments (cid:104)(X− ¯X)n(cid:105)

1
n of 1D nuclear systems. In the end, we also veriﬁed that the simulation

97

preserved the Galilean covariance, which paved the way for future study of collision of 1D

slabs.

Lastly, the work in this thesis represents a narrow niche of the development of transport

theories in the broad study of HIC physics, and certain aspects of the transport models

discussed here may come across as ad hoc and/or primitive. It is hoped that this record of

work will shed light on the issues of the incorporation of ﬂuctuations and quantum eﬀects

into transport models, and be a potential source of inspiration for any future development.

98

APPENDICES

99

APPENDIX A: Fokker-Planck equation

and Langevin equation

In this appendix, we provide a summary of the derivation of the Fokker-Planck equation

detailed in Ref. [80] and discuss the corresponding Langevin equations of diﬀerent forms.

Consider the following Boltzmann collision integral for arbitrary statistics,

Icoll =

g

h3(cid:90) d3pb(cid:90) dΩ

dσab
dΩ

vab[ ˜fa ˜fbfa(cid:48)fb(cid:48) − fafb ˜fa(cid:48) ˜fb(cid:48)].

(A.1)

fx is to be regarded as the shorthand notation for the one-body phase distribution function
f (px). ˜f = 1 + λf with λ = −1, 0, 1 corresponding to Fermi-Dirac, Boltzmann, and Bose-
Einstein statistics. g is the degeneracy factor.

For a pair of particles with momenta pa and pb, vab is the relative velocity, and dσab/dΩ

is the diﬀerential cross section, where the scattering angle Ω is deﬁned in the center-of-mass

frame of the pair. Final momenta are denoted by primed subscripts.

In the center-of-mass frame, the initial state of the colliding pair is characterized by the

relative momentum qab = pa − pb and the total momentum Pab = pa + pb. When only
elastic collisions are under consideration, Pab remains constant and qa(cid:48)b(cid:48) = qab. It follows
that the ﬁnal relative momentum qa(cid:48)b(cid:48) is completely determined by the scattering angle
Ω = (θ, φ).

For a φ-independent and forward-peaked cross section, we can ﬁrst make the following

expansion over the polar angle θ at a ﬁxed azimuthal angle φ, and then truncate it up to

100

the leading order term upon integration over φ,

˜fa ˜fbfa(cid:48)fb(cid:48)− fafb ˜fa(cid:48) ˜fb(cid:48) = ˜fa ˜fb(cid:104) ∞(cid:80)n=0

θn
n!

−fafb(cid:104) ∞(cid:80)n=0

θn
n!

∂n

∂θn (fafb)(cid:105)(cid:12)(cid:12)(cid:12)θ=0

.

∂n

∂θn ( ˜fa ˜fb)(cid:105)(cid:12)(cid:12)(cid:12)θ=0

The zero-order term vanishes conveniently. Since

qa(cid:48)b(cid:48) = qab(sinθ cosφ, sinθ sinφ, cosθ),

(A.2)

(A.3)

(A.4)

the usual chain rule gives

∂
∂θ

=(cid:88)i

∂qi

a(cid:48)b(cid:48)

∂θ (cid:12)(cid:12)(cid:12)θ=0

,

∂
∂qi

a(cid:48)b(cid:48)

where ∂qa(cid:48)b(cid:48)/∂θ|θ=0 = qab(cosφ, sinφ, 0). The integration of ∂qa(cid:48)b(cid:48)/∂θ|θ=0 over φ is zero, so
the leading order term in Eq. (A.2) is of second order.

It is easy to show that

(cid:90) 2π

0

dφ

∂2

∂θ2(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)θ=0

= π(cid:32) − 2(cid:88)i

qi ∂

∂qi + q2(cid:88)ij

∆ij

∂2

∂qi∂qj(cid:33),

(A.5)

where the subscript ab for q is suppressed for brevity and ∆ij = δij − qiqj/q2 is a projection
operator onto the plane perpendicular to qab.

Upon inserting the truncated expansion integrated over φ into the Boltzmann collision

integral Eq.

(A.1) and noting that ∂/∂qi

ab = 1

2(∂/∂pi

a − ∂/∂pi

b), one will arrive at the

following Fokker-Planck equation,

∂fa
∂t

= −∇pa ·(cid:104)1

2(cid:16) ˜Ra + ˜faRa(cid:17)fa(cid:105) + ∇2

pa(cid:16)Dafa(cid:17),

(A.6)

101

where

and

Ri
a = −
˜Ri
a = −
g
ij
a =

D

g

g

ab,

h3(cid:90) d3 pb fb Fab qi
h3(cid:90) d3pb fb ˜fb Fab qi
4h3(cid:90) d3pb fb ˜fbFab(q2

ab,
abδij− qi
abq

(A.7)

(A.8)

(A.9)

j
ab),

Fab = (πvab/2)(cid:90) 1

0

θ2 (dσab/dΩ) d cos θ.

(A.10)

By assuming that the evolution of particle momentum p is a Gaussian random process

and that the evolution of the momentum distribution f (p, t) follows the Fokker-Planck Eq.

(A.6), one can write down the corresponding single-particle Langevin equation of the Itô

form [133, 134],

dp =

1

2(cid:16) ˜R + ˜f R(cid:17)dt + σdBt,

where σ is a 3×3 positive deﬁnite matrix such that
2(cid:88)k

Dij =

1

σikσjk.

Bt denotes a Guassian random process with properties

(cid:104)dBt(cid:105) = 0,
tdB

j

t (cid:105) = dt δij.

(cid:104)dBi

(A.11)

(A.12)

(A.13)

(A.14)

Here and in what follows, some subscripts may be dropped whenever no confusion arises.

102

When the cross section dσab/dΩ is taken to be independent of qab, one can verify that

∇pa · D = 1

2

˜R. Hence, the Fokker-Planck Eq. (A.6) has then an equivalent simpliﬁed form,

∂f
∂t

= −∇p ·(cid:16)1

2

˜f R f(cid:17) + ∇p · (D∇pf ).

(A.15)

This form of the Fokker-Planck Eq. (A.15) has the corresponding Langevin equation of the

Stranovich form [133, 134],

dp =

1
2

˜f R dt + σ ◦ dBt.

(A.16)

Note that the two forms of Langevin equations are equivalent. For the purposes of numerical

integrations, in the Itô form, successive increments are evaluated at the beginnings of each

time step, while they are evaluated at the mid-points of each time step in the Stranovich

form. Readers may refer to Ref. [133,134] for a more detailed discussion of the two forms of

Langevin equations.

103

APPENDIX B: Generalized Einstein

relation in thermal equilibrium

Consider a system of particles of arbitrary statistics characterized by a momentum distribu-

tion function f (p, t), whose evolution follows the aforementioned Fokker-Planck equation,

∂f
∂t

= −∇p ·(cid:104)1

2

(1 + λf )Rf(cid:105) + ∇p · (D∇pf ),

(A.17)

where λ = −1, 0, +1 correspond to fermions, classical particles and bosons, respectively. The
thermal equilibrium momentum distribution function feq(p; T ) at constant temperature T

is stationary, and hence should be solutions to Eq. (A.17) in the absence of external ﬁelds.

The two terms on the right hand side of Eq. (A.17) can be identiﬁed with divergences of

the dissipative momentum current Jdiss = Ffricf = 1
2(1+λf )Rf and the diﬀusive momentum
current Jdiﬀ = −D∇f. In thermal equilibrium, at any arbitrary location p in momentum
space, the net momentum current must vanish: J

(eq)
diﬀ (p) = 0, i.e.,

(eq)
diss (p) + J

(eq)
diss + J

(eq)
diﬀ =

J

=

1
2
1
2

(1 + λfeq)Reqfeq − Deq∇feq
(1 + λfeq)Reqfeq

+Deqfeq(1 + λfeq)

p

mkBT

= 0,

(A.18)

104

which yields the generalized Einstein relation in thermal equilibrium,

Deq p
mkBT

1
2

= −

Req.

(A.19)

Note that the drift coeﬃcient Req and the diﬀusion tensor Deq are evaluated with the same

equilibrium distribution feq at constant temperature T , and that this generalized Einstein

relation holds true for all types of statistics. The relation can also be broken down into a

component-wise form,

for i = 1, 2, and 3.

(cid:80)j D
ij
eqpj
−Rieq

=

mkBT

2

,

(A.20)

Direct numerical integrations of the coeﬃcients by deﬁnitions in Eq. (A.7) – (A.9) have

also been performed to conﬁrm the generalized Einstein relation in the thermal equilibrium

limit.

105

BIBLIOGRAPHY

106

BIBLIOGRAPHY

[1] J. Hong. Bulk nuclear properties from dynamical description of heavy-ion collisions.

Michigan State University, 2016.

[2] M. Colonna. Collision dynamics at medium and relativistic energies. Prog. Part. Nucl.

Phys., 2020.

[3] J. Xu, L.-W. Chen, M. B. Tsang, et al. Understanding transport simulations of heavy-
ion collisions at 100A and 400A Mev: Comparison of heavy-ion transport codes under
controlled conditions. Phys. Rev. C, 93:044609, Apr 2016.

[4] M. Colonna, A. Ono, and J. Rizzo. Fragmentation paths in dynamical models. Phys.

Rev. C, 82:054613, Nov 2010.

[5] J. Colin, D. Cussol, J. Normand, N. Bellaize, R. Bougault, et al. Dynamical eﬀects in

multifragmentation at intermediate energies. Phys. Rev. C, 67:064603, Jun 2003.

[6] A. W. Steiner, M. Prakash, J. M. Lattimer, and P. J. Ellis.

nuclei and neutron stars. Phys. Rep., 411(6):325–375, 2005.

Isospin asymmetry in

[7] C. Y. Tsang, M. B. Tsang, P. Danielewicz, F. J. Fattoyev, and W. G. Lynch. Insights

on Skyrme parameters from GW170817. Phys. Lett. B, 796:1–5, 2019.

[8] A. D. Panagiotou, M. W. Curtin, H. Toki, D. K. Scott, and P. J. Siemens. Experi-
mental evidence for a liquid-gas phase transition in nuclear systems. Phys. Rev. Lett.,
52(7):496, 1984.

[9] J. Pochodzalla, T. Möhlenkamp, T Rubehn, A. Schüttauf, A. Wörner, E. Zude,
M. Begemann-Blaich, T. Blaich, H. Emling, A. Ferrero, et al. Probing the nuclear
liquid-gas phase transition. Phys. Rev. Lett., 75(6):1040, 1995.

[10] S. Das Gupta, A. Z. Mekjian, and M. B. Tsang. Liquid-gas phase transition in nuclear

multifragmentation. In Adv. Phys., pages 89–166. Springer, 2001.

[11] Ph. Chomaz, M. Colonna, and J. Randrup. Nuclear spinodal fragmentation. Phys.

Rep., 389(5-6):263–440, 2004.

107

[12] P. Danielewicz, H. Lin, J. R. Stone, and Y. Iwata. Spinodal instability at the onset of

collective expansion in nuclear collisions. arXiv preprint arXiv:1910.10500, 2019.

[13] G. D. Westfall, J. Gosset, P. J. Johansen, A. M. Poskanzer, W. G. Meyer, H. H.
Gutbrod, A. Sandoval, and R. Stock. Nuclear ﬁreball model for proton inclusive spectra
from relativistic heavy-ion collisions. Phys. Rev. Lett., 37(18):1202, 1976.

[14] J. Cugnon, T. Mizutani, and J. Vademeulen. Equilibrium in relativistic nuclear colli-

sions. a Monte-Carlo calculation. Nucl. Phys. A, 352:505–534, 1981.

[15] G.F. Bertsch and S. Das Gupta. A guide to microscopic models for intermediate energy

heavy-ion collisions. Phys. Rep., 160(4):189 – 233, 1988.

[16] R. Trockel, K. D. Hildenbrand, U. Lynen, W. Müller, H. J. Rabe, H. Sann, H. Stelzer,
W. Trautmann, R. Wada, E. Eckert, et al. Onset of multifragment emission in heavy-
ion collisions. Phys. Rev. C, 39(2):729, 1989.

[17] Y. D. Kim, M. B. Tsang, C. K. Gelbke, W. G. Lynch, N. Carlin, Z. Chen, R. Fox,
W. G. Gong, T. Murakami, T. K. Nayak, et al. Multifragment emission observed for
the reaction 36Ar+ 238U at E/A= 35 Mev. Phys. Rev. Lett., 63(5):494, 1989.

[18] E. Suraud, S. Ayik, J. Stryjewski, and M. Belkacem. The Boltzmann-Langevin equation
and its application to intermediate mass fragment production. Nucl. Phys. A, 519(1-
2):171–182, 1990.

[19] Ph. Chomaz, G. F. Burgio, and J. Randrup. Inclusion of ﬂuctuations in nuclear dy-

namics. Phys. Lett. B, 254(3-4):340–346, 1991.

[20] G. F. Burgio, Ph. Chomaz, and J Randrup. Fluctuations in nuclear dynamics from

transport theory to dynamical simulation. Nucl. Phys. A, 529(1):157–189, 1991.

[21] F. Chapelle, G. F. Burgio, Ph. Chomaz, and J. Randrup. Fluctuations in nuclear
dynamics: Comparison of diﬀerent methods. Nucl. Phys. A, 540(1-2):227–260, 1992.

[22] W. Bauer, G. F. Bertsch, and S. Das Gupta. Fluctuations and clustering in heavy-ion

collisions. Phys. Rev. Lett., 58(9):863, 1987.

[23] M. Colonna, M. Di Toro, A. Guarnera, et al. Fluctuations and dynamical instabilities

in heavy-ion reactions. Nucl. Phys. A, 642(3):449 – 460, 1998.

108

[24] P. Napolitani and M. Colonna. Bifurcations in Boltzmann-Langevin one body dynam-

ics for fermionic systems. Phys. Lett. B, 726(1):382 – 386, 2013.

[25] Y.-X. Zhang, Y.-J. Wang, M. Colonna, P. Danielewicz, A. Ono, M. B. Tsang, H. Wolter,
J. Xu, L.-W. Chen, D. Cozma, et al. Comparison of heavy-ion transport simulations:
Collision integral in a box. Phys. Rev. C, 97(3):034625, 2018.

[26] Y. M. Engel, D. M. Brink, K. Goeke, S. J. Krieger, and D. Vautherin. Time-dependent
Hartree-Fock theory with Skyrme’s interaction. Nucl. Phys. A, 249(2):215–238, 1975.

[27] M. C. Cambiaggio, G. G. Dussel, and M. Saraceno. Time-dependent Hartree-Fock-
Bogoliubov description of the pairing interaction. Nucl. Phys. A, 415(1):70–92, 1984.

[28] D. Lacroix and S. Ayik. Stochastic quantum dynamics beyond mean ﬁeld. Eur. Phys.

J. A, 50(6):95, 2014.

[29] F. V. De Blasio, W. Cassing, M. Tohyama, P. F. Bortignon, and R. A. Broglia. Non-
perturbative study of the damping of giant resonances in hot nuclei. Phys. Rev. Lett.,
68(11):1663, 1992.

[30] M. Tohyama and A. S. Umar. Fusion window problem in time-dependent Hartree-Fock

theory revisited. Phys. Rev. C, 65(3):037601, 2002.

[31] S. Ayik. A stochastic mean-ﬁeld approach for nuclear dynamics. Phys. Lett. B,

658(4):174–179, 2008.

[32] D. Lacroix and S. Ayik. Stochastic quantum dynamics beyond mean ﬁeld. Eur. Phys.

J. A, 50(6):95, Jun 2014.

[33] Y. P. Viyogi, T. J. M. Symons, P. Doll, et al. Fragmentations of 40Ar at 213 Mev/nu-

cleon. Phys. Rev. Lett., 42:33–36, Jan 1979.

[34] D. R. Bowman, G. F. Peaslee, R. T. de Souza, N. Carlin, C. K. Gelbke, W. G. Gong,
Y. D. Kim, M. A. Lisa, W. G. Lynch, L. Phair, M. B. Tsang, C. Williams, N. Colonna,
K. Hanold, M. A. McMahan, G. J. Wozniak, L. G. Moretto, and W. A. Friedman.
Multifragment disintegration of the 129Xe+197Au system at E/A=50 Mev. Phys. Rev.
Lett., 67:1527–1530, Sep 1991.

[35] W. Bauer, G. F. Bertsch, and S. Das Gupta. Fluctuations and clustering in heavy-ion

collisions. Phys. Rev. Lett., 58:863–866, Mar 1987.

109

[36] H. Kruse, B. V. Jacak, and H. Stöcker. Microscopic theory of pion production and

sidewards ﬂow in heavy-ion collisions. Phys. Rev. Lett., 54:289–292, Jan 1985.

[37] A. Ono, H. Horiuchi, T. Maruyama, and A. Ohnishi. Antisymmetrized version of
molecular dynamics with two-nucleon collisions and its application to heavy-ion reac-
tions. Prog. Theor. Phys., 87(5):1185–1206, 1992.

[38] A. Ono, H. Horiuchi, T. Maruyama, and A. Ohnishi. Momentum distribution of
fragments in heavy-ion reactions: Dependence on the stochastic collision process. Phys.
Rev. C, 47:2652–2660, Jun 1993.

[39] J. Su, F.-S. Zhang, and B.-A. Bian. Odd-even eﬀect in heavy-ion collisions at inter-

mediate energies. Phys. Rev. C, 83:014608, Jan 2011.

[40] J. Su and F.-S. Zhang. Non-equilibrium and residual memory in momentum space of
fragmenting sources in central heavy-ion collisions. Phys. Rev. C, 87:017602, Jan 2013.

[41] J. Su, K. Cherevko, W.-J. Xie, and F.-S. Zhang. Nonisotropic and nonsingle explosion
in central 129Xe + 120Sn collisions at 50–125 Mev/nucleon. Phys. Rev. C, 89:014619,
Jan 2014.

[42] C. Hartnack, Z.-X. Li, L. Neise, et al. Quantum molecular dynamics a microscopic

model from UNILAC to CERN energies. Nucl. Phys. A, 495(1):303 – 319, 1989.

[43] M. Papa, T. Maruyama, and A. Bonasera. Constrained molecular dynamics approach

to fermionic systems. Phys. Rev. C, 64:024612, Jul 2001.

[44] M. Papa. Many-body correlations in semiclassical molecular dynamics and Skyrme

interaction. Phys. Rev. C, 87:014001, Jan 2013.

[45] Y.-X. Zhang and Z.-X. Li. Elliptic ﬂow and system size dependence of transition

energies at intermediate energies. Phys. Rev. C, 74:014602, Jul 2006.

[46] Y.-X. Zhang, Z.-X. Li, and P. Danielewicz. In-medium NN cross sections determined
from the nuclear stopping and collective ﬂow in heavy-ion collisions at intermediate
energies. Phys. Rev. C, 75:034615, Mar 2007.

[47] Y.-X. Zhang, P. Danielewicz, M. Famiano, et al. The inﬂuence of cluster emission
and the symmetry energy on neutron-proton spectral double ratios. Phys. Lett. B,
664(1):145 – 148, 2008.

110

[48] Z.-Q. Feng. Momentum dependence of the symmetry potential and its inﬂuence on

nuclear reactions. Phys. Rev. C, 84:024610, Aug 2011.

[49] Z.-Q. Feng. Nuclear in-medium eﬀects and collective ﬂows in heavy-ion collisions at

intermediate energies. Phys. Rev. C, 85:014604, Jan 2012.

[50] X. G. Cao, G. Q. Zhang, X. Z. Cai, et al. Roles of deformation and orientation in
heavy-ion collisions induced by light deformed nuclei at intermediate energy. Phys.
Rev. C, 81:061603, Jun 2010.

[51] G. Q. Zhang, Y. G. Ma, X. G. Cao, et al. Uniﬁed description of nuclear stopping in
central heavy-ion collisions from 10A Mev to 1.2A Gev. Phys. Rev. C, 84:034612, Sep
2011.

[52] W. B. He, Y. G. Ma, X. G. Cao, X. Z. Cai, and G. Q. Zhang. Giant dipole resonance as a
ﬁngerprint of α clustering conﬁgurations in 12C and 16O. Phys. Rev. Lett., 113:032506,
Jul 2014.

[53] D. T. Khoa, N. Ohtsuka, M. A. Matin, et al. In-medium eﬀects in the description of
heavy-ion collisions with realistic NN interactions. Nucl. Phys. A, 548(1):102 – 130,
1992.

[54] V. S. Uma Maheswari, C. Fuchs, A. Faessler, et al. In-medium dependence and coulomb
eﬀects of the pion production in heavy-ion collisions. Nucl. Phys. A, 628(4):669 – 685,
1998.

[55] K. Shekhter, C. Fuchs, A. Faessler, M. Krivoruchenko, and B. Martemyanov. Dilepton
production in heavy-ion collisions at intermediate energies. Phys. Rev. C, 68:014904,
Jul 2003.

[56] M. D. Cozma, Y. Leifels, W. Trautmann, Q. Li, and P. Russotto. Toward a model-
independent constraint of the high-density dependence of the symmetry energy. Phys.
Rev. C, 88:044912, Oct 2013.

[57] Q. Li, Z. Li, S. Soﬀ, M. Bleicher, and H. Stocker. Probing the density dependence of
the symmetry potential at low and high densities. Phys. Rev. C, 72:034613, Sep 2005.

[58] Q.-F. Li, C.-W. Shen, C.-C. Guo, et al. Nonequilibrium dynamics in heavy-ion col-
lisions at low energies available at the GSI Schwerionen Synchrotron. Phys. Rev. C,
83:044617, Apr 2011.

111

[59] P. Napolitani and M. Colonna. Frustrated fragmentation and re-aggregation in nuclei:

A non-equilibrium description in spallation. Phys. Rev. C, 92:034607, Sep 2015.

[60] T. Gaitanos, A. B. Larionov, H. Lenske, and U. Mosel. Breathing mode in an improved

transport approach. Phys. Rev. C, 81:054316, May 2010.

[61] A. B. Larionov, T. Gaitanos, and U. Mosel. Kaon and hyperon production in

antiproton-induced reactions on nuclei. Phys. Rev. C, 85:024614, Feb 2012.

[62] O. Buss, T. Gaitanos, K. Gallmeister, et al. Transport-theoretical description of nuclear

reactions. Phys. Rep., 512(1):1 – 124, 2012.

[63] F.-S. Zhang and E. Suraud. Analysis of multifragmentation in a Boltzmann-Langevin

approach. Phys. Rev. C, 51:3201–3210, Jun 1995.

[64] W.-J. Xie, J. Su, L. Zhu, and F.-S. Zhang. Neutron-proton eﬀective mass splitting in

a Boltzmann-Langevin approach. Phys. Rev. C, 88:061601, Dec 2013.

[65] W.-J. Xie, J. Su, L. Zhu, and F.-S. Zhang. Symmetry energy and pion production in

the Boltzmann-Langevin approach. Phys. Lett. B, 718(4):1510 – 1514, 2013.

[66] P. Danielewicz. Determination of the mean-ﬁeld momentum-dependence using elliptic

ﬂow. Nucl. Phys. A, 673(1):375 – 410, 2000.

[67] B.-A. Li, L.-W. Chen, and C. M. Ko. Recent progress and new challenges in isospin

physics with heavy-ion reactions. Phys. Rep., 464(4):113 – 281, 2008.

[68] C. Fuchs and H. H. Wolter. The relativistic Landau-Vlasov method in heavy-ion

collisions. Nucl. Phys. A, 589(4):732 – 756, 1995.

[69] T. Gaitanos, M. Di Toro, S. Typel, et al. On the Lorentz structure of the symmetry

energy. Nucl. Phys. A, 732:24 – 48, 2004.

[70] G. Ferini, T. Gaitanos, M. Colonna, M. DiToro, and H. H. Wolter. Isospin eﬀects on
subthreshold kaon production at intermediate energies. Phys. Rev. Lett., 97:202301,
Nov 2006.

[71] C. M. Ko and Qi Li. Relativistic Vlasov-Uehling-Uhlenbeck model for heavy-ion colli-

sions. Phys. Rev. C, 37:2270–2273, May 1988.

112

[72] C. M. Ko and G. Q. Li. Medium eﬀects in high energy heavy-ion collisions. J. Phys.

G, 22(12):1673, 1996.

[73] T. Song and C. M. Ko. Modiﬁcations of the pion-production threshold in the nu-
clear medium in heavy ion collisions and the nuclear symmetry energy. Phys. Rev. C,
91:014901, Jan 2015.

[74] A. Guarnera, M. Colonna, and Ph. Chomaz. 3D stochastic mean-ﬁeld simulations of

the spinodal fragmentation of dilute nuclei. Phys. Lett. B, 373(4):267 – 274, 1996.

[75] M. Colonna. Fluctuations and symmetry energy in nuclear fragmentation dynamics.

Phys. Rev. Lett., 110:042701, Jan 2013.

[76] S. Ayik and C. Gregoire. Fluctuations of single-particle density in nuclear collisions.

Phys. Lett. B, 212(3):269 – 272, 1988.

[77] E. Suraud, S. Ayik, J. Stryjewski, and M. Belkacem. The Boltzmann-Langevin equa-
tion and its application to intermediate-mass-fragment production. Nucl. Phys. A,
519(1):171 – 182, 1990.

[78] P. Danielewicz and G.F. Bertsch. Production of deuterons and pions in a transport

model of energetic heavy-ion reactions. Nucl. Phys. A, 533(4):712 – 748, 1991.

[79] R. J. Lenk and V. R. Pandharipande. Nuclear mean ﬁeld dynamics in the lattice

Hamiltonian Vlasov method. Phys. Rev. C, 39:2242–2249, Jun 1989.

[80] P. Danielewicz. Nonrelativistic and relativistic Landau/Fokker-Planck equation for

arbitrary statistics. Physica A, 100(1):167 – 182, 1980.

[81] R. F. Pawula. Approximation of the linear Boltzmann equation by the Fokker-Planck

equation. Phys. Rev., 162:186–188, Oct 1967.

[82] N. L. Balazs. The Fokker-Plank and Kramers-Chandrasekhar equations for semiclas-

sical bosons and fermions. Physica A, 94(3):474 – 480, 1978.

[83] P. Ring and P. Schuck. The nuclear many-body problem. Springer Science & Business

Media, 2004.

[84] T. Hastie, R. Tibshirani, and J. Friedman. The Elements of Statistical Learning: Data

Mining, Inference, and Prediction. Elsevier, 2017.

113

[85] G. Q. Li, R. Machleidt, and Y. Z. Zhuo. Self-consistent relativistic calculation of

nucleon mean free path. Phys. Rev. C, 48:1062–1068, Sep 1993.

[86] Ph. Chomaz, M. Colonna, and J. Randrup. Nuclear spinodal fragmentation. Phys.

Rep., 389(5):263 – 440, 2004.

[87] J. Rizzo, M. Colonna, and A. Ono. Comparison of multifragmentation dynamical

models. Phys. Rev. C, 76:024611, Aug 2007.

[88] A. Ono, S. Hudan, A. Chbihi, and J. D. Frankland. Compatibility of localized wave
packets and unrestricted single particle dynamics for cluster formation in nuclear col-
lisions. Phys. Rev. C, 66:014603, Jul 2002.

[89] W. J. Dou, G. H. Miao, and J. E. Subotnik. Born-Oppenheimer dynamics, elec-
tronic friction, and the inclusion of electron-electron interactions. Phys. Rev. Lett.,
119:046001, Jul 2017.

[90] H. Lin, H. Mahzoon, A. Rios, and P. Danielewicz. Dynamics of one-dimensional cor-
related nuclear systems within non-equilibrium Green’s function theory. Ann. Phys.,
420:168272, 2020.

[91] D. J. Thouless. Stability conditions and nuclear rotations in the Hartree-Fock theory.

Nucl. Phys., 21:225–232, 1960.

[92] H.-J. Mang. The self-consistent single-particle model in nuclear physics. Phys. Rep.,

18(6):325–368, 1975.

[93] D. Gogny and P.-L. Lions. Hartree-Fock theory in nuclear physics. ESAIM: Math.

Model. Numer. Anal., 20(4):571–637, 1986.

[94] C. Ordóñez, L. Ray, and U. van Kolck. Nucleon-nucleon potential from an eﬀective

chiral Lagrangian. Phys. Rev. Lett., 72(13):1982, 1994.

[95] R. Machleidt and D. R. Entem. Towards a consistent approach to nuclear structure:
EFT of two- and many-body forces. J. Phys. G: Nucl. Part. Phys., 31(8):S1235, 2005.

[96] R. Machleidt and D. R. Entem. Chiral eﬀective ﬁeld theory and nuclear forces. Phys.

Rep., 503(1):1–75, 2011.

[97] T. Tamura, W. R. Coker, and F. Rybicki. Distorted Wave Born Approximation for

nuclear reactions. Technical report, Univ. of Texas, Austin, 1971.

114

[98] L. L. Lee Jr, J. P. Schiﬀer, B. Zeidman, G. R. Satchler, R. M. Drisko, and R. H. Bassel.
40ca(d, p) 41Ca, a test of the validity of the Distorted-Wave Born Approximation. Phys.
Rev., 136(4B):B971, 1964.

[99] T. Tamura. Coupled-channel approach to nuclear reactions. Ann. Rev. Nucl. Sc.,

19(1):99–138, 1969.

[100] I. J. Thompson. Coupled reaction channels calculations in nuclear physics. Comp.

Phys. Rep., 7(4):167–212, 1988.

[101] N. Austern, Y. Iseri, M. Kamimura, M. Kawai, G. Rawitscher, and M. Yahiro.
three-body models of

Continuum-discretized coupled-channels calculations
deuteron-nucleus reactions. Phys. Rep., 154(3):125–204, 1987.

for

[102] A. S. Umar and V. E. Oberacker. Three-dimensional unrestricted time-dependent
Hartree-Fock fusion calculations using the full Skyrme interaction. Phys. Rev. C,
73(5):054607, 2006.

[103] B. Avez, C. Simenel, and Ph. Chomaz. Pairing vibrations study with the time-

dependent Hartree-Fock-Bogoliubov theory. Phys. Rev. C, 78(4):044318, 2008.

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

[105] J. A. Maruhn, P.-G. Reinhard, P. D. Stevenson, and A. S. Umar. The TDHF code

sky3d. Comp. Phys. Comm., 185(7):2195–2216, 2014.

[106] N. Loebl, A. S. Umar, J. A. Maruhn, P.-G. Reinhard, P. D. Stevenson, and V. E. Ober-
acker. Single-particle dissipation in a time-dependent Hartree-Fock approach studied
from a phase-space perspective. Phys. Rev. C, 86:024608, Aug 2012.

[107] J. Aichelin and G. Bertsch. Numerical simulation of medium energy heavy ion reac-

tions. Phys. Rev. C, 31(5):1730, 1985.

[108] Y.-X. Zhang, Y.-J. Wang, M. Colonna, P. Danielewicz, A. Ono, et al. Comparison of
heavy-ion transport simulations: Collision integral in a box. Phys. Rev. C, 97:034625,
Mar 2018.

[109] J. Aichelin and H. Stoecker. Quantum molecular dynamics: A novel approach to

N-body correlations in heavy-ion collisions. Phys. Lett. B, 176(1-2):14–19, 1986.

115

[110] J. Aichelin. Quantum molecular dynamics: A dynamical microscopic N-body approach
to investigate fragment formation and the nuclear equation of state in heavy-ion colli-
sions. Phys. Rep., 202(5-6):233–360, 1991.

[111] Y. Kanada-En’yo, H. Horiuchi, and A. Ono. Structure of Li and Be isotopes studied

with antisymmetrized molecular dynamics. Phys. Rev. C, 52(2):628, 1995.

[112] A. Ono and H. Horiuchi. Antisymmetrized molecular dynamics for heavy-ion collisions.

Prog. Part. Nucl. Phys., 53(2):501–581, 2004.

[113] R. Wada, K. Hagel, J. Cibor, M. Gonin, Th. Keutgen, M. Murray, J. B. Natowitz,
A. Ono, J. C. Steckmeyer, A. Kerambrum, et al. Reaction mechanisms and multifrag-
mentation processes in 64Zn + 58Ni at 35 A – 79 A Mev. Phys. Rev. C, 62(3):034601,
2000.

[114] B. Borderie, G. Tăbăcaru, Ph. Chomaz, M. Colonna, A. Guarnera, M. Pârlog, M. F.
Rivet, G. Auger, Ch. O. Bacri, N. Bellaize, et al. Evidence for spinodal decomposition
in nuclear multifragmentation. Phys. Rev. Lett., 86(15):3252, 2001.

[115] J. Hong and P. Danielewicz. Subthreshold pion production within a transport descrip-

tion of central Au+ Au collisions. Phys. Rev. C, 90(2):024605, 2014.

[116] A. Rios, B. Barker, M. Buchler, and P. Danielewicz. Towards a nonequilibrium Green’s
function description of nuclear reactions: One-dimensional mean-ﬁeld dynamics. Ann.
Phys. (N.Y.), 326(5):1274–1319, 2011.

[117] P Napolitani and M Colonna. Bifurcations in Boltzmann-Langevin one body dynamics

for fermionic systems. Phys. Lett. B, 726(1-3):382–386, 2013.

[118] H. Lin and P. Danielewicz. One-body Langevin dynamics in heavy-ion collisions at

intermediate energies. Phys. Rev. C, 99(2):024612, 2019.

[119] P. Danielewicz. Quantum theory of nonequilibrium processes, I. Ann. Phys. (N.Y.),

152(2):239–304, 1984.

[120] L. P. Kadanoﬀ and G. Baym. Quantum statistical mechanics. New York, 1962.

[121] M. B. Tsang, J. R. Stone, F. Camera, P. Danielewicz, S. Gandolﬁ, K. Hebeler, C. J.
Horowitz, J. Lee, W. G. Lynch, Z. Kohley, et al. Constraints on the symmetry energy
and neutron skins from experiments and theory. Physical Review C, 86(1):015803,
2012.

116

[122] A. Rios, B. Barker, and P. Danielewicz. Towards a nonequilibrium Green’s function
description of nuclear reactions. In J. Phys.: Conf. Ser., volume 427, page 012010.
IOP Publishing, 2013.

[123] M. H. Mahzoon, P. Danielewicz, and A. Rios. Correlations within the non-equilibrium
Green’s function method. In AIP Conf. Proc., volume 1912, page 020012. AIP Pub-
lishing, 2017.

[124] H. Mahzoon, P. Danielewicz, and A. Rios. Nuclear slabs with Green’s functions: mean
ﬁeld and short-range correlations. Eur. Phys. J. Special Topics, 227(15-16):1949–1958,
2019.

[125] N. Schlünzen, S. Hermanns, M. Scharnke, and M. Bonitz. Ultrafast dynamics of
strongly correlated fermions — nonequilibrium Green functions and self energy ap-
proximations. Journal of Physics: Condensed Matter, 32(10):103001, dec 2019.

[126] A. Bulgac, M. M. Forbes, K. J. Roche, and G. Wlazłowski. Quantum friction: Cooling
quantum systems with unitary time evolution. arXiv preprint arXiv:1305.6891, 2013.

[127] H. Bruus and K. Flensberg. Many-Body Quantum Theory in Condensed Matter

Physics: An Introduction. Oxford University Press, 2004.

[128] P. Danielewicz. Quantum theory of nonequilibrium processes II. Application to nuclear

collisions. Ann. Phys. (N.Y.), 152(2):305–326, 1984.

[129] O. Hen, G. A. Miller, E. Piasetzky, and L. B. Weinstein. Nucleon-nucleon correlations,
short-lived excitations, and the quarks within. Rev. Mod. Phys., 89:045002, Nov 2017.

[130] A. Rios, A. Polls, and W. H. Dickhoﬀ. Density and isospin-asymmetry dependence of

high-momentum components. Phys. Rev. C, 89(4):044303, 2014.

[131] W. H. Dickhoﬀ and D. Van Neck. Many-Body Theory Exposed!: Propagator Descrip-
tion of Quantum Mechanics in Many-Body Systems Second Edition. World Scientiﬁc
Publishing Company, 2008.

[132] R. Chandra, L. Dagum, D. Kohr, R. Menon, D. Maydan, and J. McDonald. Parallel

Programming in OpenMP. Morgan Kaufmann, 2001.

[133] N.G. Van Kampen. Chapter IX - The Langevin Approach. Elsevier, 2007.

117

[134] J. K. Hunter. Lecture 5: Stochastic processes. lecture notes (Math 280, University of

California, Davis, Spring Quarter, 2009), 2009.

118