DISSIPATION AND DYNAMICS IN QUANTUM MANY-BODY SYSTEMS
By
Brent Wendolyn Barker

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Physics—Doctor of Philosophy
2014

ABSTRACT
DISSIPATION AND DYNAMICS IN QUANTUM MANY-BODY SYSTEMS
By
Brent Wendolyn Barker
In this thesis, we simulate the time evolution of quantum many-body systems and use
comparisons to experimental data in order to learn more about the properties of nuclear
matter and understand better the dynamical processes in central nuclear collisions. We
further advance the development of a nonequilibrium Green’s function description of both
central nuclear collisions and Bose-Einstein Condensates.
First in the thesis, we determine the viscosity of nuclear matter by adjusting the inmedium nucleon-nucleon cross section in our BUU transport model until the simulation
results match experimental data on nuclear stopping in central nuclear collisions at intermediate energies. Then we use that cross section to calculate the viscosity self-consistently.
We also calculate the ratio of shear viscosity to entropy density to determine how close the
system is to the proposed universal quantum lower limit.
Next, we use the same BUU transport model to isolate the protons emitted early in
a central nuclear collision at intermediate energy, as predicted in the model, using a filter on high transverse momentum, and we show the effect on the source function. We
predict a recontraction of protons at late times in the central collision of 112 Sn +112 Sn at
50 MeV/nucleon that results in a resurgence of emission of protons and show how to use the
transverse momentum filter and the source function to test this prediction in experiment.
Next, we develop an early implementation of a more fully quantal transport model than
the BUU equations, with our sights set on solving central nuclear collisions in 3D using

nonequilibrium Green’s functions. In our 1D, mean field, density matrix model, we demonstrate the initial state preparation and collision of 1D nuclear “slabs”. With the aim of
reducing the computational cost of the calculation, we show that we can neglect far offdiagonal elements in the density matrix without afffecting the one-body observables.
Further, we describe a method of recasting the density matrix in a rotated coordinate
system, enabling us to not only ignore the irrelevant matrix elements in the time evolution,
but also avoid computing them completely, reducing the computational cost. As an added
benefit, we find that the rotation allows us to partially decouple the position and momentum
discretization, permitting access to arbitrary regimes of kinetic energy without altering the
resolution and range of the 1D box in position space.
Finally, we exhibited the wide applicability of this density matrix approach by applying
it to a system of 2000 ultracold 87 Rb atoms in a Bose-Einstein condensate, as described by
the Gross-Pitaevskii equation, successfully acheiving a stable state in a harmonic oscillator
trap.

Copyright by
BRENT WENDOLYN BARKER
2014

Dedicated to my parents, Arlene and Phil. Without their encouragement from a young
age, I would never have come so far.

v

ACKNOWLEDGMENTS

There have been so many with whom I shared part of this journey, and who have been
instrumental in one part or another of my joy, fulfillment, and success during my time
at MSU. I can only hope to thank a small fraction of them, and even then, the list will
inevitably be biased towards those who have influenced me more recently. That said, I’d like
to acknowledge at least these folks:
• My advisor, Pawel Danielewicz, whose patience is as legendary as his insight, and whose
availability is exceedingly convenient, allowing for swift collaboration, even while he
spends half of his time on the other side of the world.
• Professional mentors: Brian O’Shea for our discussions about teaching and physics
education research. Morten Hjorth-Jensen, for his eternal optimism, compassion, and
humility. Filomena Nunes, for assisting with time management skills.
• My longterm officemate family Jun Hong and Ngoc Bic Nguyen, for providing a supportive environment in which to work. In particular, I want to thank Jun for her
unyielding, compassionate support and encouragement during our daily meetings, and
Bic for our walks and long and fruitful discussions.
• Micha Kilburn, a mentor and friend.
• My found family: Matt DeNinno, Sofia Maystrenko, and Katie Thompson. We’ve
grown so much together, and they’ve helped to keep me sane and grounded in my
durance, and have been absolutely vital to my continued success.

vi

• First year cohort: Andrew Baczewski, Brad Barquest, Marshall Bremer, Matt DeNinno,
Shea Mosby, Josh Vredevoogd, Kathy Walsh, Emi Wang. Friendships born in the crucible of graduate courses are enduring.
• Rhiannon Meharchand, for supporting and believing in me.

Her hard work and

perserverance is an inspiration.
• Collaborators: Arnau “Jacket Brother” Rios and Giuseppe Verde, for insightful discussions and sharing the load.
I feel I must be forgetting someone important. I owe a special debt of gratitude towards
that person or persons. They know who they are.

vii

TABLE OF CONTENTS

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

x

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

xi

Chapter 1 Introduction . . . . . . . . . . .
1.1 Dynamics of central collisions . . . . .
1.2 Boltmann-Uehling-Uhlenbeck equation
1.2.1 Impact parameter comparison .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

1
3
6
7

Chapter 2 Dissipation . . . . . . . . . . . . . . . . . . . . . . .
2.1 The NN cross section in the nuclear medium . . . . . . . .
2.1.1 Tempered . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 Rostock . . . . . . . . . . . . . . . . . . . . . . . .
2.1.3 Fuchs . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Stopping observables . . . . . . . . . . . . . . . . . . . . .
2.2.1 Linear momentum transfer . . . . . . . . . . . . . .
2.2.2 Rapidity variance ratio varxz . . . . . . . . . . . .
2.2.3 Isospin tracer . . . . . . . . . . . . . . . . . . . . .
2.2.4 Sensitivity to Equation of State . . . . . . . . . . .
2.2.5 Optimal cross section reduction . . . . . . . . . . .
2.3 Viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Viscosity from BUU . . . . . . . . . . . . . . . . .
2.3.2 Lower quantum limit of ratio of viscosity to entropy
2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
density
. . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

8
9
10
10
11
12
12
16
21
23
24
26
26
28
30

Chapter 3 Dynamics . . . . . . . . . . . . . . . . . .
3.1 Introduction . . . . . . . . . . . . . . . . . . . .
3.2 HBT effect . . . . . . . . . . . . . . . . . . . .
3.2.1 Correlation function . . . . . . . . . . .
3.2.2 Source function . . . . . . . . . . . . . .
3.3 Timescale of emission . . . . . . . . . . . . . . .
3.3.1 Interpreting the source function . . . . .
3.3.2 Time dependence of the source function
3.4 Late-time resurgent emission . . . . . . . . . . .
3.4.1 Analysis of second emission . . . . . . .
3.4.2 Detection in experiment . . . . . . . . .
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

31
31
34
35
36
38
40
42
45
45
51
59

viii

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

Chapter 4 Towards a fully quantal approach . . . . .
4.1 Kadanoff-Baym equations . . . . . . . . . . . . . .
4.2 Mean field approximation . . . . . . . . . . . . . .
4.3 Collision of one-dimensional nuclear slabs . . . . . .
4.3.1 Preparation of the initial state . . . . . . . .
4.3.2 Dynamics of colliding slabs . . . . . . . . . .
4.3.3 Off-diagonal structure of the density matrix
4.4 Testing importance of off-diagonal elements . . . .
4.5 Rotated density matrix . . . . . . . . . . . . . . . .
4.5.1 Fourier transform in the new coordinates . .
4.5.2 Discrete Fourier transform . . . . . . . . . .
4.5.3 Testing the rotated matrix . . . . . . . . . .
4.6 Bose-Einstein Condensates . . . . . . . . . . . . . .
4.6.1 Modeling with the Gross-Pitaevskii equation
4.6.2 Initial state . . . . . . . . . . . . . . . . . .
4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

. 61
. 64
. 66
. 67
. 68
. 72
. 74
. 77
. 82
. 83
. 84
. 85
. 92
. 92
. 95
. 100

Chapter 5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Appendix A Equations for Source Imaging Chapter . . . . . . . . . . .
Appendix B Parameterizations of the source function . . . . . . . . . .
Appendix C Existence of the BEC in one dimension . . . . . . . . . . .
Appendix D Numerical stability in the density matrix time evolution
BIBLIOGRAPHY

. . . . . . . . . . . . .

ix

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

103
104
105
108
110

. . . . . . . . . . . . . . . . . 113

LIST OF TABLES

Table 2.1

Table 3.1

Summary of cross section determination results. No cross section
reduction is favored universally, though Tempered with ν = 0.6 or
0.8 appears to be the best compromise. . . . . . . . . . . . . . . . .

26

Fit parameters for source functions. λgau and λexp are the fractions
of particle pairs that contribute to the first and second terms of the
equation. Rgau and Rexp are the widths of the Gaussian and exponential components, respectively. The fraction of particle pairs that
contribute to the Gaussian shape, fgau , is higher for the smaller system, while the larger system has a correspondingly longer exponential
tail. Uncertainties are statistical from fit and do not include uncertainties of underlying model. . . . . . . . . . . . . . . . . . . . . . .

41

x

LIST OF FIGURES

Figure 1.1

Figure 2.1

Figure 2.2

Figure 2.3

Figure 2.4

Figure 2.5

Figure 2.6

Figure 2.7

Stages of an intermediate-energy nuclear collision. From top to bottom: approach, compression, expansion (figure from Ref. [1]). . . . .

5

(For interpretation of the references to color in this and all other
figures, the reader is referred to the electronic version of this dissertation.) Schematic of asymmetric collision. Projectile transfers
momentum to the target. To measure linear momentum transfer, the
longitudinal velocity of the target-like fragment is compared to the
velocity of the center-of-mass of the collision. . . . . . . . . . . . . .

13

Linear momentum transfer for 40 Ar+Cu. Lines show the effects of different in-medium reductions of the NN cross section. The “Tempered”
reductions are marked with their tunable parameter ν. Symbols are
experimental data [2]. . . . . . . . . . . . . . . . . . . . . . . . . . .

14

Linear momentum transfer for 40 Ar+Ag. Lines show the effects of different in-medium reductions of the NN cross section. The “Tempered”
reductions are marked with their tunable parameter ν. Symbols are
experimental data [2]. . . . . . . . . . . . . . . . . . . . . . . . . . .

14

Linear momentum transfer for 40 Ar+Au. Lines show the effects of different in-medium reductions of the NN cross section. The “Tempered”
reductions are marked with their tunable parameter ν. Symbols are
experimental data [2]. . . . . . . . . . . . . . . . . . . . . . . . . . .

15

Stopping observable varxz for protons in Au + Au collisions at different beam energies at bred < 0.15. Lines show the effects of different
in-medium reductions of the NN cross section. The “Tempered” reductions are marked with their tunable parameter ν. Symbols are
experimental data from FOPI [3]. . . . . . . . . . . . . . . . . . . .

17

Same as Fig. 2.5, except that particle coalescence is enabled in the
simulation, so comparison to different particle species in experiment
is possible. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

Stopping observable varxz for Ca+Ca collisions at bred < 0.15. Lines
and symbols are the same as Fig. 2.5. . . . . . . . . . . . . . . . . .

20

xi

Figure 2.8

Figure 2.9

Figure 2.10

Figure 2.11

Figure 2.12

Stopping observable varxz for Ca+Ca collisions at bred < 0.3. Lines
and symbols are the same as Fig. 2.5. . . . . . . . . . . . . . . . . .

20

Isospin tracer observable for central collisions of 96 Zr +96 Ru (bottom)
and 96 Ru +96 Zr (top) at 400 MeV/nucleon vs. scaled center-of-mass
rapidity, where y0 = −1 is the initial projectile rapidity. Data is from
FOPI [4]. A value of Rz closer to zero indicates a higher degree of
mixing, and thus stopping. . . . . . . . . . . . . . . . . . . . . . . .

22

EOS dependence of stopping observables. S (H) refers to a soft
(stiff) compressibility, while M refers to the inclusion of momentumdependence in the mean field. SM and H are both reasonable combinations, and they give similar predictions for stopping observables
varxz and isospin tracer Rz . However, LMT is as sensitive to these
settings as it is to cross section reductions, so that observable is less
reliable for selecting the best reduction. . . . . . . . . . . . . . . . .

25

Viscosity for different temperatures, densities, and CS reductions.
Note that for ρ/ρ0 = 1, the Fuchs reduction overlaps with ν = 0.8 up
to about T = 40 MeV/kB , after which it overlaps with the Rostock
reduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

This work’s estimate for η/s (in open hexagons) with the ratio given
for other materials. RHIC has found a very low value, close to the
proposed quantum lower limit (filled, red hexagon). We find that
nuclear matter at intermediate energies approaches this value. The
lines are to guide the eye. Compilation and figure from Ref. [5]. . . .

29

Figure 3.1

Schematic of HBT effect. Two particles with momenta p1 and p2 are
emitted from the collision and observed at detectors A and B. If the
particles are identical, then the emission is a quantum superposition
of the left and right cases, leading to interference (figure from Ref. [6]). 34

Figure 3.2

Sample correlation function [7]. At large q, particles are not correlated
and the correlation function flattens to unity. . . . . . . . . . . . . .

36

Source functions from BUU transport code. Inset is a close-up in
log scale of the long-range part. The source for the smaller system
is more peaked at short range, and it falls off faster at long range,
compared to the larger system. . . . . . . . . . . . . . . . . . . . . .

39

Figure 3.3

xii

Figure 3.4

Source function S(r) from BUU transport simulation for 112 Sn +112 Sn
(triangles) and 40 Ar +45 Sc (squares). Also shown are curve fits using a short-range Gaussian and long-range exponential (seen in the
logscale plot inset). . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

Figure 3.5

Source function for 40 Ar +45 Sc constructed from only early time (green,
long dash) and only late time (blue, short dash). . . . . . . . . . . . 44

Figure 3.6

Proton emission versus time, from BUU transport simulation. Note
the clear separation between early and late emissions and rise in latetime emission in the heavier system. . . . . . . . . . . . . . . . . . .

46

Density contours at selected times after the beginning of collision
of (a) 40 Ar +45 Sc and (b) 112 Sn +112 Sn. The heavier system has a
recontraction of density at later times. Time is normalized so that
the centers of the nuclei cross at 51 fm/c (as predicted with Coulomb
effects, not including stopping). . . . . . . . . . . . . . . . . . . . .

47

2D histogram that shows where along the beam axis (z) protons are
emitted as time goes on, for a) 40 Ar +45 Sc and b) 112 Sn +112 Sn. In
the Sn collision, a second emission is quite pronounced, starting at
about 125 fm/c. At late times, the emission is seen to come from 4
distinct regions, which are the positions of high density in Fig. 3.7b.

48

Emission timeline for 112 Sn +112 112 compared to the same with an
ad hoc reduction in the NN cross section. The reduction suppresses
the late emission, while leaving the early emission mostly unaffected.

49

Density contours at selected times after the beginning of collision of
hoc reduction of σNN by a factor of 5. Unlike
the unreduced case seen in Fig. 3.7b, there is no formation of dense
fragments from the neck region at late times. . . . . . . . . . . . . .

50

Source function for 112 Sn +112 Sn calculated in a transport model with
a nominal NN cross section (triangles) and with an ad hoc reduction
in that cross section (squares). The lines represent a fit using the
Chung parameterization. . . . . . . . . . . . . . . . . . . . . . . . .

50

2D Histogram of p⊥ /m vs. t for protons emitted from collisions of
Also displayed are 1D histograms projected onto both
axes. Emphasis on early emission in an experiment is possible with a
filter on p⊥ /m 0.15 c. . . . . . . . . . . . . . . . . . . . . . . . . .

52

Same as Fig. 3.12, but for 112 Sn +112 Sn. The second surge in emission
is very visible at late times. . . . . . . . . . . . . . . . . . . . . . . .

53

Figure 3.7

Figure 3.8

Figure 3.9

Figure 3.10

112 Sn +112 Sn with an ad

Figure 3.11

Figure 3.12

40 Ar +45 Sc.

Figure 3.13

xiii

Source function for 40 Ar +45 Sc. Inset employs logscale to illuminate
the effects of the high p⊥ gate at far-away distances. . . . . . . . . .

54

Source function for 112 Sn +112 Sn. Inset is in logscale, to see the
effects of the high p⊥ gate at long range. . . . . . . . . . . . . . . .

55

Fit parameters for source functions created with filters on the emitted protons, for 40 Ar +45 Sc. The filter allowing only early t (open
squares) results in similar parameter values as for high p⊥ /m (open
triangles), showing that filtering on high p⊥ /m selects for early t.
Further, these filters enhance the Gaussian fraction fgau , showing
that the Gaussian component is primarily at early time. . . . . . . .

56

Figure 3.17

Same as Fig. 3.16, but for 112 Sn +112 Sn. . . . . . . . . . . . . . . .

57

Figure 3.18

Parameters describing the shape of the source function for 112 Sn +112 Sn,
calculated with a filter on transverse momentum, p⊥ /m > 0.15c
(squares) and without (triangles), as well as having the simulated
reaction calculated with a nominal NN cross section (open shapes)
and with an ad hoc reduction in that cross section, σ = σ0 /5 (filled
shapes), which leads to no late-time recontraction. In the σ/5 case,
the p⊥ cut does not affect the exponential fraction fexp like it does
in the nominal case. Therefore, the sensitivity of fexp to a cut in p⊥
is a signal of this recontraction. . . . . . . . . . . . . . . . . . . . . 58

Figure 3.19

Fraction of proton pairs, fexp , that contribute to the long-range exponential component of the source function, for the heavy system (left),
light system (middle), and heavy system with ad hoc reduction in
σNN (right). The first two collision systems produce a recontraction
in the neck region and resurgent emission. Here, a filter on high p⊥ /m
causes a preferential suppression of the exponential component. With
reduced σNN , the heavy system exhibits no re-emission, and the filter does not noticeably change fexp , suggesting that the filter is only
acting on 1 emission source, rather than 2. . . . . . . . . . . . . . .

59

Time evolution of the energy per particle (upper panel) and the size
of the slab (lower panel) when adiabatically switching from a starting
HO configuration to a final mean field state. Different values of the
transition time τ are considered. The inset in the top panel shows a
magnified portion of the time evolution of the energy. For reference,
mean field results from static Hartree-Fock solutions are shown as
straight solid (dark yellow) lines. . . . . . . . . . . . . . . . . . . .

71

Figure 3.14

Figure 3.15

Figure 3.16

Figure 4.1

xiv

Figure 4.2

Figure 4.3

Figure 4.4

Figure 4.5

Figure 4.6

Figure 4.7

Figure 4.8

Figure 4.9

Figure 4.10

Evolution of the center-of-mass density profile for a collision of two
A = 8 slabs at the collision energy of Ecm = 4 MeV/nucleon. . . . .

73

Intensity plots representing the real (upper panels) and the imaginary part (lower panels) of the scaled density matrix, ρ(x, x , t), for a
collision at Ecm = 4 MeV/nucleon. Off-diagonal correlations develop
as the slabs interpenetrate and fragments continue on. . . . . . . . .

75

Schematic of suppression of off-diagonal elements. White (non-shaded)
areas are heavily suppressed. The upper-left and lower-right corners
are not suppressed, to preserve periodicity in both coordinates. . . .

79

Time evolution of the energy per particle (upper panel) and of the
spatial extent (lower panel) in a collision between 1D nuclear slabs
at ECM /nucleon= 4 MeV, for different values of the cutoff parameter
x0 at fixed smoothness d0 = 2 fm. . . . . . . . . . . . . . . . . . . .

80

Schematic of density matrix in (x, x ) coordinates, with example of
periodic boundaries and the resulting physical system that is simulated. Each square contains an identical copy of the physical system.

86

With the density matrix rotated, the smallest rectangle that is periodic in both the average coordinate (left-to-right) and relative coordinate (down-to-up) is marked in dashed lines. . . . . . . . . . . . .

87

Schematic illustrating duplication of original density matrix. Regions
with the same pattern have the same exact part of the unrotated
(x, x ) density matrix in them. . . . . . . . . . . . . . . . . . . . . .

89

Rotation of matrix and creation of fictitious points. The x = x line is
shown as a dotted line, for orientation (xr = 0 line in rotated system).
When the matrix in (a) is rotated into (b) and the system is modeled
with a rectangular Fourier transform, extra grid points are created,
shown as open circles in (c). . . . . . . . . . . . . . . . . . . . . . .

90

Implementation of off-diagonal (xr -direction) suppression. Suppressed
areas are blue (shaded). Implementation is in (a), while the resulting
suppression of the original density matrix shown in (b) (compare to
Fig. 4.4). Truncation of the matrix is shown in (c), with the computed matrix bordered with a dot-dashed line and periodic extensions
shown with dashed borders. The empty space is not computed at all,
reducing the computation cost. This benefit is amplified geometrically with an increase in the number of dimensions. . . . . . . . . .

91

xv

Figure 4.11

Figure 4.12

Figure 4.13

Figure 4.14

Total energy of system while adiabatically switching from an external
HO potential to the Gross-Pitaevskii potential. A reasonably stable
energy is reached as the number of atoms in the condensate, N , increases. However, as seen in the inset, as N increases, the total energy
exhibits an increasing oscillation, suggesting higher excitations. . . .

97

Extent of system in position space while adiabatically switching from
an external HO potential to the Gross-Pitaevskii potential. The
ground state is more closely reached for a fewer number of particles
in the condensate, as seen by a slower, lower-amplitude oscillation at
late times than for larger numbers of particles. . . . . . . . . . . . .

97

Probability density at beginning (solid line, red) and end (dashed
line, green) of adiabatically switching from the harmonic oscillator
potential to the Gross-Pitaevskii potential. Since the potential energy is much stronger than the kinetic, the resulting distribution is
the Thomas-Fermi approximation, an upside-down parabola, as in
Eq. 4.22. The transition between the parabola and zero density, seen
inset, is smoother than the piecewise continuous solution and more
realistic. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

99

Density matrix for (a) the analytic harmonic oscillator ground state
before adiabatic switching and (b) after adiabatic switching to the
Gross-Piteavskii equation, yielding the Bose-Einstein Condensate. The
local density lies along the diagonal x = x , shown in Fig. 4.13. The
square support in the BEC shows a stronger correlation between
distant regions of the condensate, exhibiting stronger long-range coherence than in the HO case. . . . . . . . . . . . . . . . . . . . . . . 100

xvi

Chapter 1
Introduction
Central collisions between heavy nuclei are complex phenomena that involve many different
types of processes. These collisions are a means to learn about properties of high density, hot,
strongly interacting matter. The learning is contingent on proper modelling of the collisions.
The proper modelling is difficult because of the simultaneous effects of different processes.
Most often, it involves semiclassical approaches, and methods are often sourced from other
fields. In this work, we attempt to progress in several different directions in understanding
and modelling of the collisions.
To solve the full time-dependent, N -body wavefunction for the motion of N

30 parti-

cles interacting quantally would take far more time and resources than the world’s current
computing facilities have available. Instead, we use approximations. The quest, then, of
quantum many-body physics is to find those approximations that preserve enough of the
physics to give accurate predictions, and to use intuition and persistent random walks of
investigation to guide us to those approximations.
For central nuclear collisions, there are generally 2 approaches, microscopic and macroscopic. Macroscopically, the system can be modeled hydrodynamically, with the equation
of state and transport coefficients, like viscosity and heat conduction, as inputs to the calculation that determine how the matter interacts as it evolves forward in time. However,
the hydrodynamic approximation is strictly valid only when the particle mean free path is
much smaller than the size of the system. In nuclear systems, though, the mean free path
1

is comparable to the size of the system, and the hydrodynamic limit is not strictly valid.
Calculations have been performed using hydrodynamic models, and they do provide a coarse
description of the dynamics, but they do not provide as accurate a description as microscopic
models.
In the microscopic approach, the particles are modelled individually, and their interaction
with each other determines the dynamics. For example, in the Boltmann-Uehling-Uhlenbeck
(BUU) model, particles move in classical trajectories under the influence of a mean field
potential, and they undergo binary collisions with each other. The scattering cross sections
and mean field have parameters that are fit to experimental data, so these sorts of approaches are phenomenological. From this formalism, we can still derive, within the model,
the macroscopic quantities used as inputs for hydrodynamics. In particular, in Chapter 2,
we determine the viscosity of nuclear matter by adjusting the in-medium nucleon-nucleon
cross section within a BUU model until the simulations match data from nuclear stopping,
which is a good observable of how quickly momentum is transferred across the system, and,
hence, the viscosity.
As a macroscopic quantity, viscosity is defined for a bulk system, not for individual
particles. In particular, viscosity becomes useful for describing a system when there are
many particles interacting repeatedly in a small space. In nuclear collisions, this occurs early
in the collision, as the nuclei interpenetrate and mix with each other. To learn about the
evolution of the collision as particles spread out and start acting more like a gas than a liquid,
we must use other, microscopically-based methods. As the particles leave the high density
overlap region, if they are traveling close enough to each other, they attract or repel each
other, forming a correlation (or anti-correlation) between the momenta of those particles.
We use that correlation signature to learn about the shape of the particle emission region
2

in Chapter 3. Further, we use the shape of the resulting function, called a source function,
and compare it with the detailed dynamics in our microscopic model, to distinguish between
qualitatively different dynamics at late times in the collision. This gives us a closer look into
the actual timescale of these femtosecond processes.
The simulation mentioned in the previous few paragraphs, the BUU model, is semiclassical in nature. Since most nuclear collision simulations are semi-classical, it is difficult
to determine the validity of the approximation. In Chapter 4, we introduce the first steps of
an implementation of a fully quantal time evolution model for nuclear collisions, and we show
how we can smoothly approach the classical limit of the theory and explore the effect on the
dynamics. What’s more, since fully quantal approaches have a much greater computational
cost, in both time and storage, we demonstrate how to exploit the semi-classical nature of
the physical system to reduce that cost.
Finally, towards the end of Chapter 4, we use this new implementation to also describe
a 1D Bose-Einstein Condensate, a system of particles, in our case a gas of ultracold atoms,
all in the same, lowest energy state, and discuss some stability constraints and issues. Thus,
our model describes quantum phenomena across an energy range from GeV to neV.

1.1

Dynamics of central collisions

In a central (bred

0.3) collision1 between heavy nuclei (A

10) at intermediate energy

(Ebeam ≈ 50–2000 MeV/nucleon), the two nuclei approach each other and interpenetrate,
forming a hot, dense region in the middle. The part of each nuclei whose cross-sectional
1 The

centers of the nuclei are offset by a displacement b, called the impact parameter,
illustrated in Fig. 1.1. For comparison across different collision system sizes, we use a reduced
impact parameter, bred ≡ b/bmax , where bmax is the sum of the radii of the colliding nuclei.

3

area does not overlap with the other nucleus continues forward and is called the “spectator”, shearing off from the “participant” region2 . As the participant region compresses, the
nucleons that comprise it interact repeatedly, partially equilibrating. Then, after enough
longitudinal kinetic energy of the projectile is converted into internal degrees of freedom
(tending to randomize the velocities of the particles inside), it then expands, emitting particles and fragments. Those fragments and the spectators continue through freeze-out, when
they no longer interact with other particles. Then, if they are excited, they still will decay
or “evaporate”, emitting other particles.
Emission from nucleus-nucleus collisions is generally divided into two regimes. The first
is emission from the center of the collision very early in the time evolution, before the participant nucleons equilibrate. This is called nonequilibrium, pre-equilibrium, or dynamical
emission in the literature. The second regime is after the equilibrium, where emission occurs
as a result of decay of excited states. This is generally called slow, equilibrium, or evaporative emission. Since the entire collision event happens in ∼ 100 fm/c (1 fm/c ∼ 10−23 s),
while the typical time resolution in experiment is ∼ 1 ns, all particles from the event are
perceived to be detected simultaneously, mixing emission from both fast and slow processes.
Discerning fast versus slow decay requires comparison with a time-dependent theory, to see
which experimental observables will select for a particular time or shape of emission. We
implement a method for this in Chapter 3.

2 For

collisions with Ebeam
50 MeV/nucleon, the energy it takes to separate the peripheral, spectator nucleons from the interpenetrating ones (less than the binding energy per
nucleon, 8 MeV) is less than the kinetic energy of each nucleon in the local rest frame. As a
result, the spectator nucleons’ trajectories are not very affected by the collision.
4

Figure 1.1: Stages of an intermediate-energy nuclear collision. From top to bottom: approach, compression, expansion (figure from Ref. [1]).

5

1.2

Boltmann-Uehling-Uhlenbeck equation

To describe central nuclear reactions, we use the Boltzmann-Uehling-Uhlenbeck (BUU) set
of equations, describing the motion of a Wigner quasi-probability distribution in phase space,
f1 ≡ (r, p1 , t):
∂f1 ∂ p1 ∂f1 ∂ p1 ∂f1
+
−
=
∂t
∂p1 ∂r
∂r ∂p1

dp2 dΩ v12

dσ ˜ ˜
f f f f − f˜1 f˜2 f1 f2 ,
dΩ 1 2 1 2

(1.1)

where the subscript 1 is used only to avoid confusion with the other momenta under the
integral.
This equation starts with the single-particle Liouville equation on the l.h.s., describing the
∂ p

single-particle evolution of a phase space density in a mean field. For clarity, note that ∂p 1
1
∂ p

is the single-particle velocity, and ∂r1 is the force due to the mean field. The r.h.s. takes
into account the effect of binary collisions. The first term on the right describes particles
with p1 and p2 scattering into particles with p and p2 , thus increasing the occupancy of
f (gain). The second term describes, correspondingly, a decrease in the occupancy f . Here,
f˜a ≡ 1 − fa represents the Pauli principle blocking scattering into the final state pa . The
rate of scattering is governed by the NN cross section σ (here, a function of p1 , p2 , and
Ω, the angle between the relative momenta before and after the collision). It is this cross
section that is modified by in-medium effects in Section 2.1.
In Chapters 2 and 3, an implementation of a time-dependent solution to this equation
by Danielewicz and collaborators [8–13] is used to describe nuclear collisions. In this im˜ test particles are created in the simulation and tracked
plementation, for each nucleon, N
as particles. They move in classical trajectories under the influence of the mean field and
encounter random binary collisions with other test particles that are close to them in position
6

space. With an increase of test particle number, the simulation converges on a more average,
stable solution.

1.2.1

Impact parameter comparison

Throughout this work, we will be comparing our simulation results to experimental data. In
experiment, a range of impact parameters is selected for analysis. In the transport simulation,
the initial state is prepared with one specific impact parameter. To save computation time,
an effective impact parameter, beff , is chosen that represents the average cross-section for
the impact parameter range. For a given bmin and bmax , the effective impact parameter beff
is given by

πb2eff

πb2min + πb2max
=
2
(1.2)

beff =

1 2
+ b2max
b
2 min

For studies of central collisions, often experimental ranges set bmin = 0, so beff =
√
bmax / 2. This single impact parameter is used in many simulations.

7

Chapter 2
Dissipation
In this chapter, we focus on studying the transport of momentum through the equilibrated
region, quantified by the shear viscosity, in order to learn more about the fundamental
properties of nuclear matter.
Knowledge of the shear viscosity is important for understanding the stability of rotating
neutron stars, the formation of black holes, and the evolution of supernovae. Additionally,
there have been conjectures regarding a fundamental quantum lower limit on the ratio of
shear viscosity to entropy density (η/s) in a substance [14], thought to be approached in the
quark-gluon plasma and accessed in relativistic heavy ion collisions. The question remains
as to whether quark degrees of freedom are needed to approach such a limit. We will be able
to make a quantitative assessment of how close nuclear matter is to this “perfect liquid”.
Many groups have used giant resonances to determine the viscosity of nuclear matter
(see references in [15]). Several groups have studied the η/s ratio for nuclear collisions at
intermediate energies using different models, including statistical multifragmentation [16]
and quantum molecular dynamics (QMD) [17]. However, the group using QMD calculates
viscosity using an approach that is valid strictly within BUU transport models, and neither
group tunes their calculation to specifically reproduce observables directly related to viscosity.
In the present work, we use stopping, a phenomenom that is a measure of the degree to which the projectile nucleus transfers momentum to the target, to constrain the
in-medium nucleon-nucleon cross section in a BUU transport model. We then use this
8

(model-dependent) cross section to calculate the viscosity of nuclear matter according to
this same BUU model, thus determining this model-independent quantity in a consistent
manner. Different cross section reductions can give the same amount of stopping, so we
also test whether viscosity is, in fact, a quantity that matters for stopping, or whether the
in-medium cross section matters while the viscosity is not as correlated with the stopping
data.
Several groups (e.g. [18]) treat the in-medium cross section as half the free cross section
[19], which gives approximately the same cross section reduction as effective mass scaling [20–
22]. Gaitanos notes that this resolves some discrepancies with experimental data, but widens
other differences. In the current work, we investigate in-medium cross section reductions that
are dependent on density and/or energy.

2.1

The NN cross section in the nuclear medium

Looking ahead, it is clear that using the nucleon-nucleon (NN) cross section in the BUU
equation (Eq. 1.2) overestimates the amount of stopping found in central collisions at intermediate energy. In this document, 3 approaches are explored. In practice, it is found
that the NN cross section is reduced in the medium. For use in the BUU simulation, the
in-medium NN cross section is applied as a reduced probability of each NN collision to occur.
Thus, each reduced cross section is presented as a reduction factor multiplied by the free
cross section.

9

2.1.1

Tempered

In the Tempered cross section reduction scheme [23], a limiting approach is considered. It
assumes that some reduction happens as the local nucleon density, n, increases, and in the
limit of high density, the resulting cross section radius should not exceed the interparticle
distance. So the resulting cross section σ should follow

σ0 = νn−2/3 .

σ

(2.1)

As this is a very schematic view, a tunable parameter ν ∼ 1 is used. Practically, to smoothly
approach this limit, the following equation is used:

σ = σ0 tanh

σfree
σ0

,

(2.2)

where σ0 is defined in Eq. 2.1 (in principle, other smoothing functions could be used).

2.1.2

Rostock

The Rostock in-medium cross section, developed at Universit¨at Rostock [24], was found using
a thermodynamic T-matrix approach. In their calculation, the reduction was due to Pauli
blocking and self-energy effects [25]. To capture the essence of it without reproducing the
entire method, the resulting reduction is crudely parameterized, independent of isospin:






σ = σfree exp −0.6


ρ/ρ0
1+

10

Tc.m.
150 MeV



,
2


(2.3)

where Tc.m. is the total kinetic energy of the two interacting particles, in the frame where
the local medium is at rest. The Rostock cross section was derived assuming that the total
momentum of the particles is zero, and the velocity of the local medium is also zero, to
simplify the calculations. The cross section is used here in cases where, in the rest frame of
the local medium, the total momentum of the particles is not zero.

2.1.3

Fuchs

Another group [26] notes that in the BUU equation (Eq. 1.2), the in-medium mean field that
the matter is subject to on the l.h.s. should be derived consistent with the in-medium NN
cross section σ used in the collision integral in the r.h.s. The basis for both of these alterations
is the in-medium Dirac-Brueckner T-matrix. This cross section reduction is parameterized
here:






σnn = σnn,free exp −1.7


ρ/ρ0
1+



σnp = σnp,free exp −1.4

Tc.m.
12 MeV




3/2 

(2.4)


ρ/ρ0

.
Tc.m. 
1+
33 MeV

(2.5)

Like the Rostock reduction, this cross section was derived for two particles with total
momentum equal to zero in the rest frame of the local medium. Again, it is assumed that
this can be applied to the case where the total momentum is not equal to zero.

11

2.2
2.2.1

Stopping observables
Linear momentum transfer

In a mass-asymmetric collision of a light projectile colliding with a heavy target, one can
determine the amount of momentum per nucleon that is transferred to the target. As in
the schematic in Fig. 2.1, one finds the momentum of the largest fragment emitted by the
collision. Because of the high mass, that fragment is assumed to stem from the target (the
“target-like residue”). The higher its momentum, the more momentum was transferred from
the projectile. This corresponds to a higher degree of stopping. To compare the observable
across reaction systems, this momentum is divided by the momentum of the center of mass.
Since the velocities involved are non-relativistic, they can be used directly to find the linear
momentum transfer (LMT). This observable was originally used to distinguish between
direct and compound fission reactions in heavy nuclei [27, 28], then used more generally in
nucleus-nucleus collisions [29]. The technique relies on a clear determination of the targetlike fragment, and as the beam energy increases, there are more violent collisions and the
largest fragment produced becomes smaller. Therefore, the practical energy range for this
observable is from around the Coulomb barrier to around 150 MeV/nucleon or so. The
observable LMT is defined as
v

LMT =

vc.m.

,

(2.6)

where v is the velocity of the target-like residue in the beam direction, vc.m. is the velocity
of the reaction center of mass. A higher LMT corresponds to a higher degree of stopping.
Results of LMT for 40 Ar + Cu, + Ag, and + Au are shown in Figs. 2.2, 2.3, and 2.4,
respectively. The experimental data was collected using the K1200 cyclotron at Michigan

12

after collision

Before collision
vbeam

projectile
target
target-like fragment
Figure 2.1: (For interpretation of the references to color in this and all other figures, the
reader is referred to the electronic version of this dissertation.) Schematic of asymmetric
collision. Projectile transfers momentum to the target. To measure linear momentum transfer, the longitudinal velocity of the target-like fragment is compared to the velocity of the
center-of-mass of the collision.
State University’s National Superconducting Cyclotron Laboratory to collide the beam with
targets inside the 4π detector. A Si detector array was used to detect the time-of-flight and
energy of emitted charged heavy fragments (A ≥ 10), yielding the kinetic energy and masses
of those fragments [30]. To determine V in the equation above, a filter on just the heaviest
fragments was used. We assume that these heaviest fragments are a good estimate of the
target-like fragment.
In the BUU calculation, the target-like residue is explicitly tracked throughout the collision, and its velocity is directly calculated from the constituent particles. The various cross
section (CS) reductions described in Section 2.1 are tested, with the results shown with lines
in Figs. 2.2–2.4. The Rostock and Fuchs reductions, as well as the case with no reduction
(“free”) are labeled, while the Tempered CS is marked by its tunable parameter ν.
It is clear from the LMT figures that the free cross section overestimates the stopping in
all three reaction systems, and that the Tempered CS with ν = 0.2 underestimates it. The
Rostock and Fuchs reductions produce LMT values that are very close to each other in all
13

Data (Colin et al.)
1
Ar+Cu

V /Vc.m.

0.8
free
0.6

Fuchs
ν = 0.8

0.4

Rostock
ν = 0.6
ν = 0.4
ν = 0.2

0.2
0.02

0.04 0.06 0.08 0.1 0.12
beam energy [GeV/nucleon]

0.14

Figure 2.2: Linear momentum transfer for 40 Ar+Cu. Lines show the effects of different
in-medium reductions of the NN cross section. The “Tempered” reductions are marked with
their tunable parameter ν. Symbols are experimental data [2].

Data (Colin et al.)
1
Ar+Ag

V /Vc.m.

0.8

Rostock
Fuchs

0.6

free

0.4
ν
ν
ν
ν

0.2
0.02

= 0.8
= 0.6
= 0.4
= 0.2

0.04 0.06 0.08 0.1 0.12
beam energy [GeV/nucleon]

0.14

Figure 2.3: Linear momentum transfer for 40 Ar+Ag. Lines show the effects of different
in-medium reductions of the NN cross section. The “Tempered” reductions are marked with
their tunable parameter ν. Symbols are experimental data [2].

14

1.1

Data (Colin et al.)

1
0.9
Ar+Au
V /Vc.m.

0.8
Fuchs
0.7
0.6

free

0.5
0.4
ν = 0.2
0.3

ν=
0.4

Rostock
ν = 0.8
ν = 0.6

0.2
0.02 0.04 0.06 0.08 0.1 0.12 0.14
beam energy [GeV/nucleon]
Figure 2.4: Linear momentum transfer for 40 Ar+Au. Lines show the effects of different
in-medium reductions of the NN cross section. The “Tempered” reductions are marked with
their tunable parameter ν. Symbols are experimental data [2].

cases, with Fuchs resulting in about 7% higher values than Rostock in the 65 MeV region.
Focusing on the trends, the free CS results in a linear dependence of LMT on beam
energy in all three systems at about 27 MeV and higher, while the reductions all exhibit a
positive concavity in the 40 Ar + Cu and 40 Ar + Ag cases that more closely resembles the data.
In the case of 40 Ar + Au, all calculated lines show a roughly linear beam energy dependence,
while the experimental data shows an even larger concavity compared to the smaller systems.
Judging by eye, the cross section that best fits the 40 Ar + Cu and 40 Ar + Ag data is the
Tempered with ν = 0.4 or 0.6. In the 40 Ar + Au reaction, however, due to the confluence of
effects described in the previous paragraph, the cross section that best fits the data seems
to be the Tempered CS with ν = 0.8. Note that this determination of the in-medium cross
section is only valid within the context of this specific model, the Boltzmann transport

15

simulation described in Section 1.2. The model-independent viscosity is calculated from this
CS in Sec. 2.3.

2.2.2

Rapidity variance ratio varxz

If particles in the hot, dense region of a nuclear collision undergo many collisions (because
the particle cross section is large), the region tends to equilibrate, and particles will lose
memory of which direction they were originally travelling in. Emission from that region will
be isotropic in the reaction center of mass. The observable varxz probes this isotropy. The
FOPI Collaboration defines it thusly [31]:

varxz =

∆yx
,
∆yz

(2.7)

where ∆yx is the variance of particle rapidity1 in the direction that is transverse to the beam
and fixed in the laboratory frame, which is a random transverse direction in the reaction
center-of-mass system. ∆yz is the variance of the longitudinal rapidity (along the beam
direction).
Fig. 2.5 shows the experimental results from FOPI [31] as well as the BUU transport
simulation with the various in-medium NN CS reduction schemes used, for Au + Au, looking
at the distribution of protons. The experiment was performed using the heavy ion accelerator
SIS at GSI/Darmstadt, and charged particles were detected with a good coverage of angles
throughout the 4π region, using the FOPI detector and a set of other detectors that pro1 Rapidity,

yk , is a measure of a particle’s velocity in a certain direction k, defined by
ˆ This is convenient
tanh yk = vk , where vk is the component of particle velocity in direction k.
for relativistic calculations, for if yp/S is the rapidity of a particle p in frame S, and yS/S
is the rapidity of frame S in frame S , then the rapidity of the particle in frame S is given
simply by yp/S = yp/S + yS/S .
16

1.1
free
1

Au+Au

Fuchs
Rostock

varxz

0.9
0.8
ν = 0.8
ν = 0.6
ν = 0.4

0.7
0.6

ν = 0.2
0.1

1
beam energy [GeV/nucleon]

Figure 2.5: Stopping observable varxz for protons in Au + Au collisions at different beam
energies at bred < 0.15. Lines show the effects of different in-medium reductions of the
NN cross section. The “Tempered” reductions are marked with their tunable parameter ν.
Symbols are experimental data from FOPI [3].

vided particle tracking, energy loss determinations, time of flight determination, and charged
particle identification. The beam energies span the range from 0.09 to 1.5 GeV/nucleon, and
the impact parameter selection was limited to bred ≡ b/bmax < 0.15.
Inelastic scattering becomes more important than elastic scattering as the beam energy
increases. (The pion has a mass of 135 MeV/c2 [32], so its production become energetically
allowed in individual NN collisions when Ebeam > 300 MeV/nucleon). Since only the elastic
scattering cross section has been modified for the medium and not the inelastic channels,
comparison for the purposes of determining the in-medium cross sections should be made at
beam energies

400 MeV/nucleon, the energy at which the pion production cross section

starts to become significant in the dynamics. Given that caveat, it seems the CS reduction
that best fits the data below 300 MeV/nucleon is either ν = 0.8 or Rostock.
17

1.1

(a)

(b)

varxz

0.9
0.7
0.5

protons

1.1

deuterons
(c)

(d)

varxz

0.9
0.7
0.5

tritons

Z=1

0.1
1
beam energy [GeV/nucleon]

0.1
1
beam energy [GeV/nucleon]

Figure 2.6: Same as Fig. 2.5, except that particle coalescence is enabled in the simulation,
so comparison to different particle species in experiment is possible.

The observable varxz can also be produced for particles other than protons. The BUU
transport simulation includes composite production that was not enabled for Fig. 2.5 —
that is, no nuclear fragments heavier than A = 1 are produced in the code. If composite
production is activated in the simulation [8], then particles up to A = 3 are produced (i.e.
deuterons, tritons, helium-3). In this case, comparison with experiment for other particle
species are possible and shown in Fig. 2.6.
Looking first at the experimental data, stopping decreases with increasing A, with the
exception of a large enhancement of varxz for tritons at 400 and 600 MeV/nucleon. As

18

the beam energy decreases, especially for tritons, the experimental varxz decreases, while
the simulation shows no such behavior. This is because at low energies in experiment,
the projectile- and target-like residues have secondary decays that are not present in the
transport calculations, making them less accurate. There is an increasing sensitivity to the
NN CS with increasing A. If one focuses on those middle energies, then ν = 0.8 or Rostock
continue to best reproduce the experimental stopping (especially Rostock for the case of
tritons).
Results for varxz in Ca + Ca collisions are shown in Fig. 2.7. Looking at Ebeam
700 MeV/nucleon, one is left with one experimental point that is straddled by theory results
for ν = 0.4 and ν = 0.6. However, the uncertainty of impact parameter in experiment is quite
large (1–2 fm for Au + Au collisions [33]). If the impact parameter gate used in experiment
was twice as much (bred < 0.3 instead of < 0.15), then theory would give the results shown
in Fig. 2.8. In this case, stopping is less sensitive to the CS reduction, but still gives the best
fit for either ν = 0.4 or ν = 0.6.

19

0.8
0.75
0.7

Ca+Ca
bred < 0.15

varxz

0.65
0.6
0.55
0.5
0.45
0.4
1
beam energy [GeV/nucleon]
Figure 2.7: Stopping observable varxz for Ca+Ca collisions at bred < 0.15. Lines and symbols
are the same as Fig. 2.5.

0.8
0.75
0.7

Ca+Ca
bred < 0.3

varxz

0.65
0.6
0.55
0.5
0.45
0.4
1
beam energy [GeV/nucleon]
Figure 2.8: Stopping observable varxz for Ca+Ca collisions at bred < 0.3. Lines and symbols
are the same as Fig. 2.5.

20

2.2.3

Isospin tracer

Another observable that we use as a measure of stopping is based on the isospin tracer
technique, which is used to determine the relative yield of particles from the target and
projectile in a given region of momentum space. This is done by using collisions between
mirror nuclei and comparing to symmetric collisions of each of those two mirror nuclei.
The method is described here:

“The (N/Z)-tracer method is based on the following idea: let us assume that
we are observing the final number of protons, Z in a given cell of the momentum
space. The expected yield Z Ru measured for the Ru + Ru reaction is higher
than Z Zr of the Zr + Zr reaction since Ru has 44 protons as opposed to 40 for Zr.
Such measurements using identical projectile and target deliver calibration values
Z Ru and Z Zr for each observed cell. In the case of a mixed reaction, Ru + Zr
or Zr + Ru, the measured proton yield Z takes values intermediate between the
calibration values (Z Ru , Z Zr ). If, e.g., Z is close to Z Ru in a Ru + Zr reaction,
means that the cell is populated predominantly from nucleons of the Ru projectile
while if it is close to Z Zr it is mostly populated from nucleons of the Zr target. In
this way it is possible to trace back the relative abundance of target to projectile
nucleons contributing to a given cell.” [4]
The method yields the observable Rz , defined as

Rz =

2 × Z − Z Zr − Z Ru
.
Z Zr − Z Ru

(2.8)

In this case, Rz = 1(−1) corresponds to a momentum cell that is populated by protons
completely from the Zr (Ru) nucleus. The case of complete stopping would mean that
the protons completely mix and lose memory of which nucleus they were from. Thus, full
stopping corresponds to Rz = 0.
The experimental results, along with results from the BUU transport model, for collisions
between 96 Zr and one of its mirror nuclei, 96 Ru, are shown in Fig. 2.9, for a beam energy of
21

0.8

FOPI
free
Rostock
Fuchs
ν = 0.8
ν = 0.6
ν = 0.4

0.6
0.4

Rz

0.2
0
-0.2
-0.4
-0.6
-0.8
-1.5

-1

-0.5

0

yz0
Figure 2.9: Isospin tracer observable for central collisions of 96 Zr +96 Ru (bottom) and
96 Ru +96 Zr (top) at 400 MeV/nucleon vs. scaled center-of-mass rapidity, where y = −1
0
is the initial projectile rapidity. Data is from FOPI [4]. A value of Rz closer to zero indicates
a higher degree of mixing, and thus stopping.

22

400 MeV/nucleon and bred < 0.12. The experiment was performed at SIS/ESR-Darmstadt
using the FOPI apparatus, TOF-wall, and the HELITRON and CDC drift chambers [4].
As center-of-mass rapidity gets more negative, the momentum cells are increasingly populated by protons from the target. This makes sense, as the target’s particles are more likely
to backscatter in the center-of-mass (opposite the direction of the beam). As the bins closer
to midrapidity are examined, it is seen that those bins are populated by protons from both
colliding nuclei, as Rz is close to zero there.
It is clear from Fig. 2.9 that use of the free NN cross section overestimates the stopping
(mixing in this case). It seems that the CS best fitting the data here is either Rostock,
Fuchs, or Tempered with ν ∼ 0.8. Statistical uncertainty is a difficulty with this observable,
due to the subtraction of similar values, Z Zr − Z Ru , as well as division by the resulting small
number.

2.2.4

Sensitivity to Equation of State

The equation of state (EOS) of nuclear matter describes how the energy changes with density
and pressure, and is an important part of understanding the nuclear force and the effect on
bulk systems [34]. In particular, the compressibility of nuclear matter, K, is a key input to
calculations of neutron stars and some supernovae. It is commonly defined [35] as
∂ 2 (E/A)
.
K = 9ρ0
∂ρ2
ρ=ρ0

(2.9)

The FOPI Collaboration found dependence of the stopping on the EOS [36], specifically the
compressibility of nuclear matter, using the Isospin Quantum Molecular Dynamics (IQMD)
model [37]. In Fig. 2.10, effects of two choices of compressibility, as well as a momentum-

23

dependent mean field, are shown. ‘S’ refers to a soft compressibility, K = 210 MeV, while
‘H’ refers to a stiff or “hard” compressibility, K = 380 MeV. ‘M’ refers to the inclusion of a
momentum-dependent mean field, achieved by attributing to the nucleon an effective mass
m∗ , here using m∗ /m = 0.7.
The combinations of SM and H give very similar results with our BUU transport model.
The effects range as large as 16% for linear momentum transfer in 40 Ar + Au at high energy,
and as low as within statistical and theoretical uncertainty in the case of varxz in Au + Au
and the isospin tracer Rz . Given a momentum-dependent mean field, the case is certainly
distinguishable between stiff and soft compressibilities in stopping in Au + Au.
For the observables varxz and isospin tracer, the difference between the two plausible
choices, soft with momentum-dependent mean field and stiff with momentum-indepedent
mean field, is smaller than the difference between in-medium cross section reductions that
we investigated. However, LMT is sensitive enough to this choice to render it less reliable
for selecting the best cross section reduction.

2.2.5

Optimal cross section reduction

The determination of the in-medium NN cross section that best fits stopping data is not
especially clear. The results are summarized in Table 2.1. There is a system size dependence,
in that for all the reactions with A < 150, a stronger CS reduction is favored, while in heavier
systems, less reduction in the cross section is needed to reproduce stopping. As speculation,
this could be explained by the fact that in larger systems, each particle undergoes more
chances to collide, since it must travel through a larger interaction volume. This means that
the CS reduction is applied in more instances and has more of an effect on the stopping.

24

Data (Colin et al.)
SM
H
V /Vc.m.

V /Vc.m.

1.1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2

Ar+Ag
free
ν=0.6

1.1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2

0.02 0.04 0.06 0.08 0.1 0.12 0.14
beam energy (GeV/nucleon)

Data (Colin et al.)
SM
H
Ar + Au
free
ν=0.8

0.02 0.04 0.06 0.08 0.1 0.12 0.14
beam energy (GeV/nucleon)

(a) Linear momentum transfer (LMT)
1.2
0.5
free

free

0.8

Rz

varxz

1

ν 0.8

0

0.6
-0.5

Au + Au
0.4

0.8
0.6

ν = 0.8
0.5

free

Rz

varxz

Ca + Ca Exp.
SM
1
H
HM

ν = 0.6

0

-0.5

0.4

SM
H

-1

0.1
1
beam energy (GeV/nucleon)

-1.5

-1 (0)-0.5
y

0

(c) isospin tracer

(b) varxz

Figure 2.10: EOS dependence of stopping observables. S (H) refers to a soft (stiff) compressibility, while M refers to the inclusion of momentum-dependence in the mean field. SM and
H are both reasonable combinations, and they give similar predictions for stopping observables varxz and isospin tracer Rz . However, LMT is as sensitive to these settings as it is to
cross section reductions, so that observable is less reliable for selecting the best reduction.

25

observable reaction system energies [MeV]
best cross section reduction
40
Ar + Cu
LMT
17–115
Tempered w/ 0.4 ≤ ν ≤ 0.6
40
LMT
Ar + Ag
17–115
Tempered w/ 0.4 ≤ ν ≤ 0.6
40 Ar + Au
LMT
27–115
Tempered w/ ν = 0.8
varxz
Au + Au
90–1500
Tempered w/ ν = 0.8 or Rostock
varxz
Ca + Ca
400–1500
Tempered w/ 0.4 ≤ ν ≤ 0.6
96 Zr+96 Ru
Rz
400
Tempered w/ ν = 0.8, Rostock, or Fuchs
(and inverse)
Table 2.1: Summary of cross section determination results. No cross section reduction is
favored universally, though Tempered with ν = 0.6 or 0.8 appears to be the best compromise.

2.3

Viscosity

As a toy model to introduce the concept of viscosity, we describe laminar shear in a classical
system. Consider two plates, with fluid between them, moving in antiparallel directions, in
the steady state. The movement of one plate induces a shear stress, τ , on the layer of fluid
below it, causing that layer to have a velocity u < uplate . That layer induces a shear stress
on the layer under it, and so on. In the linear approximation, these velocities can related
using the equation τ = η(∂u/∂y), where yˆ is perpendicular to the plates. Here, η is the
coefficient of shear viscosity, which is a measure of the momentum transfer in the fluid.

2.3.1

Viscosity from BUU

The shear viscosity coefficient η can be derived from the BUU equation [38]. The result, for
symmetric matter, is reproduced here:

5T
η=
9

2

dp1 f1 p21
dσ p4 sin2 θ
dp1 dp2 dΩf1 f2 f˜1 f˜2 v12 dΩ
12

(2.10)

Note that the NN cross section, σ, enters into the denominator. Since this equation is
derived using the same assumptions as the transport model used to constrain the in-medium

26

140

ν = 0.4
ν = 0.6
ν = 0.8
Fuchs
Rostock
free

η [MeV/fm2 c ]

120
100
80
60
40

ρ/ρ0 = 0.5

20

ρ/ρ0 = 1

ρ/ρ0 = 2

0
10

30
50
70
T (MeV/kB )

10

30
50
70
T (MeV/kB )

10

30
50
70
T (MeV/kB )

Figure 2.11: Viscosity for different temperatures, densities, and CS reductions. Note that
for ρ/ρ0 = 1, the Fuchs reduction overlaps with ν = 0.8 up to about T = 40 MeV/kB , after
which it overlaps with the Rostock reduction.
NN cross section, that cross section can be inserted into this equation to find a viscosity
coefficient η that is model-independent. The calculation was performed with the effective
mass described in Section 2.2.4, which tends to increase the viscosity somewhat.
The results of such calculations are found in Fig. 2.11. In all of the densities and CS
reductions presented, the viscosity blows up at low temperature. At this temperature, the
relative velocity between particles is low. Looking at Eq. 2.10, v12 is in the denominator, so
this behavior makes sense. As the temperature increases, v12 increases while σ first drops
very quickly and then flattens out. The combination of these causes a local minimum of η
to appear. As the density increases, the NN reduction becomes stronger, directly causing an
increase in viscosity, which is observed.
We now use our newly determined viscosity to explore how close nuclear matter is to
being a “perfect fluid”.

27

2.3.2

Lower quantum limit of ratio of viscosity to entropy density

It has been theoretically found that certain strong coupling limits of gauge theories have a
constant ratio of shear viscosity to entropy density regardless of metric used,

η
=
,
s
4πkB

(2.11)

where kB is the Boltzmann constant [14]. Moreover, it is speculated that this value is a
lower limit for all relativistic, finite temperature quantum field theories with zero chemical
potential, and for single-component nonrelativistic gases of particles with spin 0 or 1/2 [14].
At Brookhaven’s Relativistic Heavy Ion Collider (RHIC), in Au + Au collisions with energy
√

sN N = 200 GeV, they found a value for η/s of 0.1 /kB , quite close to the conjectured

lower limit [39] and represented as a filled, red hexagon in Fig. 2.12. We calculate this ratio
at our intermediate energies (IE) studied in this chapter, the entropy density having been
estimated from cluster yields [8]. Represented as open hexagons in Fig. 2.12, the trend of
the nuclear matter looks to match the findings at RHIC, at a much lower energy. The
critical temperature Tc ∼ 170 MeV/kB for nuclear matter was taken from lattice QCD
calculations [40].

28

C

C

Figure 2.12: This work’s estimate for η/s (in open hexagons) with the ratio given for other
materials. RHIC has found a very low value, close to the proposed quantum lower limit
(filled, red hexagon). We find that nuclear matter at intermediate energies approaches this
value. The lines are to guide the eye. Compilation and figure from Ref. [5].

29

2.4

Conclusion

In this chapter, we investigated the viscosity of nuclear matter by adjusting the in-medium
nucleon-nucleon cross section to fit nuclear stopping data with several different stopping
observables across a wide range of nuclear mass and beam energy. We found that an inmedium reduction in the NN cross section is necessary to match experimental data. We
calculated the shear viscosity and found that it is enhanced due to the medium effects,
especially at higher densities. We also calculated the ratio of shear viscosity to entropy
density, η/s, to determine the proximity of our reaction systems to the conjectured lower
quantum limit. We found that the ratio follows a trend that matches with relativistic heavy
ion collisions, suggesting that quark degrees of freedom may not be necessary to approach
the limit of the “perfect liquid” seen at RHIC.

30

Chapter 3
Dynamics

3.1

Introduction

To explore the nature of nuclear matter at supranormal densities through observations,
one must study either neutron stars or heavy ion collisions. Experimentally, we have just
collisions between heavy nuclei, which generate brief pockets of nuclear matter at densities
2–4 times nuclear saturation density. It would be helpful to be able to focus on just the
particles that are emitted from the high-density region. Unfortunately, the spatial and
temporal resolution of current detectors are many orders of magnitude too coarse to directly
trace back the position and time of emission of particles. Instead, for nuclear collisions, the
detectors detect, among other observables, the momentum and energy spectra of emitted
particles. Fortunately, the emitted particles affect each other in transit to the detectors,
if they are close enough and move at similar velocities, developing observable correlations
between the momenta of those particles. These correlations, coupled with information gained
from transport simulations, allow us to identify particles that are emitted specifically from
these high-density regions, giving an opportunity to study how nuclear matter acts at that
scale.
Correlations between particles emitted in nuclear collisions develop in transit from the
source of emission to the detector. Those correlations depend on the particles’ relative
scattering wavefunction, as well as the phase-space distribution of emitted particles. Thus,
31

correlations between particles emitted from collisions of nuclei can be used to determine the
spacetime extent of the source of emission, giving an indication of the size and shape of
the emitting region. This is loosely analogous to using light wave interference to determine
the size and shape of a thin object, like a human hair. If we understand how the wave
interplays with itself, then we can measure the interference pattern with millimeter resolution
to discover the width of the hair at micrometer resolution.
For photons and other noninteracting, identical particles, intensity interferometry utilizes
the Hanbury Brown and Twiss (HBT) effect, illustrated in Fig. 3.1. It was first used to study
radio sources of galaxies by the eponymous authors [41]. Later, it was generalized to include
correlations induced by Coulomb and strong interactions, when it found utility in elementary
particle physics [42], relativistic heavy ion collisions [43], Bose-Einstein condensates [44], and
intermediate energy heavy ion collisions (see Refs. in [45]). For these applications, especially
heavy ion collisions, the generalized technique is now also called femtoscopy.
There are some methods of experimentally determining the timescale of the particles
emitted during a collision, other than the one we show. Generally they involve determining
which types of particle or fragments are emitted first, either using correlation functions of
different species [46] or the N/Z ratio of the fragment [47]. The lifetime of individual emission
sources has also been investigated [48]. For obvious reasons, these do not give the scale of
the time relative to some initial start of the collision.
If particles are emitted from different regions of the collision that are moving with different velocities, then it stands to reason that filtering on these laboratory velocities would
lead to isolating those different sources. For example, one group has identified several different sources of emitted particles from experimental data [49]. These are the source from
the projectile-like fragment, from the target-like fragment, and from the interacting region.
32

They can be identified by their different laboratory velocities. Further, they have identified
two sources in the interacting region: one from the emission as the fragments start to interpenetrate, and again after the participants have equilibrated and possibly formed a neck
region connecting the projectile- and target-like fragments. These have been identified using
a filter on total momentum of the particle pair in the center of mass of the collision. We use
a similar filter to isolate the early emitting source in Section 3.3.
In studying the nuclear collision with BUU, it is important to identify the regime in
which it is valid. In particular, it works fairly well for the liquid phase of nuclear matter,
where the dynamics are affected by the mean field, as well as binary collisions. This is a good
approximation for the early stage of a nuclear collision, when particles are close together,
interact fairly often, and fluctuations are stochastically washed out. At later stages, once
the interacting volume starts expanding, the particles become far enough away from each
other that the dynamics follow the individual particles, the mean field is no longer a good
description, and initial fluctuations will influence the evolution of the system. BUU is not
equipped to handle this situation. Further, fragments of all sizes are produced, while the
particular BUU model only produces fragments up to A = 3. It stands to reason that a
way of filtering the output of experiment to focus on the early stage of the reaction would
be important for accurately comparing to the BUU model.
In Section 3.2, I introduce the HBT interferometry and discussion extraction of the
emission source. In Section 3.3, I will show how filtering the particles can emphasize the
early emission. Finally, in Section 3.4, I show that the model predicts a much stronger
resurgence in emission in heavy systems compared to light systems, explore why this is the
case, and show how this can be detected in an experiment.

33

p1

A

A
+

p2

B

B

Figure 3.1: Schematic of HBT effect. Two particles with momenta p1 and p2 are emitted
from the collision and observed at detectors A and B. If the particles are identical, then the
emission is a quantum superposition of the left and right cases, leading to interference (figure
from Ref. [6]).

3.2

HBT effect

When two photons with momenta p1 and p2 , respectively, are emitted from a nuclear collision
with a small relative momentum q = (p1 − p2 )/2, they will interplay with each other on
the way to the detectors. Even for noninteracting, identical particles, their entanglement is
described by a superposition of the two product states shown in Fig. 3.1. This superposition
leads to a correlation in 2-particle intensity. Hence, the technique of learning about the size
and shape of the emission source is called intensity interferometry.
For protons and other charged particles, another effect causes most of the correlation. If
the protons are moving along with a similar velocity, and if they were emitted close enough
together in time, then they will be subject to elastic scattering between them, either being
repelled from each other by the Coulomb force or attracted by the strong nuclear force.
This will lead to a correlation (or anti-correlation) within their relative momentum, from the
perspective of the detectors. Since both this mechanism and the HBT effect give resolution
on the level of femtometers, the techniques of using the correlations in reactions to learn
about the size and shape of the emission region is called femtoscopy.

34

3.2.1

Correlation function

At the heart of observing the HBT effect is the 2-particle correlation function, which describes the correlation (or anti-correlation) between momentum states, p1 and p2 . Experimentally, it can be observed as the ratio of the 2-particle yield, dN/(dp1 dp2 ), to the
multiplication of two single-particle yields, dN/dp:

C(p1 , p2 ) ∝

dN/(dp1 dp2 )
,
(dN/dp1 )(dN/dp2 )

(3.1)

In nuclear collisions, the event mixing technique is often used, where the 2-particle coincidence yield is taken from a single collision event and the single-particle yields are taken from
different events, to avoid autocorrelation effects. At large relative momenta q = |p1 − p2 | /2,
the emission tends to be uncorrelated and C approaches a constant. It is typically normalized
such that the constant is unity.
To incorporate wider momentum volumes and thereby decrease statistical error, the angleaveraged correlation function is sometimes used. Here, C is integrated over all possible
pair total momenta. The resulting function is dependent on the relative momentum q =
(p1 − p2 )/2:
C(q) ∝

dp dΩ12

dN/(dp1 dp2 )
,
(dN/dp1 )(dN/dp2 )

(3.2)

where p = p1 + p2 and Ω12 is the direction of q.
If the particles are uncorrelated, the probability of detecting a particle with p1 should not
be affected by detecting a particle with p2 . In that case, the probability of detecting both
particles (numerator of Eq. 3.1) is proportional to the multiplication of the probabilities of
detecting each individual particle (denominator) and C(q) = 1.

35

Figure 3.2: Sample correlation function [7]. At large q, particles are not correlated and the
correlation function flattens to unity.
An example of a p-p correlation function is shown in Fig. 3.2. Below q = 20 MeV/c, it
is less than one, which represents anticorrelation as the Coulomb force and Pauli blocking
deplete low relative momenta. At large q, the protons are uncorrelated, since the probability
for them to interact is so low, and thus C(q) tends to unity. In the middle, C > 1 shows
correlation, here caused by the 1 S0 strong nuclear attraction.
In order to quantitatively analyze the correlation function and resolve the size and shape
of the emitting source, we need to use the scattering wavefunction between the two particles
and the Koonin-Pratt equation. This is described in the next section.

3.2.2

Source function

In experiment, the correlation function can be directly calculated from observables using
Eq. 3.2. We would like to use it to learn about the size and shape of the region that emitted
these particles. To resolve the actual size and shape of the emission source, the Koonin-Pratt

36

equation [50, 51] is used:
C(q) =

2

dr Ψq (r) S(r)

(3.3)

The source of correlated particles is characterized by the source function, S(r), which is
defined as the probability of two particles to be emitted at a relative distance r in their
center of mass. Ψq (r) is the scattering wave function of the two particles involved. This
is analogous to the Young double-slit experiment, where the intensity of light on a screen
(correlation function) is related to the superposition of the waves (wave function) and the
width and orientation of the slits (source function).
The corresponding angle-averaged source function S(r) is found by inverting the equation

C(q) − 1 = 4π

where K(r, q) is the scattering kernel Ψq (r)

dr r2 K(r, q) S(r),

2

(3.4)

− 1, averaged over all angles.

To find the source function, an Ansatz is used, generally consisting of a Gaussian function
2
2
e−r /(2R ) . The width R is adjusted to fit the experimental source function. This typical

Ansatz assumes that the shape of the source is Gaussian, and that there is only one source of
emission during the collision. As noted in Section 1.1, there are several different dynamical
stages of the nuclear collision, which implies that particles will be emitted in different patterns
in different dynamical stages, and thus there are several different emission sources. This
makes interpretation of the Gaussian width more difficult.
Instead, the source function can be found directly by projecting Eq. 3.3 onto basis splines
and inverting the resulting matrix. This allows the determination of the source with less prior
assumptions about its shape [52].
In BUU, the source function can be directly calculated from the relative positions and
37

momenta of emitted particles, at the time of their emission. This is advantageous, as one
need not rely on a theoretical wave function Ψq (r) to construct the source. Since correlation
is only built up between particles with similar momentum, the emitted particles are sorted
into cubic momentum bins that are 25 MeV/c wide in each direction (in the laboratory
frame). Within each bin, the distance between them, r = |r1 − r2 |, is calculated, and this
particle pair is counted as contributing to that part of the source function S(r).
To determine whether a test particle is “emitted” in our BUU model, each particle is
marked with the time of its last collision. At the end of the time evolution, if the density
of the local medium is less than 1/15th nuclear saturation density n0 = 0.16 fm−3 , then
the test particle is considered emitted (if it is still in a higher-density region, then it would
probably collide again, given more time).

3.3

Timescale of emission

To explore the dynamics of collisions, we choose one collision between high-mass nuclei,
and one between low-mass nuclei, 112 Sn +112 Sn and 40 Ar +45 Sc. Both have the same beam
energy (50 MeV/nucleon) and impact parameter selection (breduced < 0.46), as well as similar ratio of protons to neutrons (isospin asymmetry). In order to exclude emission from
evaporation of the residues and isolate the high-density region, an angular filter is used to
only accept particles that are emitted transverse to the beam direction in the center of mass
(70° ≤ θcm ≤ 110°) [53].
The source functions S(r) for the two systems, calculated from the transport theory, are
seen in Fig. 3.3, with a zoomed view of the long-range tail inset in log scale. Since S(r) is
a probability density, its 3D integral sums to unity, and thus the height at a given distance

38

0.45

40 Ar +45 Sc
112 Sn +112 Sn

0.4
10−1

S(r) [10−3 fm−3 ]

0.35
0.3
0.25
0.2

10−2

0.15
0.1

6

8

10 12 14 16 18 20

0.05
0
0

5

10
r (fm)

15

20

Figure 3.3: Source functions from BUU transport code. Inset is a close-up in log scale of
the long-range part. The source for the smaller system is more peaked at short range, and
it falls off faster at long range, compared to the larger system.
shows the relative weight of different ranges of the relative distance between 2 particles r.
In the simulation, the integral of S(r) is normalized in the range from 0 to 50 fm. Looking
at the two sources, 40 Ar +45 Sc shows a greater relative weight of the short-range component
of the source, even as 112 Sn +112 Sn has a greater number of particle pairs emitted with that
relative distance, as witnessed by the smaller statistical error bars. At long range, visible
in the inset, the source function of the smaller system falls off more rapidly, though both
exhibit roughly exponential behavior.
Next, we describe a more systematic way of characterizing the shape of the source, in
order to use it to understand the collision dynamics better.

39

3.3.1

Interpreting the source function

We want to have a quantitative method to compare source functions with each other, so
that when we apply filters to emitted particles and produce the resulting source functions,
we are able to more precisely notice any interesting effects. Additionally, experimental correlations are often fitted using analytic source functions, rather than inversion as described
in Section 3.2.2. Therefore, it can be helpful to analyze the source function produced by the
transport model in terms of such functions. To do this, we investigate different parameterizations that are fit to the source functions.
Different groups have come up with different ways of quantifying the size or shape of the
source. If all that is wanted is the source size, then r1/2 , the half width at half maximum,
could be used. This has the advantage that it is widely used and gives an estimate independent of any parametrization used. However, it loses details, especially in the behavior of the
source at longer ranges. For this, one can use various fits to the source. Some have used
generic cubic splines [54], 2-gaussian “hump fits” [55], or a superposition of Gaussian and
decaying exponential behavior. Since this last one is the behavior of the source functions in
Fig. 3.3, this is what we use1 .
However, a decaying exponential has a nonzero derivative of the form e−r at r = 0. Since
this is a one-dimensional function projected from three dimensions, the function should be
smooth across r = 0. To compensate, the following form is used [56]:

S(r) =

λexp
λgau
r2
exp − 2
+
exp −
Ngau
Nexp
2Rgau

r2 + β 2
Rexp

(3.5)

1 We also tried the convolution of two Gaussian sources with different widths and strengths,

but it did not match the simulated behavior.

40

collision

λgau

λexp

Rgau (fm) Rexp (fm)

fexp

40 Ar +45 Sc

0.100(5) 0.915(4)

3.47(5)

5.09(5)

0.783(8)

112 Sn +112 Sn

0.006(3) 1.046(4)

4.16(15)

6.18(3)

0.87(2)

Table 3.1: Fit parameters for source functions. λgau and λexp are the fractions of particle
pairs that contribute to the first and second terms of the equation. Rgau and Rexp are the
widths of the Gaussian and exponential components, respectively. The fraction of particle
pairs that contribute to the Gaussian shape, fgau , is higher for the smaller system, while
the larger system has a correspondingly longer exponential tail. Uncertainties are statistical
from fit and do not include uncertainties of underlying model.
2 /R
where β = Rgau
exp , and Ngau and Nexp are chosen to make each term have a 3D

integration equal to unity when ignoring the λs (these factors are given in Eq. B.3). Thus,
since S(r) is a probability and should integrate to unity, λgau and λexp are the fractions of
particle pairs that contribute to each term and should add to unity, if the equation accurately
represents the source. At short distances, the second term acts as a Gaussian with the
same radius as the first term, while at long distances, it behaves as e−r , thus enabling the
parameterization of the source into two components: Gaussian and exponential decay.
To find out the relative weight of the Gaussian and exponential components, it is insufficient to use the λ parameters above, since the λexp term contains part of the Gaussian
component. Instead, to find the fraction of proton pairs that contribute to the Gaussian
part of the source, the square root in the second term is expanded to the 2nd order in r, and
the resulting Gaussian is integrated over the sphere to get the exponential fraction fexp (see
Eq. B.3 for the derivation).
The results of the fitting procedure (using MINUIT [57]) for the reaction systems studied
are shown in Table 3.1. The resulting curves are shown along with the simulation results in
Fig. 3.4. The Gaussian fraction fgau is also listed in Table 3.1.
Looking at the source function, the smaller system has a smaller Gaussian width and a
41

0.45

40 Ar +45 Sc
112 Sn +112 Sn

0.4
10−1

0.35
S(r) [10−3 fm−3 ]

pBUU
pBUU

0.3
0.25
0.2
10−2

0.15
0.1

6

8

10 12 14 16 18 20

0.05
0
0

5

10
r (fm)

15

20

Figure 3.4: Source function S(r) from BUU transport simulation for 112 Sn +112 Sn (triangles)
and 40 Ar +45 Sc (squares). Also shown are curve fits using a short-range Gaussian and longrange exponential (seen in the logscale plot inset).
slightly smaller exponential tail, as expected. Overall, fexp is smaller for 40 Ar +45 Sc than
112 Sn +112 Sn.

Next, we use our time dependent model to investigate when and where these

different components are emitted, and see if we can enhance one versus the other and thereby
isolate one emission source for future study.

3.3.2

Time dependence of the source function

To determine when the protons are emitted that contribute to the various parts of the source,
we use the BUU transport model described in Section 1.2 to simulate the collisions. Then,
we use the explicit timing information given in the simulation to filter on different times and
look at the result.
42

Naively, if the source function is constructed from only particles at early times, then the
source should be smaller, since particles have less time to travel away from each other before
being emitted. This suggests that particles from early time are correlated strongly with the
short-range part of the source. Conversely, particles that are emitted from a source that has
been expanding for some time will be, on average, farther apart from one another, leading
to late emitted particles contributing primarily to longer ranges of the source function. A
caveat here is that, in the calculation of the source function from the simulation, all particles
are paired with others with similar momentum, regardless of time, so the early particles can
contribute to longer ranges by being paired with a particle from a later time.
Using the transport code, we can test this by constructing a source function from particles
that are emitted at selected times. This is done in Fig. 3.5 for 40 Ar +45 Sc. The early-time
function matches the rationale given in the previous paragraph.

43

1.8
no time filters
t < 100 fm/c
t > 100 fm/c

1.6
40 Ar +45 Sc

S(r) [10−3 fm−3 ]

1.4

100

1.2

10−1

1

10−2

0.8

10−3

0.6

10−4
10−5

0.4

6

8 10 12 14 16 18 20 22

10
r (fm)

15

0.2
0
0

5

20

Figure 3.5: Source function for 40 Ar +45 Sc constructed from only early time (green, long
dash) and only late time (blue, short dash).

44

3.4

Late-time resurgent emission

Now that we have a qualitative understanding, using a transport code, of how the emission
time affects the source function, we examine the history of particle emission directly, from
our transport model. Such a timeline of emission of protons is shown in Fig. 3.6. In both
systems, there is one early peak, and in 112 Sn +112 Sn, one broad hump is seen at later times.
This resurge in emission has not previously been predicted, to the author’s knowedge. If it
were evaporation (deexcitation of excited fragments), it should start immediately after the
initial interpenetration. Instead, there is a lull in emission followed by a burst. To learn
more, we inspect the phase distribution of these particles to find out what is happening.

3.4.1

Analysis of second emission

To start, we look at the spatial distribution of nucleons during the collision. Fig. 3.7 shows
the nucleon density for both collision systems. Initially they have a similar evolution: interpenetrating spheres lead to a large density in the middle of the collision. The collision of the
heavier system takes longer, since the nuclei are larger in volume. Around ∼ 80 fm/c, the
spectator regions have passed each other and continue on their merry way. A dilute region
of nuclei are left in the middle, the so-called “neck” region (seen around 120 fm/c). At later
times, around 230 fm/c, a qualitative difference emerges — in 112 Sn +112 Sn, the particles in
the neck region draw close together to form two sizable high-density regions, while in the
lighter system, a single recontracted higher-density region is barely visible in the middle of
the neck. Comparing with Fig. 3.6, this recontraction happens at the same time that proton
emission increases.
To confirm that the protons are being emitted from these recontracting fragments, one

45

0.08

40 Ar +45 Sc

0.07

112 Sn +112 Sn

dN /dt [c/fm]

0.06
0.05
0.04
0.03
0.02
0.01
0
0

50

100

150
200
time [fm/c]

250

300

Figure 3.6: Proton emission versus time, from BUU transport simulation. Note the clear
separation between early and late emissions and rise in late-time emission in the heavier
system.
can see where the emitted particles are coming from, spread out over time. Fig. 3.8 shows
this. Indeed, the emission is coming from the higher density region at early time, and at
later times, an increase in emission occurs at the same time and location as the recontraction
of the 2 fragments in the neck region.
To determine what causes this recontraction and emission, we compare the two collisions we are studying. The most significant difference between the two (for central collision
dynamics, at least) is the volume of the interaction region and number of nucleons participating. One effect is that there is more volume in which NN collisions occur in the larger-A
collision. Perhaps the nucleons collide so much that enough of them are pulled out of the
projectile- and target-like fragments and into the neck region, where the mean field then
causes them to seek higher densities (closer to the nuclear saturation density) and recontract
into fragments.
To test this theory, the nucleon-nucleon cross section (σNN ) in the heavier system can
46

-20-10 0 10 20
z [fm]

40.4 fm/c

120.4 fm/c

230.4 fm/c

-20-10 0 10 20
z [fm]

-20-10 0 10 20
z [fm]

-20-10 0 10 20
z [fm]

ρ/ρ0 [fm−3 ]

t = 5.4 fm/c

2
1.6
1.2
0.8
0.4
0

2
1.6
1.2
0.8
0.4
0

ρ/ρ0 [fm−3 ]

x [fm]

20
10
0
-10
-20

x [fm]

(a) 40 Ar +45 Sc
20
10
0
-10
-20

t = 8.4 fm/c

-20-10 0 10 20
z [fm]

43.4 fm/c

128.4 fm/c

238.4 fm/c

-20-10 0 10 20
z [fm]

-20-10 0 10 20
z [fm]

-20-10 0 10 20
z [fm]

(b) 112 Sn +112 Sn
Figure 3.7: Density contours at selected times after the beginning of collision of (a)
40 Ar +45 Sc and (b) 112 Sn +112 Sn. The heavier system has a recontraction of density at
later times. Time is normalized so that the centers of the nuclei cross at 51 fm/c (as predicted with Coulomb effects, not including stopping).

47

0

50 100 150 200 250
t [fm/c]

40
30
20
10
0
-10
-20
-30
-40

112 Sn +112 Sn

0

(a)

45
40
35
30
25
20
15
10
5
0

dN /(dzdt) [c/fm2 ]

z [fm]

40 Ar +45 Sc

45
40
35
30
25
20
15
10
5
0

dN /(dzdt) [c/fm2 ]

z [fm]

40
30
20
10
0
-10
-20
-30
-40

50 100 150 200 250
t [fm/c]
(b)

Figure 3.8: 2D histogram that shows where along the beam axis (z) protons are emitted as
time goes on, for a) 40 Ar +45 Sc and b) 112 Sn +112 Sn. In the Sn collision, a second emission
is quite pronounced, starting at about 125 fm/c. At late times, the emission is seen to come
from 4 distinct regions, which are the positions of high density in Fig. 3.7b.
be reduced artificially, to mimic the effects of fewer collisions. To get a qualitative look, we
reduce it by 5 times. The reduction causes the resurgent particle emission to vanish, as seen
in the timeline of emission in Fig. 3.9, while not significantly altering the early emission. In
Fig. 3.10, we see that there is no recontraction of particles in the neck region at later times,
and in Fig. 3.11, the source function shows a slight reduction in the long-range exponential
part, which may be consistent with the behavior of the source function of the 40 Ar +45 Sc
collision. Thus, we learn that the formation of recontracted fragments in the neck region
is due to the certain ratio of nucleon mean free path to the size of the interacting region2 .
Now we show how to detect whether this recontraction occurs in reality, since we cannot
2 Production

of composite particles (deuterons, tritons, and 3 He) was turned off for these
simulations, as the model used for composite production [8] is derived for higher energies
and might not be valid at this energy [58]. Turning it on changes some of the dynamics,
including causing some recontraction in the lighter system. It still shows and even reaffirms
the effect of late-time recontraction causing emission of particles.
48

0.08

(0)

σNN = σNN
(0)
σNN = σNN /5

0.07
dN /dt [c/fm]

0.06

112 Sn+112 Sn

0.05
0.04
0.03
0.02
0.01
0
0

50

100

150
200
time [fm/c]

250

300

Figure 3.9: Emission timeline for 112 Sn +112 112 compared to the same with an ad hoc
reduction in the NN cross section. The reduction suppresses the late emission, while leaving
the early emission mostly unaffected.
manually adjust this ratio in experiment.

49

t = 8.7 fm/c

-20-10 0 10 20
z [fm]

43.7 fm/c

128.7 fm/c

238.7 fm/c

-20-10 0 10 20
z [fm]

-20-10 0 10 20
z [fm]

-20-10 0 10 20
z [fm]

2
1.6
1.2
0.8
0.4
0

ρ/ρ0 [fm−3 ]

x [fm]

20
10
0
-10
-20

Figure 3.10: Density contours at selected times after the beginning of collision of
112 Sn +112 Sn with an ad hoc reduction of σ
NN by a factor of 5. Unlike the unreduced
case seen in Fig. 3.7b, there is no formation of dense fragments from the neck region at late
times.

0.0004
BUU
BUU, σNN /5

0.00035
0.0003

S(r)

(fm−3 )

0.0001
0.00025
0.0002

1e-05

0.00015
0

0.0001

5

10

15

20

5e-05
0
0

5

10
r (fm)

15

20

Figure 3.11: Source function for 112 Sn +112 Sn calculated in a transport model with a nominal
NN cross section (triangles) and with an ad hoc reduction in that cross section (squares).
The lines represent a fit using the Chung parameterization.

50

3.4.2

Detection in experiment

To determine whether this recontraction occurs in experiment, we need an observable that
is based on the momentum structure, since that is all we have access to, and sensitive to the
shape of the emitting source. As mentioned in the introduction, Ghetti and collaborators
have used filters on laboratory velocity and total cms momentum of the pair to focus on
sources of emission coming from different spacetime regions of the collision. Thus, it stands
to reason to look at the spectrum of particles in transverse (to the beam) momentum. Such a
2D histogram is shown in Figs. 3.12 and 3.13 for 40 Ar +45 Sc and 112 Sn +112 Sn, respectively.
The figures show how many protons are emitted at a given cell of transverse momentum
(p⊥ /m) and time (t). The two emission regions at early and late time are clearly seen. In
112 Sn +112 Sn,

the second re-emission is visible, showing that it only occurs at low p⊥ /m

(below ∼ 0.2 c), while the early emission extends to ∼ 0.4 c for both systems.3 We learn
that for both systems, a filter on high p⊥ /m suppresses the later emission, emphasizing the
protons emitted during the initial interpenetration. Thus, we correlate a theoretical quantity,
time, with an experimental observable, transverse momentum — by filtering on high p⊥ , we
can access the particle sample dominated by the early emitted protons.
The effects of this filter in the source function S(r) are seen in Figs. 3.14 and 3.15
for 40 Ar +45 Sc and 112 Sn +112 Sn, respectively. Also included are best-fit lines using the
Chung fit parameterization, Eq. 3.3.1, including a short-range Gaussian and a longer-range
exponential. In both cases, the relative proportion of protons contributing to the short-range
Gaussian component is increased, as reflected by an increased value (since the source function
3A

feature for future study is the 2 distinct regions at early time in 112 Sn +112 Sn, 1 at
low p⊥ /m and 1 at high. Turning composite production on shifts the lower one later in time
relative to the more energetic region

51

dN/dt [c/fm]

d2 N
[fm−1 ]
d(p⊥ /m)dt
0.075

0.3
0.2
0.1
0

0.05
0.025
0

p⊥ /m [c]

0.4
0.3
0.2
0.1
0
0

50

100

150
t [fm/c]

200

10 20
dN
[c−1 ]
d(p⊥ /m)

250

0

Figure 3.12: 2D Histogram of p⊥ /m vs. t for protons emitted from collisions of 40 Ar +45 Sc.
Also displayed are 1D histograms projected onto both axes. Emphasis on early emission in
an experiment is possible with a filter on p⊥ /m 0.15 c.

52

dN/dt [c/fm]

d2 N
[fm−1 ]
d(p⊥ /m)dt
0.4
0.3
0.2
0.1
0

0.08
0.06
0.04
0.02
0

p⊥ /m [c]

0.4
0.3
0.2
0.1
0
0

50

100

150
t [fm/c]

200

0 15 30 45 60
dN
[c−1 ]
d(p⊥ /m)

250

Figure 3.13: Same as Fig. 3.12, but for 112 Sn +112 Sn. The second surge in emission is very
visible at late times.

53

0.003
no filters
p⊥ /m > 0.15c

40 Ar+45 Sc

0.0025

0.01
0.001

S(r) (fm−3 )

0.002

0.0001
1e-05

0.0015

1e-06
0.001

1e-07
0

5

10

15

20

0.0005
0
0

5

10
r (fm)

15

20

Figure 3.14: Source function for 40 Ar +45 Sc. Inset employs logscale to illuminate the effects
of the high p⊥ gate at far-away distances.
is a probability distribution). To more precisely see this effect, the best-fit Gaussian width
and exponential width are given in Figs. 3.16 and 3.17.
In other work, it was found that the short-range size of the source, characterized by r1/2 ,
the half width at half maximum, was different by about 0.2 fm, depending on the momentum
filter used [59]. Looking more closely into the shape of the source using the Chung fit, we
see that the width of the Gaussian source, Rgau , changes very slightly (The error bars only
reflect the uncertainty and goodness of the fit, not the model uncertainties). Meanwhile,
the width, or shallowness, of the exponential tail is affected greatly, suggesting that this
region was primarily populated by protons with low p⊥ /m. This notion is verified by the
calculation of the fraction of protons contributing to the exponential component, fexp , also
given in Figs. 3.16 and 3.17 — when we filter on high transverse momentum, we greatly

54

0.0003
no filters
p⊥ /m > 0.15c

0.00025
0.0001

S(r) (fm−3 )

0.0002
0.00015

1e-05

0.0001
0

5

10

15

20

5e-05
0
0

5

10
r (fm)

15

20

Figure 3.15: Source function for 112 Sn +112 Sn. Inset is in logscale, to see the effects of the
high p⊥ gate at long range.

55

9
8

no cuts
p⊥ /m > 0.15c
t < 100 fm/c

1

7

0.8

0.6

5
4

0.4

fraction

width [fm]

6

3
0.2

2
1

0

0
Rgau

Rexp

fexp

Figure 3.16: Fit parameters for source functions created with filters on the emitted protons,
for 40 Ar +45 Sc. The filter allowing only early t (open squares) results in similar parameter
values as for high p⊥ /m (open triangles), showing that filtering on high p⊥ /m selects for
early t. Further, these filters enhance the Gaussian fraction fgau , showing that the Gaussian
component is primarily at early time.
decrease the contribution to the exponential component.
To verify that this is really representing the early part of the collision, we also construct
the source function including only protons emitted before 100 fm/c and show the fit parameters as well. For all the parameters but Rexp for 112 Sn +112 Sn, it shows the same trend as for
high-p⊥ /m, indicating that this momentum filter indeed does allow a dependence on early
time. For Rexp , they have diverging behaviors. However, the uncertainty for this width for
the early time filter shows that the there is still suppression of the exponential component,
as there is less of an exponential trend to fit precisely.
According to Figs. 3.16 and 3.17, the source functions of both collisions under consider-

56

9
8

no cuts
p⊥ /m > 0.15c
t < 100 fm/c

1

7

0.8

0.6

5
4

0.4

fraction

width [fm]

6

3
0.2

2
1

0

0
Rgau

Rexp

fexp

Figure 3.17: Same as Fig. 3.16, but for 112 Sn +112 Sn.
ation are affected qualitatively similarly by the filter p⊥ /m > 0.15c. Since 40 Ar +45 Sc does
not have any late-time recontraction and 112 Sn +112 Sn does, then this behavior of the source
function is not a good indicator of this recontraction that we can search for in an experiment.
Instead, we look at the collision in which we made an ad hoc reduction in the NN cross
section in Section 3.4. The effect of the p⊥ /m filter on the source function is shown in
Fig. 3.18. The width of the Gaussian and exponential components of the source, Rgau and
Rexp respectively, are quite similar to that of the calculation with a nominal cross section, as
is the effect of the p⊥ /m filter on Rexp . Thus, the shape of each source is about the same as in
the nominal case, which gives credence to using this calculation as a prediction of dynamical
behavior, even though it has an unphysical NN cross section. The major difference is in
the fraction of proton pairs that contribute to the long-range exponential component, fexp ,
with the calculation with the reduced NN cross section having a substantially suppressed
57

width [fm]

8

1

no cuts
p⊥ cut
σ0 /5
σ0 /5, p⊥ cut

0.9

7

0.8

6

0.7

5

0.6

4

0.5

3

fraction

9

0.4
Rgau

Rexp

fexp

Figure 3.18: Parameters describing the shape of the source function for 112 Sn +112 Sn, calculated with a filter on transverse momentum, p⊥ /m > 0.15c (squares) and without (triangles),
as well as having the simulated reaction calculated with a nominal NN cross section (open
shapes) and with an ad hoc reduction in that cross section, σ = σ0 /5 (filled shapes), which
leads to no late-time recontraction. In the σ/5 case, the p⊥ cut does not affect the exponential fraction fexp like it does in the nominal case. Therefore, the sensitivity of fexp to a cut
in p⊥ is a signal of this recontraction.
exponential fraction. This is consistent with our earlier conclusion that the substantial
exponential component is caused by the late-time recontraction, which is not seen in the
time evolution in this calculation.
The summary of fexp results is shown in Fig. 3.19. In the simulation with the ad hoc
σNN reduction, this fraction is not affected by the p⊥ /m filter, as it is in the nominal
σ case. This suggests that in the nominal case, there are two emission sources and one is
suppressed by the filter more than the other, where in the σ/5 case, there is one source that is
suppressed universally. Thus, if in experiment, the exponential fraction in the source function

58

1

no cuts
p⊥ /m > 0.15c

0.8

fexp

0.6
0.4
0.2
0

112 Sn+112 Sn

40 Ar+45 Sc

112 Sn+112 Sn,

σNN /5

Figure 3.19: Fraction of proton pairs, fexp , that contribute to the long-range exponential
component of the source function, for the heavy system (left), light system (middle), and
heavy system with ad hoc reduction in σNN (right). The first two collision systems produce a
recontraction in the neck region and resurgent emission. Here, a filter on high p⊥ /m causes a
preferential suppression of the exponential component. With reduced σNN , the heavy system
exhibits no re-emission, and the filter does not noticeably change fexp , suggesting that the
filter is only acting on 1 emission source, rather than 2.
is significantly suppressed by a filter on high p⊥ /m, then there is probably a recontraction
of nucleons occurring after the initial compression and expansion stages of the reaction.

3.5

Conclusion

In this chapter, we investigated the dynamics of central nuclear collisions at intermediate
energies. In particular, we used the source function to learn about the size and shape of the
proton emission region and provided an observable that correlates with late-time dynamics
of the collision.
We learned that the formation of recontracted fragments in the neck region is connected
with a certain ratio of nucleon mean free path to the size of the interacting region. To study

59

this effect, we demonstrated the correlation betweeen a theoretically accessible quantity,
time, with an experimental observable, transverse momentum.
Finally, we predicted a second emission source, coming from a recontraction of particles
in the neck region at late times. We showed that with a transverse momentum filter, we
see this late-time source being affected differently than the early-time source, as evidenced
by the exponential fraction changing significantly with the application of the filter. This
is unlike many analyses of experiments at intermediate energies, where, in the end, only
one source size was seen, no matter the transverse momentum filter. We showed that to
test the theory of a late-time recontraction, one can use a high-p⊥ filter and compare the
relative contributions to the short- and long-range parts of the source function. If the relative
contribution is the same, then only one source was present; if it changes, then two sources
are present in the collision.

60

Chapter 4
Towards a fully quantal approach
In this work so far, we have simulated central nuclear collisions with a BUU transport model.
While it has been and continues to be very useful and retains much predictive power, it has
some drawbacks that we would like to overcome. The approach is inherently semi-classical,
with an arbitrary divide between mean field and collision effects, and it is disconnected from
approaches to nuclear structure. In this chapter, we take some early steps in implementing
a transport model that does not have the above drawbacks.
Nonequilibrium Green’s function (NGF) techniques [58, 60–62] represent a powerful tool
to describe the evolution of correlated quantum many-body systems. The study of the
Kadanoff-Baym (KB) dynamics in different quantum many-body systems is an ongoing
challenge for theorists. The first direct implementations of NGF in homogeneous nuclear
systems saw light almost 30 years ago [58]. Later advances include, among others, the
extension to homogeneous electronic systems [63], small atoms [64], small molecules [65],
quantum dots [66], and, more recently, double excitations in electronic systems [67, 68]. The
advancement of dynamical quantum many-body methods is necessary for our understanding
of systems that range from ultracold atoms [69] to nuclei [70] and atoms [64], a range of
phenomena spanning 15 orders of magnitude in energy (neV to MeV).
Within nuclear physics, NGF techniques have been mostly used for derivations rather
than exploited directly [62, 71, 72]. With so few quantitative studies on the impact of correlations in the time evolution of nuclear systems, one might wonder about the validity of
61

nuclear reactions simulations that are based on the mean-field picture or on semiclassical
approaches. To clarify these issues, we would like to develop nuclear reaction simulations
beyond the Hartree-Fock approximation. Nonequilibrium Green’s functions, their time evolution described by the KB equations, represent a large step in that direction. However, with
a naive algorithm, the KB equations are too computationally costly; therefore, we present
a method to reduce this cost. We then apply a simplified model to the problem of 1D Bose
Einstein Condensates (BEC), as they are an example of a physical quantum system in one
dimension.
Historically, just a handful of methods have been developed to describe central nuclear
reactions. Few of these are general enough to be employed in generating practical predictions. The time-dependent Hartree-Fock (TDHF) method has been exhaustively employed
in describing low-energy reactions [73]. Nowadays, TDHF simulations can be performed in
full 3D and involve nuclei as heavy as uranium [74]. However, the validity of TDHF requires
correlations to play a negligible role in the dynamics [72]. At low energies, one expects that
the role of correlations is minimal due to antisymmetrization effects. Conversely, one would
expect correlations to dominate at higher beam energies, where the Pauli principle plays less
of a role in the reaction dynamics.
In the semiclassical limit, the two-time KB equations reduce to a single-time Boltzmann
equation, like the one used in the previous chapters. The Boltzmann equation approach
deals with the increasing complexity of the reaction by using simplifying assumptions with
an average, single-particle approach. Because of its inherently semiclassical nature, however,
this kind of description remains genuinely disconnected from the quantum methods employed
for nuclear structure [75–77]. More importantly, there is no systematic way of extending
Boltzmann equation-type approaches to include more quantal effects.
62

Correlations can lead to a fast thermalization of the occupation of single-particle states
and to enhanced stopping, compared to what can be found in TDHF [78–80]. Further,
quantum effects may dominate in high-density regions and thus affect reaction observables.
Therefore, if one is aiming at increasing incident energies, it is important to develop a
quantal approach to central nuclear reactions that, besides the mean field effects, can also
incorporate correlations. Unfortunately, even at the single-particle level, NGF techniques
require handling vast amounts of information that can easily overwhelm the capabilities of
computing systems, rendering this approach impractical in the naive approach.
In this chapter, we first introduce briefly the Kadanoff-Baym equations, which describe
the time dependence of nonequilibrium Green’s functions, followed by the mean field approximation that we use in the rest of the chapter. Next, in Section 4.3, we implement the mean
field approximation in 1D, describing the time evolution of a collision of nuclear “slabs”,
including the solution of the initial ground state using the same framework as the time evolution and a discussion of the off-diagonal elements of the density matrix. In Section 4.4, we
test the importance of those far off-diagonal elements, in anticipation of truncating them to
reduce the computational cost.
Next we describe a coordinate rotation of the density matrix that allows us both to partially decouple the position and momentum domain, increasing the efficiency of accessing
arbitrary momentum regimes, as well as straightforwardly truncating the density matrix in
the off-diagonal direction to reduce the computational cost. Finally, we implement the density matrix evolution to describe the initial state of a Bose-Einstein Condensate of ultracold
atoms.

63

4.1

Kadanoff-Baym equations

We aim to describe the dynamics of correlated nuclear systems using the KB equations.
They are, in general, derived without any particular assumption on the physical system under
consideration. They describe the time evolution of the single-particle Green’s functions for
fermions or bosons,

G< (x1 , t1 ; x1 , t1 ) ≡

ˆ(x1 , t1 ) ,
i a
ˆ† (x1 , t1 ) a

(4.1)

G> (x1 , t1 ; x1 , t1 ) ≡ −i a
ˆ(x1 , t1 ) a
ˆ† (x1 , t1 ) ,
where a
ˆ† and a
ˆ are the single-particle creation and annihilation operators, respectively, and
where the expectation values, · ≡ Tr{ˆ
ρi ·}, are taken with respect to an initial density
matrix, ρˆi , at t = t0 . We consider initial states that are uncorrelated, even though the
theory can incorporate correlated initial states [58].
The KB equations for 1D systems, governing the evolution of Green’s functions in their
arguments, follow from considerations of the equations of motion for creation and destruction
operators [60]:

i

2 ∂2
∂
+
∂t1 2m ∂x21

G≶ (11 ) =
+

dx2 ΣHF (x1 , x2 ; t1 ) G≶ (x2 , t1 ; 1 )
t1

d2 Σ+ (12) G≶ (21 ) +

t0

−i

2 ∂2
∂
+
∂t1
2m ∂x2

G≶ (11 ) =

t

1

d2 Σ≶ (21 ) G− (21 ) ,

t0

dx2 G≶ (1; x2 , t1 ) ΣHF (x2 , x1 , t1 )

1

+

t1

d2 G+ (12) Σ≶ (21

t0

t

)+

1

d2 G≶ (12) Σ− (21 ) ,

t0

(4.2)

64

where the notation 1 ≡ (x1 , t1 ) has been introduced. The retarded (+) and advanced (−)
Green’s functions are defined according to

G± (1, 2) = Gδ (1, 2) ± Θ [±(t1 − t2 )] G> (1, 2) − G< (1, 2) ,

(4.3)

where Gδ stands for a possible singular contribution at t1 = t2 . The self-energies Σ±
are defined similarly, except that the singular contribution is instead explicitly included
in the first term on the r.h.s. of Eq. 4.2 as ΣHF , to emphasize its mean field nature and
to foreshadow the approximation in the next section. The generalized self-energy Σ(1, 2)
introduces interaction effects on the time evolution and also describes excitation processes
within the system [58, 60].
The self-energy in the previous equations has been separated into two different components. The first one involves the Hartree-Fock (HF) contribution, ΣHF (x1 , x2 ; t), which
accounts for the instantaneous, one-body interaction of the considered particle with the
mean-field produced by the other particles of the system. The term involving Σ≶ describes
time-dependent excitation processes, beyond the mean-field changes. Such contributions account for the effect of correlations on the dynamics and need to be included for a complete
description of nuclear reactions.
1D calculations are rather economical in terms of CPU and storage. Let us consider,
though, the difficulties associated to storing density matrices (or Green’s functions) in D
dimensions without explicit spin or isospin degrees of freedom. A uniform mesh of size Nx
in each direction will yield matrices of size Nx2D . Fairly sparse meshes with Nx ∼ 10 are
already computationally demanding for D > 2. For 1D this is not an issue, and we present
results with Nx ∼ 100, which can be carried out straightforwardly, so that we may study a

65

tractable situation and see how to reduce the computational cost before moving to 3D.
In addition to spatial variables, NGF requires a suitable infrastructure to keep track of
the time degree of freedom, as memory effects come into play. From Eq. 4.2, one can see that
the time derivative of the Green’s function G< at times t1 and t1 depends on the Green’s
functions and self-energies at all the previous times t, t0 < t < t1 and t0 < t < t1 , via
the time integrals on the r.h.s. Consequently, to find a solution of the KB equations, one
must keep track of all the previous time-steps. For a uniform mesh of Nt points in time, the
G≶ functions thus require handling Nx2D × Nt2 matrices. For D > 1, this becomes a major
concern in the numerical implementation of the KB equations [81].

4.2

Mean field approximation

The KB equations simplify substantially when the correlation effects, described in terms of
Σ≶ , are neglected. In that case, the evolution equations for G< and G> can be decoupled.
Since the single-particle observables are more straightforwardly expressed in terms of G< ,
we concentrate on the evolution of that function. Thanks to the instantaneous nature of
ΣHF , the set of equations for the time-diagonal elements of the Green’s function, t1 = t1 ,
can also be closed and there is no need to consider explicitly time off-diagonal terms. Note,
in particular, that this reduces the dimensionality of the problem to Nx2D , as one can update
G< at every timestep, and there is no need to keep explicit track of expectation values
involving the past.
The evolution equations simplify even more when assuming a negligible range for the

66

particle-particle interactions. In that case, the evolution equation becomes

i

2
∂ <
G (x1 , x1 ; t) = −
∂t
2m

∂2
∂2
−
∂x21 ∂x2

+ U (x1 , t) − U (x1 , t) G< (x1 , x1 ; t) , (4.4)

1

where U is related to ΣHF , in the local limit, with ΣHF (x1 , x2 ; t) = δ(x1 − x2 )U (x1 , t). This
can easily be linked to TDHF by expanding the Green’s functions in terms of single-particle
states [70]. For the same initial state and mean-field parametrization, both approaches yield
identical dynamics. Our goal here is not so much to provide an accurate description of 1D
nuclear mean-field dynamics, as this has been a path well tread by others in the past, but
rather to explore its implementation in terms of Green’s functions.

4.3

Collision of one-dimensional nuclear slabs

As a first example of implementation, we investigate the collision of two one-dimensional slabs
of spin-isospin saturated nuclear matter. Thus, each shell (energy state) contains 4 nucleons,
a combination of spin-up and spin-down neutrons and protons. Consider the evolution of
the one-body density matrix, found via the Green’s functions,

ρ(x, x , t) = −iG< (x, t; x , t)

(4.5)

In the 1D interpretation, the nucleon density is given by

n1 (x, t) = νρ(x, x; t),

(4.6)

where ν is the nucleon degeneracy, ν = 4, multiplied by the diagonal of the density matrix.

67

In this mean field approximation, we would like to interpret the results in a 3D situation,
even though it is a 1D simulation. The 3D density, n, can be related to the 1D density
by [70]
n(x, t) =

5
3

πn20
6ν 2

1/3

n1 (x, t),

(4.7)

where n0 is the nuclear saturation density, n0 = 0.16 fm−3 . This relation allows us to employ
well-known 3D mean field parameterizations for 1D calculations. We choose a simple nuclear
mean field parameterization for our calculations,

U (t) =

2+σ
3
t0 n(x, t) +
t3 [n(x, t)]1+σ ,
4
16

(4.8)

with the free parameters, t0 , t3 , and σ, fitted to the saturation properties of nuclear matter [70] (Here, σ has no relation to the cross section. It is merely an unfortunate homonym).
For such a local mean field, the time evolution of the density matrix can be implemented
numerically in a rather straightforward way using the Split Operator Method [70,82]. This involves a repeated switching between position and momentum space, optimally accomplished
using Fast Fourier Transforms (FFT) [82], implemented here with the FFTW library [83].

4.3.1

Preparation of the initial state

To describe the time evolution of a system, one needs a well-defined initial state. To describe a
collision between nuclei, the nuclei start in their ground state (besides the necessary center-ofmass motion). However, the potential described in Eq. 4.8 does not have an analytic solution
for the ground state, so some numerical method needs to be used to produce it. In the BUU
implementation used in Chapters 2 and 3, the Thomas-Fermi approximation is used. Instead,

68

we would like to use the same computational framework to initialize the ground state as we
do for the time evolution. This becomes important for describing correlated initial states,
and it reduces development time due to code reuse.
In particular, the initial state for a reaction (or any other dynamical process) should
be a ground state within the same approximation scheme as used in the dynamics. This
guarantees that observables for the ground state do not change with time, and thus it is
actually a ground state.
Since it is more reliable to start with an analytic ground state, we slowly switch from a
potential with such a state to our desired potential as we evolve forward in time, a method
known as adiabatic switching [84]. If switched smoothly and slowly enough, then the density
matrix will move to be, at every step, in the ground state of the current potential. Therefore,
at the end of the switching, the density matrix represents the ground state of the desired,
non-analytic potential.
For the collision of 1D nuclear slabs, we use a harmonic oscillator (HO) initial state
[70]. With the adiabatic switching, the single-particle potential acquires an explicit time
dependence and we adopt the switching potential

Uswitching = F (t)U0 + [1 − F (t)]U ,

(4.9)

where U0 is the initial HO potential, U0 = 12 mΩ2 x2 , U is our mean field potential, and F (t)
is the switching function. Ω is chosen to minimize the shape difference between the ground
state of the two potentials, to speed up the transition. For the switching time evolution from
time ti to tf , we choose F (t) such that it transitions monotonically and smoothly from 1 to
0 with a transition timescale τ , which should be longer than any characteristic time of the

69

system.
Slower adiabatic switchings should provide better approximations to the mean-field ground
states. In the top panel of Fig. 4.1, we show the evolution of the energy per nucleon for a
HO slab with initial time t = −1000 fm/c. We concentrate on a single case, with Ns = 2
HO shells filled (i.e. A = 8 nucleons in the 1D interpretation). We consider different transition times, τ , that characterise the adiabaticity of the transition. For any of the employed
transition times τ ≥ 5 fm/c, the energy evolves to a value very close to that for the static
Hartree-Fock solution (see inset). Judging the quality of the adiabatic transition on the
basis of the energy alone is potentially treacherous, though, as the energy is quadratic in the
deviation of wave function from the ground state. In principle, final states for the adiabatic
evolution might be found, with a wavefunction poorly approximating that of the mean-field
ground-state, but giving an energy close to the ground-state value.
Alternatively, the total density, n, is linear in the deviation of the wavefunction from the
ground state, so it may provide a better measure of the wavefunction quality. With this in
mind, in the bottom panel of Fig. 4.1, we show the evolution of the size of the slab, defined
as twice its average extent in 1D. Indeed, for τ = 5 and 10 fm/c, slab sizes exhibit significant
oscillations in the final state, indicating that the ground state has not been yet satisfactorily
reached. Instead, the final state in the dynamics is an excited state of the final hamiltonian
and thus fluctuates over time. For τ

30 fm/c, the oscillations become insignificant, with

the slab sizes practically coinciding with that of the static Hartree-Fock solution.
Now we have the initial state of one nuclear slab. To simulate a collision, they must be
made to collide with each other. In practice, this is acheived with a momentum translation
or “boost” operator [85]:
ρK (x, x ) = eiKx ρ0 (x, x )e−iKx
70

(4.10)

2<|x|> [fm]

E/A [MeV]

Adiabatic switching
20
10
0
-10
-20
-30

-28.5

τ=5 fm/c
τ=10 fm/c
τ=20 fm/c
τ=30 fm/c
τ=40 fm/c
Static HF

-29
-29.5

-550

-500

-450

3.2
3.1
3
-1000

-800

-600
-400
t [fm/c]

-200

0

Figure 4.1: Time evolution of the energy per particle (upper panel) and the size of the slab
(lower panel) when adiabatically switching from a starting HO configuration to a final mean
field state. Different values of the transition time τ are considered. The inset in the top
panel shows a magnified portion of the time evolution of the energy. For reference, mean
field results from static Hartree-Fock solutions are shown as straight solid (dark yellow) lines.

71

where ρK corresponds to a moving slab with momentum per nucleon K, while ρ0 is that
of a stationary slab. The complete initial state of the colliding system consists of two such
slabs moving towards each other.

4.3.2

Dynamics of colliding slabs

In order to explore the importance of the off-diagonal elements of the density matrix, we
take as an example the case of colliding nuclear slabs. Here, we examine the dynamics of this
collision to provide firm ground to investigate the off-diagonal elements in the next section.
The relatively simple 1D model of nuclear collisions discussed here demonstrates a surprisingly rich range of phenomena [70, 86, 87]. Qualitatively different physical processes are
observed within the model when changing the center-of-mass (CM) energy for the reactions. At low collision energies (ECM /nucleon ∼ 0.1 – 0.5 MeV), the nuclear slabs fuse
into one compound slab that remains excited for long times. For intermediate energies
(ECM /nucleon ∼ 0.5 – 15 MeV), a fusion process is observed, followed by a break-up into a
number of smaller slabs. Higher reaction energies (ECM /nucleon > 15 MeV) yield a pile-up
of density at the system center, followed by a violent break-up phase with the formation of a
low-density neck. The process is reminiscent of multifragmentation in nuclear reactions and
the results of the 3D transport model used in Chapter 3.
We do not explore this wide phenomenology here, but rather focus on one particular case
that will highlight the importance of off-diagonal matrix elements in ρ(x, x ), since we will
be truncating those elements later. Specifically, we concentrate on the collision energy range
immediately above fusion reactions. When the two colliding slabs fuse, the compound system
is too excited in energy to remain bound and the system undergoes a fission process into
three fragments. Because of the symmetry of this specific collision, the final state is formed
72

0.1

0

0
t=75 fm/c

0.2

t=100 fm/c

t=125 fm/c

0.1

0.1

0

0
t=150 fm/c

0.2

t=175 fm/c

t=200 fm/c

0.1
0

0.2

0.2
0.1

-20 -10 0 10 20
x [fm]

-20 -10 0 10 20
x [fm]

-20 -10 0 10 20
x [fm]

0

-3

Density [fm ]

0.1

0.2

-3

t=50 fm/c

Density [fm ]

t=25 fm/c

-3

t=0 fm/c

0.2

Density [fm ]

-3

Density [fm ]

-3

Density [fm ]

-3

Density [fm ]

ECM/A=4 MeV

Figure 4.2: Evolution of the center-of-mass density profile for a collision of two A = 8 slabs
at the collision energy of Ecm = 4 MeV/nucleon.
by a central fragment and two outgoing, identical slabs. Fig. 4.2 provides the snapshots
for the time evolution of one such collision at ECM /nucleon = 4 MeV. Following the fusion
of the original slabs at t ∼ 50 fm/c, two density peaks form at the edges of the compound
system and subsequently emerge from the central region to opposite sides. The central region
eventually recontracts and oscillates, as it is a highly excited state of an A ∼ 8 system. The
A ∼ 4 fragments remain excited too, as evidenced in the small changes with time of their
central density.

73

4.3.3

Off-diagonal structure of the density matrix

In order to determine whether we can neglect the far off-diagonal elements in the density
matrix, we first investigate how the structures there are created.
At equal arguments in a specific basis, the density matrix ρ of Eq. (4.5) yields the singleparticle density in that representation. For example, in the reciprocal basis, ρ(k, k = k, t)
yields the single-particle momentum occupation. At different arguments (e.g. k = k ), the
off-diagonal matrix elements reflect correlations between the density at different locations
in the domain. As we have discussed earlier, the task of following all the elements in 3D
is likely to overwhelm present computer storage capabilities. Nevertheless, the quantities of
most direct physical interest, including densities, tend to be associated with either diagonal
or near-diagonal elements of the matrix [60]. If the physical system is such that we can safely
neglect far off-diagonal elements, NGF calculations would become practical, due to a drastic
reduction in dimensionality. We argue that central nuclear reactions are good examples of
this, as the late stages involve limited phase coherence between different nuclear fragments.
To examine and quantify the importance of elements far from the diagonal in space
representation, let us explore the off-diagonal structure of ρ(x, x ) in the same reaction as
presented in Fig. 4.2. The upper (lower) panels of Fig. 4.3 show intensity plots for the real
(imaginary) parts of ρ(x, x , t) for three different times in the reaction, chosen to represent
the initial, overlap and late (“freeze-out”) stages. The diagonal of Re{ρ} (upper panels),
x = x , coincide with the density profiles displayed in Fig. 4.2. The complex character of
the density matrix at initial time reflects the fact that the slabs have been boosted to have
a momentum directed toward the center of the system.
To understand off-diagonal structures, it is useful to consider the expansion of the density

74

20

t=0 fm/c

t=75 fm/c

t=200 fm/c

0.25
0.15

x’ [fm]

10

0.05
0

-0.05

-10

-0.15

-20

-0.25

20

0.15
0.10

x’ [fm]

10

0.05
-3
0.00 [fm ]

0

-0.05

-10
-20
-20 -10

[fm-3]

-0.10
0 10
x [fm]

20

-10

0 10
x [fm]

-20 -10

0 10
x [fm]

20

-0.15

Figure 4.3: Intensity plots representing the real (upper panels) and the imaginary part (lower
panels) of the scaled density matrix, ρ(x, x , t), for a collision at Ecm = 4 MeV/nucleon. Offdiagonal correlations develop as the slabs interpenetrate and fragments continue on.

75

matrix in terms of wavefunctions,

nα φα (x, t) φ∗α (x , t) ,

ρ(x, x , t) =

(4.11)

α

where nα is the occupation fraction of the state α. The wavefunctions of the initial slabs
are confined to within the size of the slab and, consequently, ρ will have a limited, square
support in the (x, x ) plane. The left panels of Fig. 4.3 show this structure around the two
initial slabs.
At this collision energy, a compound slab is formed (central panel, 75 fm/c) in the centerof-mass. The two outgoing fragments and the central compound system are clearly identified
as blobs near the diagonal, x = x . Cross-correlations develop between individual fragments,
as signalled by the patches of significant values of the density matrix far away from the x = x
diagonal. As an example, the structure at (x ∼ 0, x ∼ 15) fm for t = 200 fm/c in Fig. 4.3,
both for the real and imaginary parts of the matrix, signal the overlap of single-particle states
between the central fragment and the fragment that is leaving in the positive direction. That
is, a state that was localized within one of the fragments is now spread between these two
fragments, carrying with it its correlated phase.
The discussed correlation patches may be understood in terms of the fragmentation of
natural orbitals. Orbitals can have components in each fragment rather than being split
among them. Consequently, if those two fragments are in positions x1 and x2 , when taking
the product in Eq. (4.11), one will find non-zero elements of ρ around (x1 , x2 ). Physically,
this is a reflection of the fact that, during the breakup process, the nucleons from the original
orbitals have finite probabilities of ending up in different fragments. The amplitudes for those
probabilities maintain a phase relationship leading to correlations in the density matrix. The

76

entanglement of the internal wavefunctions is thus the sole reason for the persistence of the
far-away off-diagonal patches in the density matrix. The real and imaginary parts for those
structures are comparable in magnitude, as seen in Fig. 4.3. This points to the involvement of
significant relative phases, as expected from the large difference in momentum and position.
Note that these entanglement correlations persist for fragments that are 30 fm apart for the
t = 200 fm/c panels of Fig. 4.3.
The entanglement between these far-away fragments, represented by the nonzero offdiagonal elements, will only affect the dynamics if the fragments interact again — a consequence of the irrelevance of global phase. In the context of a nuclear collision, this is unlikely,
as particles tend to experience a Hubble-like expansion, no longer interacting with particles
that have dissimilar momenta. This suggests that perhaps the off-diagonal elements can be
safely ignored, leading to truncating them and reducing the computational cost. We test
this in the next section.

4.4

Testing importance of off-diagonal elements

Physics beyond the mean-field dynamics (short-range correlations, particle and gamma decays, etc.) are likely to introduce additional decoherence between the separating fragments,
beyond that stemming only from different mean-field orbitals. Correspondingly, higher values for the off-diagonal elements are likely to persist more in a mean-field approach than in
any simulation with correlations. As a corollary, if one finds out that the far off-diagonal
elements may be disregarded within the mean-field dynamics, then it should be even more
possible to disregard such elements in more realistic approaches.
To quantify the importance of off-diagonal elements, we want to force them to be zero

77

at all times and see how it affects the dynamics. This suppression must be done smoothly,
however, to avoid the density matrix elements reflected from the suppression and creating
spurious oscillations. To acheive this in the evolution, we employ an ad hoc procedure based
on super-operators1 that eliminates off-diagonal elements in the time evolution [70]. Fig. 4.4
sketches the spatial structure of the super-operator. White areas, beyond x − x > x0 + d0 ,
are artificially eliminated at each time step. The remaining blue areas are left unaffected.
Note that because we are using a fourier transform to access the momentum space, the
structure is periodic in position space. Any operation we perform should preserve this
periodicity, thus the lack of suppression at the (−L, L) and (L, −L) corners. Fixing the
smoothness factor d0 = 2 fm, one can progressively eliminate more and more off-diagonal
matrix elements by decreasing x0 . For the break-up reaction explored here, the results of
this procedure are shown in Fig. 4.5. Neither the total energy nor its kinetic or potential
components (upper panel) are significantly affected by these cuts. The spatial extent of the
system, presented in the bottom panel, is also mostly unaffected, pointing towards a good
description of the density distribution. Overall, the elimination procedure we have devised
indicates that off-diagonal elements play a very small feedback effect on the dynamics.
The off-diagonal suppression procedure we have just presented, though, does not directly
involve any savings in terms of storage or CPU time, since though these off-diagonal elements
are now zero, we are still using them for computation. The most substantial savings can
be achieved by using an alternative spatial discretization scheme in a rotated coordinate
frame. With this, we can also control the spareness and lengths of the mesh in the relative
coordinate while keeping the average (physical) direction unaffected. This is investigated in
1A

super-operator is an operator that acts on other operators. In this case, it acts on the
density operator.

78

x

x

=

x

L

0

x0

d0

L − x0 − d0

−L

x
−L

0

L

Figure 4.4: Schematic of suppression of off-diagonal elements. White (non-shaded) areas are
heavily suppressed. The upper-left and lower-right corners are not suppressed, to preserve
periodicity in both coordinates.

79

ECM/A=4 MeV
E/A [MeV]

20
0
-20

Kinetic

Total
Potential

-40

2<|x|> [fm]

25
x0=25 fm
x0=20 fm
x0=15 fm
x0=10 fm

20
15
10
5
0

50

100
t [fm/c]

150

200

Figure 4.5: Time evolution of the energy per particle (upper panel) and of the spatial extent
(lower panel) in a collision between 1D nuclear slabs at ECM /nucleon= 4 MeV, for different
values of the cutoff parameter x0 at fixed smoothness d0 = 2 fm.

80

the next section.

81

4.5

Rotated density matrix

Now that we’ve shown that the far off-diagonal elements of the density matrix can safely be
suppressed, we describe a technique to practically truncate the matrix in the off-diagonal
direction by recasting it into a rotated coordinate frame. We also show how this rotation
partially decouples the position and momentum discretization, permitting access to arbitrary
regimes of kinetic energy without altering the resolution and range of the 1D box in position
space. First we provide justification for why it is helpful to rotate the coordinate system,
rather than simply rewrite the Fourier transforms.
Naively, in order to not compute the unnecessary off-diagonal regions of the density
matrix, we could rewrite the Fourier transform code to simply not consider them. If we used
an off-the-shelf FFT library, then each row would also have a different periodicity, since
each row has a different range. However, much work has gone into making FFT libraries
very efficient, especially for rectangular multi-dimensional transforms. Therefore, our first
step is to set up a rectangular density matrix, with one of the directions being the direction
that we truncate. The next step, not explored here, would be to use parallel processing FFT
implementations. Adaptive Mesh Refinement (AMR) should also be explored, so as to only
compute important regions, irrespective of where they lie in the density matrix.
To be able to truncate in this way, we need one direction of the grid to be in the offdiagonal direction of the (x, x ) matrix. As one travels in this direction away from the x = x
diagonal, we notice that x − x increases, while (x + x ) remains constant. This suggests
that a convenient coordinate system would be the average and relative coordinates, which we
call xa and xr , respectively. Note that this is not a transform into a center-of-mass system,
since x and x refer to different positions of the same object. In this new representation, the

82

line that gives the probability density is now xr = 0, compared to x = x in the previous
coordinate system.

4.5.1

Fourier transform in the new coordinates

In the unrotated, (x, x ) system, the transform to momentum space takes the form

ρ(k) (k, k ) =

1
2π

dx dx ρ(x, x )e−ikx eik x

(4.12)

The second exponential factor has no negative sign in the argument, tied to the fact that the
density matrix is hermitian-symmetric under exchange of position or momentum arguments,
ρ(x, x ) = ρ† (x , x). The 1/(2π) prefactor is in both the transform and inverse transform, in
order to produce the suitably normalized probability density in momentum space.
In the rotated, (xa , xr ) system, we make the substitutions

qa =

q+q
,
2

qr = q − q ,

(4.13)

where q represents x or k. Since these are not actually center-of-mass and relative coordinates
of a 2-particle system, we do not need to ensure that ka is the conjugate of xa . Therefore, we
are free to choose 1/2 as a factor in both of those coordinates. This allows a straightforward
interpretation of xa as the “physical coordinate” — that is, if qr = 0, then qa gives the
physical position or momentum, where q stands in for x or k in turn.
Substituting these coordinates in Eq. 4.12, we arrive at the transform

ρ(k) (kr , ka ) =

1
2π

dxa dxr ρ(xa , xr )e−ika xr e−ikr xa .

83

(4.14)

This substitution gives a new transform that is identical in structure to the previous one,
making it amenable to direct use of FFT libraries. Note that xr (xa ) transforms into ka
(kr ).

4.5.2

Discrete Fourier transform

Since we are assuming periodic boundary conditions in position space, this creates a discrete
spectrum of allowed momenta, the resolution ∆k being inversely proportional to the length
of the periodic box in position space, xmax . Further, as we are approximating a continuous
space with a discrete position basis with a finite number of positions, there is a maximum
momentum, kmax , that can be described by the calculation, inversely proportional to the
resolution in position space, ∆x (the closer together the positions are, the more-quicklyvarying density can be described, which corresponds to a higher momentum). These relations
are true for both x and x . In particular, for a box that has a length 2xmax , the relations are

∆k =

π
xmax

,

kmax =

π
,
∆x

(4.15)

and similar for ∆k and kmax . Note that this is for the choice of position range −xmax to
xmax and momentum range −kmax to kmax .
In the rotated system, ka (kr ) is conjugate to xr (xa ) in the Fourier transform, and the
relations become
∆ka(r) =

π
(max)

(max)

,

ka(r)

xr(a)

=

π
.
∆xr(a)

(4.16)

Thus, we can independently vary the resolution and range of xa and ka , by modifying the
characteristics of the relative coordinates. In the unrotated system, for a given size in position

84

space, for increasing the maximum momentum, the computational storage scales as (DN )2 ,
where D is the number of physical dimensions and N is the number of grid points in each
x and x direction (for an N × N density matrix). In this rotated system, one can vary the
number of points in the relative direction independently, making the computational storage
scale as (DN ) for an increase in maximum momentum described by the system.

4.5.3

Testing the rotated matrix

To validate the method, we first obtain exactly the same results as in the unrotated case.
In order to perform exactly the same operations, we must match the periodic of the original
system. This periodicity, illustrated in Fig. 4.6, is such that ρ(x, x ) = ρ(x + 2xmax , x ) =
ρ(x, x + 2xmax ). In this schematic example, we give a physical system that consists of a
large and small slab, represented as filled circles centered on the diagonal. The circles could
represent quantum harmonic oscillator states. There is a correlation between them (a nonzero phase factor), as evidenced by the fainter ellipses in the far off-diagonal region. We show
this assymetric, correlated system to provide orientation throughout the following rotation
and subsequence truncation in the relative direction.
When rotated into (xa , xr ) coordinates, we select an area that is a) rectangular, to be
convenient for the FFT, and b) is periodic in both the average and relative directions.
Fig. 4.7 shows our choice.
However, this area is twice the original area, hence twice the computational storage of
the original matrix (and over twice the computation time for the FFT). This is because
the orginal (x, x ) matrix appears twice in the rotated one. The maximum extent in xr is
(max)

xr

(max)

= 2xa

. In Fig. 4.8, regions with the same pattern are representing the same region

of the unrotated matrix. Using this duplication, we rewrite the Discrete Fourier Transform
85

x

x

Figure 4.6: Schematic of density matrix in (x, x ) coordinates, with example of periodic
boundaries and the resulting physical system that is simulated. Each square contains an
identical copy of the physical system.

86

x
x

xr
xa

Figure 4.7: With the density matrix rotated, the smallest rectangle that is periodic in
both the average coordinate (left-to-right) and relative coordinate (down-to-up) is marked
in dashed lines.

87

(DFT) to span only the innermost 4 rectangles, reducing the maximum extent in the relative
(max)

direction to xr

(max)

= 2xa

, producing the rectangle in Fig. 4.8 that is bounded in the xr

direction by dot-dashed lines.
However, this configuration leads to a surprising result — in momentum space, ρk (kr , ka ),
every grid point where the sum of the indices, ka + kr , is odd, the density is identically zero,
creating a checkerboard pattern. In the derivation, this results from a factor in the DFT,
[1 + (−1)ka +kr ], that appears when the duplicated density matrix from Fig. 4.8 is rewritten
as a non-duplicated sum. Conceptually, this can be seen to arise from the rotation itself, as
seen in Fig. 4.9.
We remove this duplication by separating the Fourier transform into two — one for even
kr and ka , and one for odd. These transforms do not communicate between each other, and
so the N × N Fourier transform is reduced to two N × N transforms, regaining the original
order of computational cost.
With this work done, we now have a rotated density matrix that can be computed with
the same computational cost as the original. Now we can perform the work we set out to do
at the beginning of this chapter — reduce the computational cost by truncating the matrix
in the xr direction.
Fig. 4.10 illustrates the implementation. In the left-most panel, the suppression of faroff-diagonal elements is displayed. Again, a smooth transition to total suppression is used,
seen by the shading (blue) gradient in the figure. The resulting suppression of the original
(x, x ) density matrix in Fig. 4.10b (due to the duplicated regions described earlier) can be
compared to Fig. 4.4, showing that the suppression is identical, leaving the corners of large
x − x = |xr | unsuppressed. To perform the actual truncation, we simply reduce the extent
of the matrix in the xr direction, producing the matrix bordered with a dot-dashed line in
88

x
x

xr
xa

Figure 4.8: Schematic illustrating duplication of original density matrix. Regions with the
same pattern have the same exact part of the unrotated (x, x ) density matrix in them.

89

(a) unrotated

(b) rotated

(c) extra zeroes

Figure 4.9: Rotation of matrix and creation of fictitious points. The x = x line is shown
as a dotted line, for orientation (xr = 0 line in rotated system). When the matrix in (a) is
rotated into (b) and the system is modeled with a rectangular Fourier transform, extra grid
points are created, shown as open circles in (c).
the right-most panel. The duplicated regions, bordered with dashed lines, are still modeled,
since they are duplications of the computed matrix.
Here, there is a complication regarding the use of stock FFT libraries. The standard
DFT is

N/2−1

fj e−i2πjk/N

f˜k =

(4.17)

j=−N/2

If we reduce the extent in xr without changing the number of points, then the DFT remains
the same and we get the exact result compared to the suppressed, untruncated case. Since
(max)

this increases the resolution in xr , the maximum momentum ka

is increased and we can

access higher energy regimes without increasing the computational cost. However, if we want
to change the number of points while keeping the DFT the same, we run into a problem.
If we reduce the number of points by a factor ν and only change the range of indices in the

90

x

xr
xa
x

xr

xr
xa

(a) suppressed regions

xa
(b) periodic view

(c) matrix truncation

Figure 4.10: Implementation of off-diagonal (xr -direction) suppression. Suppressed areas
are blue (shaded). Implementation is in (a), while the resulting suppression of the original
density matrix shown in (b) (compare to Fig. 4.4). Truncation of the matrix is shown in (c),
with the computed matrix bordered with a dot-dashed line and periodic extensions shown
with dashed borders. The empty space is not computed at all, reducing the computation
cost. This benefit is amplified geometrically with an increase in the number of dimensions.
transform, while keeping ∆x the same, we get
νN/2−1

fj e−i2πjk/N

f˜k =

(4.18)

j=−νN/2

However, if we feed these νN points to a stock DFT algorithm, it will compute
νN/2−1

fj e−i2πjk/(νN ) ,

f˜k =

(4.19)

j=−νN/2

thus raising each phase factor to the power 1/ν. Until we learn how to rewrite this in the
standard form, we resort to using our own Fourier transform algorithms. Since Fortran 95
introduced built-in, optimized matrix multiplication routines, there is not a large sacrifice in
computation time.
We have shown how to practically truncate the density matrix to reduce the computational cost. This method also results in the ability to tune independently the resolution and

91

range of position and momentum. In the next section, we show a practical application of
this 1D density matrix time evolution in modeling the Bose-Einstein condensate (BEC).

4.6

Bose-Einstein Condensates

Before moving to a 3D simulation or adding correlations (both are beyond the scope of this
thesis), we demonstrate the use of this code in a physical 1D quantum system, in particular
the Bose-Einstein condensate (BEC).
A BEC is a collection of bosons, the majority of which are in the same ground state. In
magneto-optical trap (MOT) experiments, there are typically

104 atoms in this ground

state, their wavefunction spread out across ∼ 0.1 mm. An ideal BEC has all of the particles
in the ground state. This occurs below a certain critical temperature associated with Bose
statistics.
In this section, we first describe the mean field model for BECs. Then, we use the density
matrix formalism developed in this chapter to find the ground state of that model to find
the density distribution of the BEC ground state.

4.6.1

Modeling with the Gross-Pitaevskii equation

The mean field description of the BEC is given by the Gross-Pitaevskii equation (GPE) [88].
This models a dilute gas of atoms interacting with a contact potential, with particles moving
slowly enough to only scatter with binary s-wave interactions. The resulting equation is
the Schr¨odinger equation with a potential term that is linearly dependent on the density.
And since the condensate is all in one quantum state, a single wave function can be used to

92

describe all the atoms,
2

−

2m

∇2 + Vext (r) + g |φ(r)|2 φ(r) = µφ(r) ,

(4.20)

where Vext (x) is an external potential, typically a harmonic oscillator trapping potential,
and g is the interaction strength. g is dependent on the s-wave scattering length as and the
mass m of the particle, g = 4π 2 as /m.
For a 1D condensate, an elongated trap is used, where a steep potential is created in 2
perpendicular directions, called transverse or radial (HO with frequency 2πω⊥ ), and a much
shallower potential is used in the third, longitudinal direction (frequency 2πωz ). Due to
this anisotropy, excitations in the transverse direction are limited to zero-point fluctuations,
and thus the dynamics in the transverse direction are frozen.2 When these directions are
integrated over, we arrive at the 1D GPE,
2

−

∂2
+ Vext (x) + g1D |φ(x)|2 φ(x) = µφ(x) ,
2m ∂ 2 x

(4.21)

where g1D = 2as ω⊥ [89]. Thus, the interaction potential is dependent on the number of
particles in the trap, unlike in the 3D case.
In general, this equation cannot be solved analytically. In the Thomas-Fermi approximation (the limit where potential energy is much greater than the kinetic energy) and with a
harmonic oscillator external potential Vext (x) = (1/2)mΩ2 x2 , the kinetic term is neglected,
2 Strictly

speaking, this is a so-called quasi-1D condensate, as the critical temperature for
condensation in 1D tends to zero (see Appendix C).

93

and the probability density of the ground state is an inverted parabola,




1


µ − 12 mω 2 x2

g
|ψ(x)|2 =





 0

,

,

|x| <

2µ
mω 2

,

(4.22)

otherwise

the piecewise continuous function being necessary to preserve the positive-definiteness of the
probability density. Note that this sudden change in the wave function in position will cause
an extended population of momentum states. The practical effects of this are addressed in
a following section.
For the case where the Thomas-Fermi equation does not appply, a numerical method
to arrive at the ground state must be used. We use the adiabatic switching technique
described in Section 4.3.1. For this technique, we start with a wave function and potential
that is analytically solvable. Since most BECs are produced and contained in a harmonic
oscillator potential produced by a magneto-optical trap (MOT), we use this as our external
potential. An intuitive choice for the initial, analytic potential would be the bare harmonic
oscillator potential, with no density dependent interaction term. To minimize the number
of switching steps required, we select by eye the initial HO frequency that results in a
similar physical extent of the system as the final, density-dependent potential, as measured
by |x| . Another method would be to variationally solve for the frequency of the harmonic
oscillator potential whose ground state (a Gaussian) has the lowest energy in the final,
density-dependent potential. The resulting Gaussian width, σ, can be solved iteratively with

σ= 4

1
mωz

2
2
g1D σ +
π
2m

94

(4.23)

and the HO frequency is given by ΩHO = /(σ 2 m).
To do this, we need to choose a physical system, preferably one that has been experimentally realized. We choose for our system a gas of ultracold 87 Rb atoms, placed in an
elongated trap so as to be nearly one-dimensional and thus describable by our 1D code [90].
87 Rb

has an s-wave scattering length of as = 5.29 nm. Since it is positive, the interaction is

repulsive, and thus we need a trapping potential to contain them and actually form the condensate. The trap frequencies used in the experiment were 840 Hz in the transverse (radial)
direction and 14 Hz in the longitudinal (axial) direction. With a trap potential anisotropy
ratio of ωr2 /ωa2 = 3600, this is a clear candidate for the use of the 1D simulation — there
will be a large spectrum of excitations in the axial direction before any radial excitations are
produced.

4.6.2

Initial state

To produce the ground state of the BEC, we considered several different initial, analytic
potentials for which we could calculate the ground state. These included the infinite square
well, finite square well, harmonic oscillator, and Kronig-Penney model. The finite square
well and harmonic oscillator (HO) have the disadvantage that they are not periodic, while
our system is assumed to be periodic through our use of the Fourier transform. The infinite
square well, while not strictly violating periodicity, does not lend itself easily to adiabatically
switching to a different, finite potential. The Kronig-Penney model is perhaps the potential
that gives a ground state best suited to the periodic nature of our system. This potential
is a periodic system of rectangular potential barriers. Thus, the ground state is that of the
infinite periodic system that the Fourier transform assumes, and thus the ground state of
the Kronig-Penney model should be closer to our true ground state than the ground state of
95

other, nonperiodic potentials.
In the end, however, the Kronig-Penney ground state is complicated and not as intuitive
as that of the HO, so we use the HO potential, in order to manipulate it by hand more
easily. Its nonperiodic nature is not so relevant, since its ground state density is effectively
zero near the edges of the box.
The system of 87 Rb atoms described in the previous section was realized experimentally
with 5 × 105 atoms. To discover what the resolution and range for position and momentum
should be to maintain a stable evolution, we start simulating the system with just 1 particle,
then increase the number of particles in subsequent trials, since the interaction potential
is dependent on the number of particles N , as seen in Eq. 4.21. To maintain a similar
shape of the resulting condensate across different numbers of particles, the trap frequency is
proportionally reduced, keeping constant the ratio ωz2 /N and thus the shape of the potential.
To verify that the system smoothly transitioned to the new ground state, we examine
how the energy and spatial extent vary over time, shown in Figs. 4.11 and 4.12. For a very
small number of particles, N = 256, the resulting system is quite stationary, the energy
varying by about 0.6% over 1 s and the spatial extent varying by about 2% over 1 s at the
end of the adiabatic switching. Given that the relevant timescale of dynamics is ∼ 1 ms, this
is a very good approximation of the ground state. For greater numbers of particles, with
a matching increase in the trapping potential, the resulting ground state is less stationary,
with the spatial extent varying by about 3% over 0.5 s. As N increases and the trapping
potential frequency increases, we must adjust the position and time discretization to obey
stabilty constraints (see Appendix D).
We satisfy the momentum-space stability constraint in Eq. D.4 for the BEC ground state
found in the previous section. For that system, we also satisfy the position-space constraint,
96

18

N
N
N
N

16
Etotal /N [neV]

14

= 256
= 512
= 1024
= 2048

12

17.7

10

17.6

8

17.5

6

17.4

4

17.3
63.2

2

63.6

64

0
0

10

20

30

40
tadiabatic [s]

50

60

70

|x| [µm]

Figure 4.11: Total energy of system while adiabatically switching from an external HO
potential to the Gross-Pitaevskii potential. A reasonably stable energy is reached as the
number of atoms in the condensate, N , increases. However, as seen in the inset, as N
increases, the total energy exhibits an increasing oscillation, suggesting higher excitations.

65
60
55
50
45
40
35
30
25
20
15
10

N
256
512
1024
2048

60.5
60
59.5
59
58.5
58
63.2

0

10

20

30

40
tadiabatic [s]

50

63.6
60

64
70

Figure 4.12: Extent of system in position space while adiabatically switching from an external
HO potential to the Gross-Pitaevskii potential. The ground state is more closely reached for
a fewer number of particles in the condensate, as seen by a slower, lower-amplitude oscillation
at late times than for larger numbers of particles.

97

Eq.D.5, for N

50. As N increases, we increase the trap frequency Ω to maintain the same

shape of the potential. If Ω and N were to be increased to match the experiment described,
with Ω = 2π × 14 Hz and N = 5 × 105 , the position-space constraint would be exceeded
by 7 orders of magnitude, and satisfying this constraint would require a substantial change
in xmax and ∆t, potentially to degrees that unacceptably increase the computational cost.
This is migitated, somewhat, by the fact that the BEC never reaches the edge of the grid,
and thus it never experiences the maximum potential. If it only approaches, for example,
half the distance to the edge, then this would relax the position-space constraint by a factor
of 4. Still, more work needs to be done to allow the increase in particle number to the typical
size of experimentally realized condensates.
Fig. 4.13 shows the probability density ρ(x) ≡ ρ(x, x = x) before and after the adiabatic switching procedure. During the switching, the external HO transitions to a higher
frequency, and the density dependent potential turns on, pushing the density to spread out.
In this stationary system, with large N , the potential energy is much higher than the kinetic,
resulting in the applicability of the Thomas-Fermi approximation, with the solution for the
density given in Eq. 4.22 — an inverted parabola. In fact, this is the shape that is seen in
the figure, except for the edges of the parabola. Analytically, the Thomas-Fermi solution is
piecewise continuous here, with an undefined first derivative. Such a sharp change in the
density in position space would be represented as an infinite sum of waves in momentum
space, akin to the Fourier series of a square wave. Because of our finite grid, and a limited
maximum momentum, sharp boundaries cannot be represented. Instead, the edge becomes
softened, as seen in the inset in Fig. 4.13.
The 2D density matrix ρ(x, x ) is shown in Fig. 4.14, both before and after the adiabatic
switching. During the switching, the amplitude spreads in both the average direction, (x +
98

0.25

before adiabatic switching: HO
after: GPE

ρ(x) mm−1

0.2

0.016
0.012

0.15

0.008
0.004

0.1

0
0.12

0.05
0
-0.2

-0.15

-0.1

-0.05

0
0.05
x [mm]

0.16

0.1

0.15

0.2

0.2

Figure 4.13: Probability density at beginning (solid line, red) and end (dashed line, green) of
adiabatically switching from the harmonic oscillator potential to the Gross-Pitaevskii potential. Since the potential energy is much stronger than the kinetic, the resulting distribution is
the Thomas-Fermi approximation, an upside-down parabola, as in Eq. 4.22. The transition
between the parabola and zero density, seen inset, is smoother than the piecewise continuous
solution and more realistic.
x )/2, and relative direction, x − x . The circular support for the HO state has shifted
to a square-shaped support for the GPE state. The square shape means that the more
distant regions of the density (the upper right and lower left regions of density matrix) have
a stronger correlation with each other than in the HO case. This makes sense for a BEC
since the whole system is one coherent macroscopic state.

99

0.2
0.15

0

0.1

-0.2

0.05

0.4

0.05

0.2

0.04
0.03

0

0.02

-0.2

0.01

-0.4
0
-0.4 -0.2 0 0.2 0.4
x [mm]

-0.4
0
-0.4 -0.2 0 0.2 0.4
x [mm]

(a) HO

(b) GPE

ρ(x, x ) mm−1

0.2

x [mm]

0.25
ρ(x, x ) mm−1

x [mm]

0.4

Figure 4.14: Density matrix for (a) the analytic harmonic oscillator ground state before adiabatic switching and (b) after adiabatic switching to the Gross-Piteavskii equation, yielding
the Bose-Einstein Condensate. The local density lies along the diagonal x = x , shown in
Fig. 4.13. The square support in the BEC shows a stronger correlation between distant
regions of the condensate, exhibiting stronger long-range coherence than in the HO case.

4.7

Conclusion

With the goal of eventually describing 3D nuclear collisions including correlation effects
beyond the mean field with the Kadanoff-Baym equations, we have developed the apparatus
to simulate 1D systems in the mean field limit, showing that we can successfully reduce the
computational cost of the calculation by truncating the density matrix in the off-diagonal
direction. We used this model to describe 1D Bose-Einstein Condensates, thus using the
model to describe physical phenomena across 15 orders of magnitude in energy.

100

Chapter 5
Conclusion
In this thesis, we have simulated the time evolution of quantum many-body systems and
used comparisons to experimental data in order to learn more about the properties of nuclear
matter, understand better the dynamical processes in central nuclear collisions, and advance
the development of a nonequilibrium Green’s function description of both central nuclear
collisions and Bose-Einstein Condensates.
In Chapter 2, we determined the viscosity of nuclear matter by adjusting the in-medium
(med)

nucleon-nucleon cross section (σNN ) in a BUU transport model, using theory-motivated
reductions, until the simulation results matched experimental data on nuclear stopping.
Then we used that cross section to calculate the viscosity self-consistently. We confirmed
(med)

that it is necessary to reduce the σNN

to agree with experiment. Using our new viscosity

determination, we showed that the ratio of shear viscosity to entropy density, η/s, of nuclear
matter in intermediate energy collisions follows a trend that is qualitatively consistent with
results from collisions at relativistic energies.
In Chapter 3, we used the same BUU transport model as in the previous chapter to isolate
the protons emitted early in a central nuclear collision at intermediate energy, as predicted
in the model, using a filter on high transverse momentum. We predicted that there is a
detectable effect that this filter has on the experimentally observable source function. Then
we explored the prediction of a recontraction of protons at late times in the central collision
of 112 Sn +112 Sn at 50 MeV/nucleon that results in a secondary emission of protons. We
101

showed that the existence of this reemission would be detected by how the source function
changed with the filter on high transverse momentum.
In Chapter 4, we developed an early implementation of a more fully quantal transport
model, with our sights set on solving central nuclear collisions in 3D using nonequilibrium
Green’s functions. In our 1D, mean field approximation, density matrix model, we demonstrated the initial state preparation and collision of 1D nuclear “slabs”, showing the successful
ground state calculation and qualitative features of low-energy nuclear collisions. With the
aim of reducing the computational cost of the calculation, in anticipation of more dimensions and costly correlations, we showed that we can neglect far off-diagonal elements in the
density matrix, exploiting the proximity of the nuclear system to the classical limit, without
afffecting the one-body observables.
Further, we described a method of recasting the density matrix in a rotated coordinate
system, enabling us to not only ignore the irrelevant matrix elements in the time evolution,
but also avoid computing them completely, reducing the computational cost. As an added
benefit, we found that the rotation allows us to partially decouple the position and momentum discretization, permitting access to arbitrary regimes of kinetic energy without altering
the resolution and range of the 1D box in position space. Finally, we exhibited the wide
applicability of this density matrix approach by applying it to a system of 2000 ultracold
87 Rb

atoms in a Bose-Einstein condensate, as described by the Gross-Pitaevskii equation,

successfully acheiving a stable state in a harmonic oscillator trap.

102

APPENDICES

103

Appendix A
Equations for Source Imaging
Chapter
Momentum of a projectile nucleon in the reaction center of mass, p1,cm , given the kinetic
energy in the laboratory frame Tlab , mass of a nucleon m0 , and mass number of projectile
and target A1 and A2 is given by
p21,cm = −

1 A21 + A22 2
m0
2 A21
p21,cm

−

+

1
2

+ m20

p21,cm

+

2
A2
m
A1 0

(A.1)

2
A1 + A2
A
m0 + 2 Tlab m0 .
A1
A1

This equation is solved iteratively. The momentum of a target nucleon is correspondingly

p2,cm =

A1
p
.
A2 1,cm

104

(A.2)

Appendix B
Parameterizations of the source
function

B.1

Half-width at half-maximum

The most leading-order approximation for describing the source function is r1/2 , the halfwidth at half-maximum. As the source function is generally seen to be Gaussian at short
ranges, this value could be used to find the width of the Gaussian (if there are no other
components super-imposed at short range). The other benefit is that it is agnostic to the
actually details of the source function, and provides a coarse, understandable description.
The primary question is then how to determine the maximum of the source function.
Fitting the source with a single Gaussian function is problematic, as the source falls off
exponentially at large distances (see Fig. 3.4 inset). The HiRA group has fit the value of
the source between 0.5 and 2.5 fm to a parabola symmetric around r = 0,

S(r) = a + br2 ,

(B.1)

using a as the maximum [7]. They did not use the value of S(r) near zero, as that value had
a large uncertainty.

105

B.2

Hump fit

To fit the entire source function, one can also assume that it the combination of two Gaussian
sources, one short- and one long-range. An example of this is the so-called hump fit, used in
RHIC physics [55] (here abbreviated for one dimension):

S(r) = exp −Fs

r 2
× exp −Fl
2rs

r 2
,
2rl

(B.2)

where Fs = 1/(1 + (r/r0 )2 ) and Fl = 1 − Fs . At short range, the first exponential dominates
4
as a Gaussian and the second goes as e−r . At a longer range, the first exponential reduces

to a constant while the second acts as a Gaussian.

B.3

Chung fit

However, looking at Fig. 3.3, the long-range behavior of S(r) is exponential. Naively, one
would use a function that is the superposition of a Guassian and an exponential. However,
since this is a one-dimensional function projected from three dimensions, the function should
be smooth across r = 0, and an exponential of the form e−r has a nonzero first derivative
at the origin. Thus, one could use a functional form like [56],
λgau
λexp
r2
− r2 + β 2
S(r) = √
exp
−
+
exp
,
3
2
N (Rgau , Rexp )
Rexp
2Rgau
2πRgau
N (Rgau , Rexp ) = 4πβ 3

K0 (z) 2K1 (z)
+
z
z2

2
Rgau
β
,β=
,z=
Rexp
Rexp

(B.3)

fgau is the fraction of proton pairs that contribute to the Gaussian part of the source
within this source function parametrization, found by expanding the second term of Eq. B.3
106

in small r and integrating over all space.
∞
4π
dr r2 Sshort−range (r)
λgau + λexp 0
∞
λgau
λexp
β + r2 /2β
4π
r2
exp
−
=
dr √
exp
−
+
3
2
λgau + λexp 0
N (Rgau , Rexp )
Rexp
2Rgau
2πRgau

fgau =

=

=

4π
λgau + λexp
1
λgau + λexp

2
Rgau
λexp
+
exp
−
√
3
2
N (Rgau , Rexp )
Rexp
2πRgau
√
3
2
Rgau
2πRgau
λgau + λexp
exp − 2
N (Rgau , Rexp )
Rexp

∞

λgau

0

dr r2 exp −

r2
2
2Rgau

(B.4)

107

Appendix C
Existence of the BEC in one
dimension
The occupation probability, f ( ), for a system of bosons in thermal equilibrium is

f( ) =

where β is the inverse temperature,

e−β( −µ)
1 − e−β( −µ)

,

(C.1)

is the energy of the state in question, and µ is the

chemical potential. Thus the occupation probability approaches infinity as µ approaches
from below — that is, as it costs less and less energy to increase the occupation of this state.
This divergence leads to the phenomenon of BECs [91].
Therefore, if there exists a critical density that makes the chemical potential go to zero,
then we have a condensate. This density ρ, in a volume of dimension D (in the non-relativistic
case), is
2

ρ(µ →

0− , T )

= (2s + 1)

dD p
e−p /(2mT )
(2π )D 1 − e−p2 /(2mT )
∞

= (2s + 1)

= (2s + 1)

2
dD p
e−np /(2mT )
D
(2π )
n=1
√
D
2
mT
dD x e−x /2
2π

108

(C.2)
∞

n−D/2
n=1

For D = 1 and D = 2, the sum over n diverges — hence, strictly speaking, an ideal BEC
cannot form in 1 or 2 dimensions.

109

Appendix D
Numerical stability in the density
matrix time evolution
Since we are describing continuous space with a finite set of discrete points, we must choose
suitable parameters for the resolution and range of position and momentum such that the
numerical approximations do not affect the physics predictions. We also use a finite timestep,
which leads to its own constraints.
We approximate continuous time with a series of discrete timesteps, and assume that the
Hamiltonian is constant over that timestep. Thus the time evolution operator for the wave
function is approximated,
t+∆t

i
Uˆ (t + ∆t, t) = exp −
t

ˆ )
dτ H(τ
(D.1)

i ˆ
≈ exp − H(t)∆t

Since there is a density-dependent mean field potential, and the density changes with
time, the Hamiltonian changes with time. So, to maintain the validity of the approximation
in Eq. D.1, the density should change slowly with time, and thus, the time evolution operator
that the density matrix is multiplied by should be not too far away from unity. In addition,
the Split Operator Method (SOM) that we are using for the time evolution (described in
Section 4.3) is accurate to O (∆t)3 , so ∆t must itself be kept small, as we shouldn’t
110

evolve too much with the part of the Hamiltonian that is diagonal in the position basis
before evolving with the momentum-diagonal part. Since we are employing the SOM, we
can address the momentum- and position-dependent terms of the time evolution operator
separately to determine the constraints on the resolution, range, and timestep parameters.
The factor that must be kept close to unity is the one multiplied by the density matrix.
The time dependence of the density matrix can be related to the its definition in terms of
wave functions (Eq. 4.11). Combining this with the wave function evolution in Eq. D.1, we
arrive at
ρˆ(t + ∆t) = Uˆ ρˆ(t) Uˆ †
(D.2)
nα |φα φα | Uˆ † .

= Uˆ
α

In momentum space, using just the kinetic term of the Hamiltonian, this becomes

ρ(k) (k, k

k2
∆t
, t + ∆t) = exp −i
2m

nα φα (k, t) φ∗α (k
α

2

k
∆t .
, t) exp i
2m

(D.3)

Thus, the most significant phase factor that a grid point in momentum space will be multiplied by is exp[i

2
kmax
∆t]. This maximum value for k is given by Eq. 4.15 (or Eq. 4.16 for
2m

the rotated system), which is in terms of ∆x. The argument of the exponential should be
small compared to π/2, since that would make the phase factor equal to i. This yields the
constraint
π 2 ∆t
2m (∆x)2

111

π
.
2

(D.4)

A similar result is attained for the position-basis evolution. In general, the constraint is

Vmax − Vmin

π
.
2

∆t

(D.5)

For the case when the largest part of the potential is an external harmonic oscillator, which
is the case in the trapped BEC system we study in Section 4.6, and when it is much larger
than the minimum potential at any position, the constraint becomes
mΩ2 2
xmax ∆t
2

112

π
.
2

(D.6)

BIBLIOGRAPHY

113

BIBLIOGRAPHY
[1] Norbert Herrmann, Johannes P. Wessels, and Thomas Wienold, “Collective flow in
heavy-ion collisions,” Annual Review of Nuclear and Particle Science 49, 581–632 (1999),
http://www.annualreviews.org/doi/abs/10.1146/annurev.nucl.49.1.581
[2] E. Colin, Rulin Sun, N. N. Ajitanand, John M. Alexander, M. A. Barton, P. A. DeYoung,
A. Elmaani, C. J. Gelderloos, E. E. Gualtieri, D. Guinet, S. Hannuschke, J. A. Jasma,
L. Kowalski, Roy A. Lacey, J. Lauret, E. Norbeck, R. Pak, G. F. Peaslee, M. Stern,
N. T. B. Stone, S. D. Sundbeck, A. M. Vander Molen, G. D. Westfall, and J. Yee,
“Splintering central collisions: Systematics of momentum and energy deposition for
(17–115) A MeV 40 Ar,” Physical Review C 57, R1032–R1036 (Mar. 1998), http://
link.aps.org/doi/10.1103/PhysRevC.57.R1032
[3] W. Reisdorf, A. Andronic, R. Averbeck, M.L. Benabderrahmane, O.N. Hartmann,
N. Herrmann, K.D. Hildenbrand, T.I. Kang, Y.J. Kim, M. Kiˇs, P. Koczo´
n, T. Kress,
Y. Leifels, M. Merschmeyer, K. Piasecki, A. Sch¨
uttauf, M. Stockmeier, V. Barret, Z. Basˇ
rak, N. Bastid, R. Caplar,
P. Crochet, P. Dupieux, M. Dˇzelalija, Z. Fodor, P. Gasik,
Y. Grishkin, B. Hong, J. Kecskemeti, M. Kirejczyk, M. Korolija, R. Kotte, A. Lebedev,
X. Lopez, T. Matulewicz, W. Neubert, M. Petrovici, F. Rami, M.S. Ryu, Z. Seres,
B. Sikora, K.S. Sim, V. Simion, K. Siwek-Wilczy´
nska, V. Smolyankin, G. Stoicea,
Z. Tymi´
nski, K. Wi´sniewski, D. Wohlfarth, Z.G. Xiao, H.S. Xu, I. Yushmanov, and
A. Zhilin, “Systematics of central heavy ion collisions in the 1 A GeV regime,” Nuclear
Physics A 848, 366–427 (Dec. 2010), ISSN 0375-9474, http://www.sciencedirect.
com/science/article/pii/S0375947410006810
[4] F. Rami, Y. Leifels, B. de Schauenburg, A. Gobbi, B. Hong, J. P. Alard, A. Andronic,
R. Averbeck, V. Barret, Z. Basrak, N. Bastid, I. Belyaev, A. Bendarag, G. Berek,
ˇ
R. Caplar,
N. Cindro, P. Crochet, A. Devismes, P. Dupieux, M. Dˇzelalija, M. Eskef,
C. Finck, Z. Fodor, H. Folger, L. Fraysse, A. Genoux-Lubain, Y. Grigorian, Y. Grishkin, N. Herrmann, K. D. Hildenbrand, J. Kecskemeti, Y. J. Kim, P. Koczon, M. Kirejczyk, M. Korolija, R. Kotte, M. Kowalczyk, T. Kress, R. Kutsche, A. Lebedev, K. S.
Lee, V. Manko, H. Merlitz, S. Mohren, D. Moisa, J. M¨osner, W. Neubert, A. Nianine,
D. Pelte, M. Petrovici, C. Pinkenburg, C. Plettner, W. Reisdorf, J. Ritman, D. Sch¨
ull,
Z. Seres, B. Sikora, K. S. Sim, V. Simion, K. Siwek-Wilczy´
nska, A. Somov, M. R.
Stockmeier, G. Stoicea, M. Vasiliev, P. Wagner, K. Wi´sniewski, D. Wohlfarth, J. T.
Yang, I. Yushmanov, A. Zhilin, and (FOPI Collaboration), “Isospin tracing: A probe of
nonequilibrium in central heavy-ion collisions,” Physical Review Letters 84, 1120–1123
(Feb. 2000), http://link.aps.org/doi/10.1103/PhysRevLett.84.1120

114

[5] Roy A. Lacey, N. N. Ajitanand, J. M. Alexander, P. Chung, W. G. Holzmann, M. Issah,
A. Taranenko, P. Danielewicz, and Horst St¨ocker, “Has the QCD critical point been
signaled by observations at the BNL relativistic heavy ion collider?.” Physical Review
Letters 98, 092301 (Mar. 2007), http://link.aps.org/doi/10.1103/PhysRevLett.
98.092301
[6] David Brown, “An introduction to source imaging in heavy-ion reactions,” Website (May
2002), https://adg.llnl.gov/Research/source_imaging/intro.html
[7] Micha A. Kilburn, Proton-Proton Correlation Functions as a Probe to Reaction Dynamics, Ph.D., Michigan State University, East Lansing, MI, USA (2011), http:
//www.nscl.msu.edu/ourlab/publications/download/Kilburn2011_281.pdf
[8] P. Danielewicz and G.F. Bertsch, “Production of deuterons and pions in a transport model of energetic heavy-ion reactions,” Nuclear Physics A 533, 712–748
(Nov. 1991), ISSN 0375-9474, http://www.sciencedirect.com/science/article/
pii/037594749190541D
[9] P. Danielewicz, “Formation of composites emitted at large angles in intermediate and
high energy reactions,” Nuclear Physics A 545, 21–34 (Aug. 1992), ISSN 0375-9474,
http://www.sciencedirect.com/science/article/pii/037594749290443N
[10] Pawel Danielewicz and Qiubao Pan, “Blast of light fragments from central heavy-ion
collisions,” Physical Review C 46, 2002 (Nov. 1992), http://link.aps.org/doi/10.
1103/PhysRevC.46.2002
[11] Qiubao Pan and Pawel Danielewicz, “From sideward flow to nuclear compressibility,”
Physical Review Letters 70, 2062–2065 (Apr. 1993), http://link.aps.org/doi/10.
1103/PhysRevLett.70.2062
[12] P. Danielewicz, “Effects of compression and collective expansion on particle emission
from central heavy-ion reactions,” Physical Review C 51, 716–750 (Feb. 1995), http:
//link.aps.org/doi/10.1103/PhysRevC.51.716
[13] P. Danielewicz, Roy A. Lacey, P.-B. Gossiaux, C. Pinkenburg, P. Chung, J. M.
Alexander, and R. L. McGrath, “Disappearance of elliptic flow: A new probe for
the nuclear equation of state,” Physical Review Letters 81, 2438–2441 (Sep. 1998),
http://link.aps.org/doi/10.1103/PhysRevLett.81.2438
[14] P. K. Kovtun, D. T. Son, and A. O. Starinets, “Viscosity in strongly interacting quantum
field theories from black hole physics,” Physical Review Letters 94, 111601 (Mar. 2005),
http://link.aps.org/doi/10.1103/PhysRevLett.94.111601
115

[15] N. Auerbach and S. Shlomo, “η/s ratio in finite nuclei,” Physical Review Letters 103,
172501 (Oct. 2009), http://link.aps.org/doi/10.1103/PhysRevLett.103.172501
[16] Subrata Pal, “Shear viscosity to entropy density ratio in nuclear multifragmentation,” Physical Review C 81, 051601 (May 2010), http://link.aps.org/doi/10.1103/
PhysRevC.81.051601
[17] Chenglong Zhou, Yugang Ma, and Deqing Fang, “Shear viscosity to entropy density
ratio in Au+Au central collisions,” Plasma Science and Technology 14, 585–587 (Jul.
2012), ISSN 1009-0630, http://iopscience.iop.org/1009-0630/14/7/04
[18] Cai Yanhuang and M. Di Toro, “Semiclassical description of isovector giant multipole
resonances,” Physical Review C 39, 105–113 (Jan. 1989), http://link.aps.org/doi/
10.1103/PhysRevC.39.105
[19] T. Gaitanos, M. Colonna, M. Di Toro, and H.H. Wolter, “Stopping and isospin equilibration in heavy ion collisions,” Physics Letters B 595, 209–215 (Aug. 2004), ISSN 03702693, http://www.sciencedirect.com/science/article/pii/S0370269304008792
[20] V. R. Pandharipande and Steven C. Pieper, “Nuclear transparency to intermediateenergy nucleons from (e,e’p) reactions,” Physical Review C 45, 791–798 (Feb. 1992),
http://link.aps.org/doi/10.1103/PhysRevC.45.791
[21] Declan Persram and Charles Gale, “Elliptic flow in intermediate energy heavy ion
collisions and in-medium effects,” Physical Review C 65, 064611 (May 2002), http:
//link.aps.org/doi/10.1103/PhysRevC.65.064611
[22] Bao-An Li and Lie-Wen Chen, “Nucleon-nucleon cross sections in neutron-rich matter
and isospin transport in heavy-ion reactions at intermediate energies,” Physical Review
C 72, 064611 (Dec. 2005), http://link.aps.org/doi/10.1103/PhysRevC.72.064611
[23] P. Danielewicz, “Hadronic transport models,” Acta Physica Polonica B 33, 45 (2002),
acta Phys.Polon. B33 (2002) 45-64, http://arxiv.org/abs/nucl-th/0201032
[24] T. Alm, G. R¨opke, and M. Schmidt, “Critical enhancement of the in-medium nucleonnucleon cross section at low temperatures,” Physical Review C 50, 31–37 (Jul. 1994),
http://link.aps.org/doi/10.1103/PhysRevC.50.31
[25] T. Alm, G. R¨opke, W. Bauer, F. Daffin, and M. Schmidt, “The in-medium nucleonnucleon cross section and BUU simulations of heavy-ion reactions,” Nuclear Physics
A 587, 815–827 (May 1995), ISSN 0375-9474, http://www.sciencedirect.com/
science/article/pii/037594749500026W
116

[26] C. Fuchs, Amand Faessler, and M. El-Shabshiry, “Off-shell behavior of the in-medium
nucleon-nucleon cross section,” Physical Review C 64, 024003 (Jul. 2001), http://link.
aps.org/doi/10.1103/PhysRevC.64.024003
[27] W. J. Nicholson and I. Halpern, “Direct-interaction effects in medium-energy fission of
uranium,” Physical Review 116, 175–177 (Oct. 1959), http://link.aps.org/doi/10.
1103/PhysRev.116.175
[28] Torbjørn Sikkeland, Eldon L. Haines, and Victor E. Viola, “Momentum transfer in
heavy-ion-induced fission,” Physical Review 125, 1350–1357 (Feb. 1962), http://link.
aps.org/doi/10.1103/PhysRev.125.1350
[29] V. E. Viola, B. B. Back, K. L. Wolf, T. C. Awes, C. K. Gelbke, and H. Breuer, “Linear
momentum transfer in nonrelativistic nucleus-nucleus collisions,” Physical Review C 26,
178–188 (Jul. 1982), http://link.aps.org/doi/10.1103/PhysRevC.26.178
[30] E. Colin, Rulin Sun, N. N. Ajitanand, John M. Alexander, M. A. Barton, P. A. DeYoung,
A. Elmaani, C. J. Gelderloos, E. E. Gualtieri, D. Guinet, S. Hannuschke, J. A. Jasma,
L. Kowalski, Roy A. Lacey, J. Lauret, E. Norbeck, R. Pak, G. F. Peaslee, M. Stern,
N. T. B. Stone, S. D. Sundbeck, A. M. Vander Molen, G. D. Westfall, and J. Yee,
“Splintering central collisions: Systematics of momentum and energy deposition for
(17–115)A MeV 40 ar,” Physical Review C 57, R1032–R1036 (Mar. 1998), http://
link.aps.org/doi/10.1103/PhysRevC.57.R1032
[31] W. Reisdorf, M. Stockmeier, A. Andronic, M.L. Benabderrahmane, O.N. Hartmann,
N. Herrmann, K.D. Hildenbrand, Y.J. Kim, M. Kiˇs, P. Koczo´
n, T. Kress, Y. Leifels,
ˇ
X. Lopez, M. Merschmeyer, A. Sch¨
uttauf, V. Barret, Z. Basrak, N. Bastid, R. Caplar,
P. Crochet, P. Dupieux, M. Dˇzelalija, Z. Fodor, Y. Grishkin, B. Hong, T.I. Kang,
J. Kecskemeti, M. Kirejczyk, M. Korolija, R. Kotte, A. Lebedev, T. Matulewicz,
W. Neubert, M. Petrovici, F. Rami, M.S. Ryu, Z. Seres, B. Sikora, K.S. Sim, V. Simion,
K. Siwek-Wilczy´
nska, V. Smolyankin, G. Stoicea, Z. Tymi´
nski, K. Wi´sniewski, D. Wohlfarth, Z.G. Xiao, H.S. Xu, I. Yushmanov, and A. Zhilin, “Systematics of pion
emission in heavy ion collisions in the regime,” Nuclear Physics A 781, 459–508
(Jan. 2007), ISSN 0375-9474, http://www.sciencedirect.com/science/article/
pii/S0375947406007676
[32] Particle Data Group, J. Beringer, J. F. Arguin, R. M. Barnett, K. Copic, O. Dahl,
D. E. Groom, C. J. Lin, J. Lys, H. Murayama, C. G. Wohl, W. M. Yao, P. A. Zyla,
C. Amsler, M. Antonelli, D. M. Asner, H. Baer, H. R. Band, T. Basaglia, C. W. Bauer,
J. J. Beatty, V. I. Belousov, E. Bergren, G. Bernardi, W. Bertl, S. Bethke, H. Bichsel, O. Biebel, E. Blucher, S. Blusk, G. Brooijmans, O. Buchmueller, R. N. Cahn,
M. Carena, A. Ceccucci, D. Chakraborty, M. C. Chen, R. S. Chivukula, G. Cowan,
G. D’Ambrosio, T. Damour, D. de Florian, A. de Gouvˆea, T. DeGrand, P. de Jong,
117

G. Dissertori, B. Dobrescu, M. Doser, M. Drees, D. A. Edwards, S. Eidelman, J. Erler,
V. V. Ezhela, W. Fetscher, B. D. Fields, B. Foster, T. K. Gaisser, L. Garren, H. J.
Gerber, G. Gerbier, T. Gherghetta, S. Golwala, M. Goodman, C. Grab, A. V. Gritsan, J. F. Grivaz, M. Gr¨
unewald, A. Gurtu, T. Gutsche, H. E. Haber, K. Hagiwara,
C. Hagmann, C. Hanhart, S. Hashimoto, K. G. Hayes, M. Heffner, B. Heltsley, J. J.
Hern´andez-Rey, K. Hikasa, A. H¨ocker, J. Holder, A. Holtkamp, J. Huston, J. D. Jackson, K. F. Johnson, T. Junk, D. Karlen, D. Kirkby, S. R. Klein, E. Klempt, R. V.
Kowalewski, F. Krauss, M. Kreps, B. Krusche, Yu. V. Kuyanov, Y. Kwon, O. Lahav, J. Laiho, P. Langacker, A. Liddle, Z. Ligeti, T. M. Liss, L. Littenberg, K. S.
Lugovsky, S. B. Lugovsky, T. Mannel, A. V. Manohar, W. J. Marciano, A. D. Martin,
A. Masoni, J. Matthews, D. Milstead, R. Miquel, K. M¨onig, F. Moortgat, K. Nakamura, M. Narain, P. Nason, S. Navas, M. Neubert, P. Nevski, Y. Nir, K. A. Olive,
L. Pape, J. Parsons, C. Patrignani, J. A. Peacock, S. T. Petcov, A. Piepke, A. Pomarol,
G. Punzi, A. Quadt, S. Raby, G. Raffelt, B. N. Ratcliff, P. Richardson, S. Roesler,
S. Rolli, A. Romaniouk, L. J. Rosenberg, J. L. Rosner, C. T. Sachrajda, Y. Sakai, G. P.
Salam, S. Sarkar, F. Sauli, O. Schneider, K. Scholberg, D. Scott, W. G. Seligman, M. H.
Shaevitz, S. R. Sharpe, M. Silari, T. Sj¨ostrand, P. Skands, J. G. Smith, G. F. Smoot,
S. Spanier, H. Spieler, A. Stahl, T. Stanev, S. L. Stone, T. Sumiyoshi, M. J. Syphers,
F. Takahashi, M. Tanabashi, J. Terning, M. Titov, N. P. Tkachenko, N. A. T¨ornqvist,
D. Tovey, G. Valencia, K. van Bibber, G. Venanzoni, M. G. Vincter, P. Vogel, A. Vogt,
W. Walkowiak, C. W. Walter, D. R. Ward, T. Watari, G. Weiglein, E. J. Weinberg,
L. R. Wiencke, L. Wolfenstein, J. Womersley, C. L. Woody, R. L. Workman, A. Yamamoto, G. P. Zeller, O. V. Zenin, J. Zhang, R. Y. Zhu, G. Harper, V. S. Lugovsky,
and P. Schaffner, “Review of particle physics,” Physical Review D 86, 010001 (Jul. 2012),
http://link.aps.org/doi/10.1103/PhysRevD.86.010001
[33] A. Andronic, J. Lukasik, W. Reisdorf, and W. Trautmann, “Systematics of stopping and
flow in Au+Au collisions,” in Dynamics and Thermodynamics with Nuclear Degrees of
Freedom, edited by Philippe Chomaz, Francesca Gulminelli, Wolfgang Trautmann, and
Sherry J. Yennello (Springer Berlin Heidelberg, 2006) pp. 31–46, ISBN 978-3-540-464969, http://www.springerlink.com/content/r254w7472431kn85/abstract/
[34] Pawel Danielewicz, Roy Lacey, and William G. Lynch, “Determination of the equation of
state of dense matter,” Science 298, 1592–1596 (Nov. 2002), ISSN 0036-8075, 1095-9203,
http://www.sciencemag.org/content/298/5598/1592
[35] J.P. Blaizot, J.F. Berger, J. Decharg´e, and M. Girod, “Microscopic and macroscopic determinations of nuclear compressibility,” Nuclear Physics A 591, 435–457
(Aug. 1995), ISSN 0375-9474, http://www.sciencedirect.com/science/article/
pii/037594749500294B
[36] W. Reisdorf, A. Andronic, R. Averbeck, M.L. Benabderrahmane, O.N. Hartmann,
N. Herrmann, K.D. Hildenbrand, T.I. Kang, Y.J. Kim, M. Kiˇs, P. Koczo´
n, T. Kress,
118

Y. Leifels, M. Merschmeyer, K. Piasecki, A. Sch¨
uttauf, M. Stockmeier, V. Barret, Z. Basˇ
rak, N. Bastid, R. Caplar, P. Crochet, P. Dupieux, M. Dˇzelalija, Z. Fodor, P. Gasik,
Y. Grishkin, B. Hong, J. Kecskemeti, M. Kirejczyk, M. Korolija, R. Kotte, A. Lebedev,
X. Lopez, T. Matulewicz, W. Neubert, M. Petrovici, F. Rami, M.S. Ryu, Z. Seres,
B. Sikora, K.S. Sim, V. Simion, K. Siwek-Wilczy´
nska, V. Smolyankin, G. Stoicea,
Z. Tymi´
nski, K. Wi´sniewski, D. Wohlfarth, Z.G. Xiao, H.S. Xu, I. Yushmanov, and
A. Zhilin, “Systematics of central heavy ion collisions in the 1 a GeV regime,” Nuclear
Physics A 848, 366–427 (Dec. 2010), ISSN 0375-9474, http://www.sciencedirect.
com/science/article/pii/S0375947410006810
[37] C. Hartnack, Rajeev K. Puri, J. Aichelin, J. Konopka, S. A. Bass, H. St¨ocker, and
W. Greiner, “Modelling the many-body dynamics of heavy ion collisions: Present status
and future perspective,” The European Physical Journal A - Hadrons and Nuclei 1, 151–
169 (Feb. 1998), ISSN 1434-6001, 1434-601X, http://link.springer.com/article/
10.1007/s100500050045
[38] Pawel Danielewicz, Brent Barker, and Lijun Shi, “From stopping to viscosity in nuclear
reactions,” AIP Conference Proceedings 1128, 104–111 (May 2009), ISSN 0094243X,
http://proceedings.aip.org/resource/2/apcpcs/1128/1/104_1
[39] Roy A. Lacey, “Is there a sonic boom in the little bang at RHIC?.” Nuclear Physics
A 785, 122–127 (Mar. 2007), ISSN 0375-9474, http://www.sciencedirect.com/
science/article/pii/S0375947406009109
[40] F. Karsch, E. Laermann, and A. Peikert, “Quark mass and flavour dependence of the
QCD phase transition,” Nuclear Physics B 605, 579–599 (Jul. 2001), ISSN 0550-3213,
http://www.sciencedirect.com/science/article/pii/S0550321301002000
[41] R. Hanbury Brown and R.Q. Twiss, “A new type of interferometer for use in radio
astronomy,” Philosophical Magazine Series 7 45, 663–682 (1954), ISSN 1941-5982, http:
//www.tandfonline.com/doi/abs/10.1080/14786440708520475
[42] Gerson Goldhaber, Sulamith Goldhaber, Wonyong Lee, and Abraham Pais, “Influence of
bose-einstein statistics on the antiproton-proton annihilation process,” Physical Review
120, 300–312 (Oct. 1960), http://link.aps.org/doi/10.1103/PhysRev.120.300
[43] Michael Annan Lisa, Scott Pratt, Ron Soltz, and Urs Wiedemann, “FEMTOSCOPY IN
RELATIVISTIC HEAVY ION COLLISIONS: two decades of progress,” Annual Review
of Nuclear and Particle Science 55, 357–402 (2005), http://www.annualreviews.org/
doi/abs/10.1146/annurev.nucl.55.090704.151533

119

[44] A. Perrin, R. B¨
ucker, S. Manz, T. Betz, C. Koller, T. Plisson, T. Schumm, and
J. Schmiedmayer, “Hanbury brown and twiss correlations across the bose-einstein condensation threshold,” Nature Physics 8, 195–198 (2012), ISSN 1745-2473, http://www.
nature.com.proxy1.cl.msu.edu/nphys/journal/v8/n3/full/nphys2212.html
[45] D. O. Handzy, M. A. Lisa, C. K. Gelbke, W. Bauer, F. C. Daffin, P. Decowski, W. G.
Gong, E. Gualtieri, S. Hannuschke, R. Lacey, T. Li, W. G. Lynch, C. M. Mader,
G. F. Peaslee, T. Reposeur, S. Pratt, A. M. Vander Molen, G. D. Westfall, J. Yee,
and S. J. Yennello, “Two-proton correlation functions for ˆ{36}ar+ˆ{45}sc at E/A=80
MeV,” Physical Review C 50, 858–870 (Aug. 1994), http://link.aps.org/doi/10.
1103/PhysRevC.50.858
[46] D. Ardouin, “Recent light particle correlation data from heavy ion collisions at intermediate and low energies,” International Journal of Modern Physics E 06, 391–435 (Sep.
1997), ISSN 0218-3013, 1793-6608, http://www.worldscientific.com/doi/abs/10.
1142/S021830139700024X
[47] E. De Filippo, F. Amorini, L. Auditore, V. Baran, I. Berceanu, G. Cardella, M. Colonna,
E. Geraci, S. Gian`ı, L. Grassi, A. Grzeszczuk, P. Guazzoni, J. Han, E. La Guidara,
G. Lanzalone, I. Lombardo, C. Maiolino, T. Minniti, A. Pagano, M. Papa, E. Piasecki,
S. Pirrone, G. Politi, A. Pop, F. Porto, F. Rizzo, P. Russotto, S. Santoro, A. Trifir`o,
M. Trimarchi, G. Verde, M. Vigilante, J. Wilczy´
nski, and L. Zetta,“Correlations between
isospin dynamics and intermediate mass fragments emission time scales: a probe for the
symmetry energy in asymmetric nuclear matter,” arXiv:1209.6461(Sep. 2012), http:
//arxiv.org/abs/1209.6461
[48] M. A. Lisa, C. K. Gelbke, P. Decowski, W. G. Gong, E. Gualtieri, S. Hannuschke,
R. Lacey, T. Li, W. G. Lynch, G. F. Peaslee, S. Pratt, T. Reposeur, A. M. Vander Molen, G. D. Westfall, J. Yee, and S. J. Yennello, “Observation of lifetime effects
in two-proton correlations for well-characterized sources,” Physical Review Letters 71,
2863–2866 (Nov. 1993), http://link.aps.org/doi/10.1103/PhysRevLett.71.2863
[49] R. Ghetti and J. Helgesson,“Isospin effects in two-particle correlation functions: A probe
for the isospin dependence of the nuclear equation of state,” Nuclear Physics A 752,
480–483 (Apr. 2005), ISSN 0375-9474, http://www.sciencedirect.com/science/
article/pii/S0375947405002034
[50] Steven E. Koonin, “Proton pictures of high-energy nuclear collisions,” Physics Letters
B 70, 43–47 (Sep. 1977), ISSN 0370-2693, http://www.sciencedirect.com/science/
article/pii/0370269377903409

120

[51] S. Pratt, T. Cs¨org˝o, and J. Zim´anyi, “Detailed predictions for two-pion correlations in
ultrarelativistic heavy-ion collisions,” Physical Review C 42, 2646–2652 (Dec. 1990),
http://link.aps.org/doi/10.1103/PhysRevC.42.2646
[52] D. A. Brown and P. Danielewicz, “Observing non-gaussian sources in heavy-ion reactions,” Physical Review C 64, 014902 (Jun. 2001), http://link.aps.org/doi/10.
1103/PhysRevC.64.014902
[53] Michael Youngs, personal communication (Nov. 2012)
[54] V. Henzl, M. A. Kilburn, Z. Chajecki, D. Henzlova, W. G. Lynch, D. Brown, A. Chbihi,
D. D. S. Coupland, P. Danielewicz, R. T. deSouza, M. Famiano, C. Herlitzius, S. Hudan,
Jenny Lee, S. Lukyanov, A. M. Rogers, A. Sanetullaev, L. G. Sobotka, Z. Y. Sun, M. B.
Tsang, A. Vander Molen, G. Verde, M. S. Wallace, and M. Youngs, “Angular dependence
in proton-proton correlation functions in central ˆ{40}ca + ˆ{40}ca and ˆ{48}ca +
ˆ{48}ca reactions,” Physical Review C 85, 014606 (Jan. 2012), http://link.aps.org/
doi/10.1103/PhysRevC.85.014606
[55] Paul Chung, “3D pion & kaon source imaging from run4 and run7 200 AGeV
Au+Au collisions,” (2010), http://wpcf2010.bitp.kiev.ua/presentations/sept18/
Kiev10_Chung_STAR_PionKaonSourceImaging.ppt
[56] P. Chung, A. Taranenko, R. Lacey, W. Holzmann, J. Alexander, and M. Issah, “Extraction of source images in Au+Au and d+Au collisions at RHIC,” Nuclear Physics A 749,
275–278 (Mar. 2005), ISSN 0375-9474, http://www.sciencedirect.com/science/
article/pii/S0375947404012680
[57] F. James and M. Roos, “Minuit - a system for function minimization and analysis of the parameter errors and correlations,” Computer Physics Communications
10, 343–367 (Dec. 1975), ISSN 0010-4655, http://www.sciencedirect.com/science/
article/pii/0010465575900399
[58] P. Danielewicz, “Quantum theory of nonequilibrium processes, i,” Annals of Physics
(N.Y.) 152, 239 (1984)
[59] G. Verde, P. Danielewicz, D. A. Brown, W. G. Lynch, C. K. Gelbke, and M. B. Tsang,
“Probing transport theories via two-proton source imaging,” Physical Review C 67,
034606 (Mar. 2003), http://link.aps.org/doi/10.1103/PhysRevC.67.034606
[60] L. P. Kadanoff and G. Baym, Quantum Statistical Mechanics (Benjamin, New York,
1962)

121

[61] Wim Botermans and Rudi Malfliet, “Quantum transport theory of nuclear matter,” Physics Reports 198, 115–194 (Dec. 1990), ISSN 0370-1573, http://www.
sciencedirect.com/science/article/pii/037015739090174Z
[62] H. S. K¨ohler, “Memory and correlation effects in nuclear collisions,” Physical Review C
51, 3232–3239 (Jun. 1995), http://link.aps.org/doi/10.1103/PhysRevC.51.3232
[63] N. H. Kwong, M. Bonitz, R. Binder, and H. S. K¨ohler, “Semiconductor
kadanoff-baym equation results for optically excited Electron–Hole plasmas in
quantum wells,” physica status solidi (b) 206, 197–203 (1998), ISSN 15213951, http://onlinelibrary.wiley.com/doi/10.1002/(SICI)1521-3951(199803)
206:1<197::AID-PSSB197>3.0.CO;2-9/abstract
[64] Nils Erik Dahlen and Robert van Leeuwen, “Solving the kadanoff-baym equations for
inhomogeneous systems: Application to atoms and molecules,” Physical Review Letters 98, 153004 (Apr. 2007), http://link.aps.org/doi/10.1103/PhysRevLett.98.
153004
[65] K. Balzer, S. Bauch, and M. Bonitz, “Time-dependent second-order born calculations
for model atoms and molecules in strong laser fields,” Physical Review A 82, 033427
(Sep. 2010), http://link.aps.org/doi/10.1103/PhysRevA.82.033427
[66] K. Balzer, M. Bonitz, R. van Leeuwen, A. Stan, and N. E. Dahlen, “Nonequilibrium
green’s function approach to strongly correlated few-electron quantum dots,” Physical
Review B 79, 245306 (Jun. 2009), http://link.aps.org/doi/10.1103/PhysRevB.79.
245306
[67] K. Balzer, S. Hermanns, and M. Bonitz, “Electronic double excitations in quantum
wells: Solving the two-time kadanoff-baym equations,” EPL (Europhysics Letters) 98,
67002 (Jun. 2012), ISSN 0295-5075, http://iopscience.iop.org.proxy1.cl.msu.
edu/0295-5075/98/6/67002
[68] N. S¨akkinen, M. Manninen, and R. van Leeuwen, “The Kadanoff–Baym approach to
double excitations in finite systems,” New Journal of Physics 14, 013032 (Jan. 2012),
ISSN 1367-2630, http://iopscience.iop.org.proxy1.cl.msu.edu/1367-2630/14/
1/013032
[69] T. Gasenzer, “Ultracold gases far from equilibrium,” The European Physical
Journal Special Topics 168, 89–148 (Mar. 2009), ISSN 1951-6355, 1951-6401,
http://epjst.epj.org/index.php?option=com_article&access=doi&doi=10.
1140/epjst/e2009-00960-5

122

[70] Arnau Rios, Brent Barker, Mark Buchler, and Pawel Danielewicz, “Towards a
nonequilibrium green’s function description of nuclear reactions: One-dimensional
mean-field dynamics,” Annals of Physics 326, 1274–1319 (May 2011), ISSN 00034916, arXiv:1009.0215, http://www.sciencedirect.com/science/article/pii/
S000349161000223X
[71] P Danielewicz, “Quantum theory of nonequilibrium processes II. application to nuclear
collisions,” Annals of Physics 152, 305–326 (Feb. 1984), ISSN 0003-4916, http://www.
sciencedirect.com/science/article/pii/0003491684900939
[72] M. Tohyama, “Application of quantum theory of particle collisions to ˆ{16}o + ˆ{16}o
reactions,” Physical Review C 36, 187–192 (Jul. 1987), http://link.aps.org/doi/10.
1103/PhysRevC.36.187
[73] J. W. Negele, “The mean-field theory of nuclear structure and dynamics,” Reviews
of Modern Physics 54, 913–1015 (Oct. 1982), http://link.aps.org/doi/10.1103/
RevModPhys.54.913
[74] C´edric Golabek and C´edric Simenel, “Collision dynamics of two ˆ{238}u atomic nuclei,”
Physical Review Letters 103, 042701 (Jul. 2009), http://link.aps.org/doi/10.1103/
PhysRevLett.103.042701
[75] S.K. Bogner, R.J. Furnstahl, and A. Schwenk, “From low-momentum interactions to nuclear structure,” Progress in Particle and Nuclear Physics 65, 94–147
(Jul. 2010), ISSN 0146-6410, http://www.sciencedirect.com/science/article/
pii/S0146641010000347
[76] G. Hagen, T. Papenbrock, D. J. Dean, and M. Hjorth-Jensen, “Medium-mass nuclei
from chiral nucleon-nucleon interactions,” Physical Review Letters 101, 092502 (Aug.
2008), http://link.aps.org/doi/10.1103/PhysRevLett.101.092502
[77] W.H. Dickhoff and C. Barbieri, “Self-consistent green’s function method for nuclei and nuclear matter,” Progress in Particle and Nuclear Physics 52, 377–496
(Apr. 2004), ISSN 0146-6410, http://www.sciencedirect.com/science/article/
pii/S0146641004000535
[78] H.S. K¨ohler, “TDHF with two-body dissipation,” Nuclear Physics A 343, 315–
332 (Jul. 1980), ISSN 0375-9474, http://www.sciencedirect.com/science/article/
pii/0375947480906557
[79] P. Grang´e, J. Richert, G. Wolschin, and H.A. Weidenm¨
uller, “Influence of a collision term on finite-size one-dimensional TDHF,” Nuclear Physics A 356, 260–268
123

(Mar. 1981), ISSN 0375-9474, http://www.sciencedirect.com/science/article/
pii/0375947481901263
[80] M. Tohyama, “One-dimensional solution of an extended time-dependent hartree-fock
theory,” Physics Letters B 163, 14–20 (Nov. 1985), ISSN 0370-2693, http://www.
sciencedirect.com/science/article/pii/0370269385901820
[81] H.S. K¨ohler, N.H. Kwong, and Hashim A. Yousif, “A fortran code for solving the
Kadanoff–Baym equations for a homogeneous fermion system,” Computer Physics Communications 123, 123–142 (Dec. 1999), ISSN 0010-4655, http://www.sciencedirect.
com/science/article/pii/S001046559900260X
[82] M.D Feit, J.A Fleck Jr., and A Steiger, “Solution of the schr¨odinger equation by a spectral method,” Journal of Computational Physics 47, 412–433 (Sep. 1982), ISSN 00219991, http://www.sciencedirect.com/science/article/pii/0021999182900912
[83] M. Frigo and S.G. Johnson, “The design and implementation of FFTW3,” Proceedings
of the IEEE 93, 216–231 (2005), ISSN 0018-9219
[84] Albert Messiah, Quantum Mechanics (Dover, New York, 1999) ISBN 0-486-40924-4
[85] D.J. Thouless and J.G. Valatin, “Time-dependent hartree-fock equations and rotational
states of nuclei,” Nuclear Physics 31, 211–230 (Mar. 1962), ISSN 0029-5582, http:
//www.sciencedirect.com/science/article/pii/0029558262907411
[86] P. Bonche, S. Koonin, and J. W. Negele, “One-dimensional nuclear dynamics in the
time-dependent hartree-fock approximation,” Physical Review C 13, 1226–1258 (Mar.
1976), http://link.aps.org/doi/10.1103/PhysRevC.13.1226
[87] H.S. K¨ohler and B.S. Nilsson, “Effect of two-body collisions on nuclear dynamics,” Nuclear Physics A 417, 541–563 (Apr. 1984), ISSN 0375-9474, http://www.
sciencedirect.com/science/article/pii/0375947484904123
[88] Franco Dalfovo, Stefano Giorgini, Lev P. Pitaevskii, and Sandro Stringari, “Theory of
bose-einstein condensation in trapped gases,” Reviews of Modern Physics 71, 463–512
(Apr. 1999), http://link.aps.org/doi/10.1103/RevModPhys.71.463
[89] A. Mu˜
noz Mateo, V. Delgado, and Boris A. Malomed, “Gap solitons in elongated geometries: The one-dimensional gross-pitaevskii equation and beyond,” Physical Review
A 83, 053610 (May 2011), http://link.aps.org/doi/10.1103/PhysRevA.83.053610

124

[90] H. Ott, J. Fortagh, G. Schlotterbeck, A. Grossmann, and C. Zimmermann, “Boseeinstein condensation in a surface microtrap,” Physical Review Letters 87, 230401 (Nov.
2001), http://link.aps.org/doi/10.1103/PhysRevLett.87.230401
[91] Scott Pratt, Lecture Notes on Statistical Mechanics (East Lansing, MI, USA, 2007)

125