_ LIBRARY
Michigan State
University

 

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 cJClRC/DateDue.p65-p.15

 

MODELING RELATIVISTIC HEAVY ION COLLISIONS

By

Sen Cheng

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Physics and Astronomy

2002

ABSTRACT

MODELING RELATIVISTIC HEAVY ION COLLISIONS

By

Sen Cheng

Experiments with relativistic heavy-ion collisions offer the prospect of studying a range
of novel phenomena. One such phenomenon is the restoration and rebreaking of chiral
symmetry, which might lead to detectable consequences in the fluctuations of charged vs.
neutral pions. Alternative sources for isospin fluctuations such as total isospin conservation,
Bose-Einstein symmetrization, and hadronic resonance decays are examined in this work.
Exact and calculable expressions are derived for multiplicity distributions and isospin fluc-
tuations within a canonical ensemble, in which additive quantum numbers and total isospin

are strictly conserved.

Most attention in the research of relativistic heavy-ion collisions is directed toward
creating and studying the quark-gluon plasma, a new state of matter predicted by lattice
Quantum Chromodynamics calculations. Discerning the consequences of the phase tran-
sition from the collision debris requires extensive numerical modeling. To this end, a new
framework for modeling relativistic heavy-ion collisions, Gromit, is developed in this work
on the basis of the Boltzmann transport equation. Gromit overcomes the limitations of
previous numerical descriptions through its modular construction and improved collision

finding meshes and collision algorithms.

The new framework is employed to study some general properties of microscopic trans-

port models such as the role of viscosities due to finite-range interactions and the sensitivi—

ties of observables like elliptic flow and particle spectra to the sampling factor A. Gromit is
also used to explore the influence of secondary hadronic interactions on balance functions,
a new observable proposed as a possible signal of a new state of matter. As previous defini-
tions of balance functions were applicable only for vanishing net charge, a new formulation

for extracting balance functions in the presence of net charges is introduced and studied.

Furthermore, a new technique is presented for calculating two-pion correlations from
microscopic models. The technique accounts for quantum corrections due to pions orig-
inating from resonance decays. Considering a simple thermal model, the importance of

such quantum corrections is quantified by comparing to semi-classical methods.

To Vini and my family.

iv

ACKNOWLEDGMENTS

This thesis could only be completed with the help of many people and organizations,
to whom I am very much indebted. In particular, I would like to thank the “Studienstiftung
des deutschen Volkes” and Wolfgang Bauer for the initial financial support that brought
me to Michigan State University; the National Superconducting Cyclotron Laboratory for
providing an excellent workplace; Marc-Jan van Goethem for proof-reading my thesis; my
adviser Scott Pratt for academic and professional support, particularly during the last year
when times were rough for both of us; and Malini Vinita Samarasinghe and my family for

much and often needed moral support, which kept me going.

Table of Contents

List of Tables ix
List of Figures x
1 Introduction 1
2 Isospin Fluctuations in Ensembles with Exact Quantum Number and Isospin

Conservation 10
2.1 Chiral Symmetry in Nuclear Physics ..................... 12
2.2 Recursion Relations for Partition Functions ................. 15
2.2.1 Nondegenerate Systems ....................... 15
2.2.2 Degenerate Systems ......................... 17
2.2.3 Monte Carlo Algorithm for Particle Production ........... 19
2.3 Recursion Relations for Multiplicity Distributions and Isospin Fluctuations 21
2.3.1 Nondegenerate Systems ....................... 21
2.3.2 Degenerate Systems ......................... 23
2.3.3 Single-level Multiplicity Distribution Conserving Total Isospin . . . 25
2.3.4 Isospin Fluctuations ......................... 27
2.4 Numerical Results and Comparisons ..................... 29
2.4.1 Total Isospin Conservation ...................... 30
2.4.2 Including Symmetrization and Resonances ............. 32
The Microscopic Evolution Model Gromit 37
3.1 Transport Equations .............................. 38
3.1.1 Monte Carlo Solution ......................... 40
3.1.2 Scaling Properties of the Boltzmann Equation ............ 41
3.2 The Framework ................................ 43
3.2. l Hadrons ................................ 44
3.2.2 Initial Particle Distribution ...................... 44

vi

3.2.3 Resonances ..............................
3.2.4 Cross Sections ............................
3.2.5 Collision Handling ..........................

4 Results and Tests of Microscopic Models
4.1 Particle Spectra ................................
4.2 Elliptic Flow .................................
4.2.1 Sensitivity to the Sampling Factor A .................
4.2.2 Comparison to Experimental Data ..................
4.2.3 Elliptic Flow in the Initial State ...................
4.3 Regeneration of p Mesons ..........................

5 Viscosities from Finite-Range Interactions
5.1 Connecting Viscosities to Finite-Range Interactions .............
5.1.1 The Role of the Collision Kernel in Boltzmann Descriptions . . . .
5.1.2 Effective Viscosities from Finite-Range Interactions ........
5.2 Analyzing the Collision Kernel ........................
5.2.1 Viscous Heating in a Bjorken Expansion ...............
5.2.2 Comparison to Numerical Results ..................

6 Correlation Function for Particles from Resonance Decays
6.1 Derivation of Corrections to the Correlation Function ............
6.1.1 Correlations from Direct Sources ..................
6.1.2 Correlations from Resonant Sources .................
6.1.3 Correlation Functions from Monte Carlo Simulations ........
6.2 The Importance of Proposed Modification ..................
6.2.1 Formalism for Semi-Classical Models ................
6.2.2 Numerical Results ..........................

6.2.3 Limit of a Narrow Resonance ....................

7 Balance Functions — A QGP Signal?
7.1 Balance Functions for Relativistic Heavy Ion Collisions ...........
7.1.1 Theory ................................
7.1.2 Experiment ..............................
7.2 Model Calculations with Hadronic Rescattering ...............
7.2.1 Accounting for Experimental Acceptance ..............
7.3 Corrections for Net Charge ..........................

vii

46
48
55

66
68
72
73
76
78
81

87
88
88
89
93
94
98

103
105
105
107
110
113
114
115
120

122
123
123
125
127
130
132

8 Conclusions

Appendix

A Parameter Definitions for Numerical Simulations

A.l Initialization File for Pionwind Study

A.2 Initialization File for Balance Function Study ................

Bibliography

viii

140

145

146
146
148

150

List of Tables

3.1 Parameter A for NN-excitation matrix elements ................ 52
7.1 Summary of balance function widths without acceptance cuts. ....... 128
7.2 Summary of balance function widths with STAR acceptance cuts. ..... 131

ix

List of Figures

2.1 Potential of the linear sigma model. ..................... 14
2.2 Probability of observing no neutral pions for a system of 12 pions, sym-

metrization and resonances are neglected. .................. 31
2.3 Fluctuation 02 for nondegenerate systems as a function of system size. . . . 32
2.4 Probability of producing no neutral pions in a system of A = 24 pions at

T = 125 MeV for ensembles restricted to an overall isosinglet and with

symmetrization included ............................ 34
2.5 Isospin fluctuations as a function of density scaled by the width of a random

distribution, with symmetrization and isospin conservation .......... 35
3.1 Illustration of attractive vs. repulsive scattering. ............... 63
4.1 Transverse Mass Spectra from Gromit-r are displayed for four cross sections. 70
4.2 Spectra for pions and protons resulting from Gromit-I with 40 mb cross

sections for four sampling factors, A = l, 2, 8, 32. .............. 71
4.3 Spectra for pions and protons resulting from four models run with o = 40

mb. ...................................... 72
4.4 Elliptic flow of pions and protons as a function of transverse energy for

several sampling factors, scaled cross section. ................ 74
4.5 Elliptic flow of pions and protons as a function of transverse energy for

several sampling factors, rejection algorithm. ................ 75
4.6 Relative multiplicity distribution as a function of impact parameter accord-

ing to a Glauber-type model. ......................... 77
4.7 Comparison of elliptic flow from a Gromit model to experimental data

from Au+Au at 130 A GeV measured by the STAR collaboration for some

choices of cross sections. The experimental results are shown without sys-

tematic errors. The calculation were performed with a sampling factor of

A = 32. .................................... 78
4.8 Anisotropy in the distribution of the distance vector in the initial condition

according to a Glauber—type model of a relativistic heavy-ion collision. . . 79
4.9 Elliptic flow in the final state as a function of centrality with and without

elliptic flow in the initial state. ........................ 80
4.10 Distribution of “It’s at which the p’s decayed with and without collisions. . . 83

4.11 Distribution of p masses integrated over the entire collision history with
and without collisions. ............................ 84

4.12 Transverse mass distribution for pions and combined proton and anti-protons
with and without collisions ........................... 85

5.1 Heating due to nonlocal interactions as calculated numerically from the
collision kernel ................................ 97

5.2 Mean transverse energy in an expanding system as a function of proper time. 98

5.3 Mean transverse energy in an expanding system as a function of proper
time, forward-peaked cross section. ..................... 101

6.1 Reduced correlation function for a source of p mesons, T = 150 MeV,

P = 200 MeV/c ................................. 116
6.2 Same as Figure 6.1, except that P = 800 MeV/c ................ 117
6.3 Reduced correlation function for a source of (r) mesons, T = 150 MeV,

P = 200 MeV/c ................................. 118
6.4 Same as Figure 6.3, except that P = 800 MeV/c ................ 119
7.1 Width of the electric charge balance function for gold on gold collisions at

‘ /sNN = 200 GeV as a function of the number of participants ........ 126
7.2 Pion balance function in relative rapidity from Hij ing and Hij ing + Gromit

without acceptance cuts. ........................... 128
7.3 Pion balance function in qiznv from Hij ing and Hij ing + Gromit without

acceptance cuts. ............................... 129
7.4 Pion balance function in relative rapidity from Hij ing and Hi j ing + Gromit

with STAR acceptance cuts. ......................... 130
7.5 Pion balance function in q]?nv from Hij ing and Hij ing + Gromit with

STAR acceptance cuts. ............................ 131

7.6 Pion balance function in relative pseudo-rapidity from Hij ing and Hij ing
+ Gromit with STAR acceptance cuts. .................... 132

7.7 Uncorrected balance function in relative pseudo-rapidity from Hij ing and
Hij ing + Gromit for all electric charges. .................. 133

7.8 Uncorrected balance function in relative pseudo-rapidity from Hij ing and
Hij ing + Gromit for protons and anti-protons. ............... 134

7.9 Corrected pp’ balance function in relative pseudo-rapidity from Hij ing and
Hij ing + Gromit. .............................. 137

7.10 Corrected electric charge balance function in relative pseudo-rapidity from
Hij ing and Hij ing + Gromit. ........................ 138

7.11 Comparison between corrected and original n+1t‘ balance function in rel-
ative rapidity from Hij ing. ......................... 139

7.12 Comparison between corrected and original n+rt‘ balance function in rel-
ative rapidity from Hij ing + Gromit. .................... 139

xi

Chapter 1

Introduction

All interactions in the universe can be described in terms of four fundamental interactions:
gravitation, electromagnetism, the weak force, and the strong force. Macroscopically, grav-
itation is well explained by the theory of general relativity but it is generally neglected at the
scale of elementary particles because its strength is 30 or more orders of magnitude smaller
than that of the weak force. The electromagnetic and weak forces, however, were success-
fully combined into the electro-weak force and theoretically described by an SU (2) (8 U ( 1)
gauge theory [76, 153]. The strong force is represented by Quantum Chromodynamics

(QCD), a nonAbelian SU (3) gauge theory, as interaction between colored quarks [73, 87].

In the theory of QCD, the exchange bosons, gluons, carry color charge and therefore
interact with each other. The interaction between gluons has a profound impact on the
behavior of QCD and is responsible for its peculiar properties like confinement of quarks
and gluons into hadrons, asymptotic freedom [83, 82, 84], and a running coupling constant
as that decreases with energy transfer. This property of the strong coupling is the reason
why the perturbative expansion order-by—order in as can only yield satisfactory results
for high-energy scattering experiments, such as those that first observed the t0p quark at
Fermilab [1, 2]. Perturbative QCD (pQCD) cannot be applied to describe phenomena at

length scales above 0.1 fm. This prevents pQCD from describing the binding of quarks

into hadrons or describing the condensation of gluonic and quark degrees of freedom into
the vacuum condensates. These condensates couple to the quarks that comprise hadrons
and are responsible for most of the mass of the proton and therefore most of the mass of

the universe.

A nonperturbative approach to deal with the complexity of QCD is lattice QCD, where
space is discretized into cells and partition functions are calculated by Monte Carlo sam-
pling (see for example [51]). Lattice QCD can be used to calculate the phase diagram for
quark-gluon matter, which has recently become feasible even for finite chemical potential
[68, 67]. The calculated QCD phase diagrams show a deconfinement phase transition from
hadronic matter to a plasma of free quarks and gluons. The order of this transition depends
on the degrees of freedom used in the lattice calculation and remains a theoretical uncer-
tainty [40, 98]. Originally, the quark-gluon plasma (QGP) had been postulated by extrapo-
lating the notion of asymptotic freedom to a hot and dense environment [48, 69, 46, 140].
It was thought that if the strength of the interaction between quarks decreased with smaller

distances, then quarks should form a liberated plasma at high temperatures and densities.

The critical temperature for the phase transition from ordinary hadronic matter to the
QGP is predicted to be on the order of 150—200 MeV, varying between different lattice cal-
culations. The extreme conditions necessary to create a QGP existed at around 5 us after the
big bang and might also exist within the interiors of massive stars. There is also a prospect
of creating and studying the QGP in accelerator experiments by colliding relativistic heavy
ions, which has led to several experiments at either fixed target machines, such as the Al-
ternating Gradient Synchrotron (AGS) at the Brookhaven National Laboratory (BNL) and
the Super Proton Synchrotron (SPS) at the European Organization for Nuclear Research
(CERN), or colliders such as the Relativistic Heavy Ion Collider (RI-11C) at BNL and the
planned Large Hadron Collider (LHC) at CERN.

One aspect of the potential phase transition at high temperature and/or density is the
restoration of chiral symmetry, which is expected to occur concurrently with the deconfine—
ment transition. Chiral symmetry is a spontaneously broken symmetry of QCD in the limit
of zero quark masses with the pion as the associated Goldstone boson. The spontaneous
breaking of chiral symmetry leads to a nonzero expectation value for the scalar quark—anti-
quark condensate at low temperatures and densities. This condensate is melted as chiral
symmetry is restored at high temperature and/or density. As the system cools off and chi-
ral symmetry is spontaneously broken again, all directions along the pion fields and the
scalar quark—anti-quark field are equally probable. The chosen direction at recondensation
might not coincide with the usual scalar quark—anti-quark condensate, which is referred
to as disoriented chiral condensate (DCC) [16, 17, 33, 133, 132]. As chiral symmetry is
also explicitly broken, the DCC would eventually decay, possibly leading to measurable

signatures.

Nucleus-nucleus collisions create a phase of hot and dense matter that subsequently
cools quickly. Therefore, irrespective of whether a QCD phase transition took place during
the evolution of the collision process, only ordinary photons, leptons, and hadrons escape
the collision region and reach the detectors. The challenge in relativistic heavy-ion colli-
sion experiments is to reconstruct the hot and dense phase, which lasted on the order of
10‘23 s and had an extent on the order of 10’15 m, from its remnants. Although lattice
QCD predicts the existence of a QGP, it offers little guidance on this reconstruction as it
can neither describe hadronization nor the hadronic interactions thereafter. Quantum me-
chanical treatments such as solving the Klein—Gordon or Dirac Equations are not feasible in
this case because of the large number of independent particles involved (at RHIC typically
around 6000). They also might not be required because the thermal wavelength is small

compared to the size of the system.

Semi-classical models have provided the basis to gaining insight into the collision dy-

namics of a relativistic heavy-ion collision and have proven to be essential to designing
experiments, analyzing data, and making detailed predictions. In a simplified picture, most
models for relativistic heavy-ion collisions can be fit into three major categories: hydrody-

namical or statistical models, Glauber-type models, and microscopic models.

Hydrodynamical models were employed in the 1950’s to describe the dynamics of
pp collisions because the underlying processes were not well understood and the mod-
els relied mainly on the basic physical principles of energy and momentum conservation
[66, 122, 104]. These models assume local thermal equilibrium and neglect viscosities,
surface energy, Coulomb energy, and single-particle effects, when modeling the collision
system as a nuclear fluid. The first numerical solution of hydrodynamics for collisions of
heavy ions became feasible in 1975 [11, 12]. Hydrodynamical models have evolved into
a workhorse for predicting and understanding measurements of observables at RHIC en-
ergies, e.g., [141, 102]. However, despite the maturity of these models, some challenges

remain, such as the inclusion of viscosities in the numerical treatments [112].

Glauber-type models exploit the well-understood dynamics of nucleon-nucleon colli-
sions. Nucleus-nucleus (AA) collisions are modeled by dividing the nuclei into nucleons
with some density distribution and describing the entire reaction as a superposition of bi-
nary nucleon-nucleon collisions [14, 138, 151, 41, 94]. These models provide invaluable
benchmarks to compare with experimental results since they only contain dynamics from
nucleon-nucleon collisions, in which the number of interactions between quarks and gluon
is insufficient to create the temperature and energy density needed for the phase transition

to a QGP.

Microscopic transport models describe the dynamics of a relativistic heavy-ion col-
lision by modeling its microscopic constituents and their interactions. Such models are

currently available in a wide variety of flavors, ranging from quantum molecular dynamics

[143, 23] to solutions of the Boltzmann Transport Equation (BTE) in its various incarna—
tions [120, 109, 107, 10, 62, 114, 159, 72, 154]. Some of these models are based on trans-
port theory that has a long history dating back to the works of Boltzmann and Maxwell on
the kinetic theory of gases in the 1870’s. Over time, the original Boltzmann equation has
been modified to incorporate new physical insights like relativity, quantum statistics, and
the nuclear mean field. In addition to choosing the equations that they solve, microscopic
models select the degrees of freedom, i.e. partons, hadrons, Lund string phenomenology
[15, 13], or subsets thereof; the interactions between those particles; and the numerical
solution techniques. This flexibility might account for the wide variety of currently avail-
able models that fall into this category. The main appeal of microscopic transport models
derives from their ability to describe nonequilibrium phenomena [85, 109], as opposed to

hydrodynamical models, which have to rest on the assumption of local thermal equilibrium.

It should be mentioned that some models transcend the categories mentioned above by
either following different ideas like solving the equations of motion of classical QCD color
fields [58] or by combining several models into a hybrid model, where different stages of
the relativistic heavy-ion collision are treated with different models. Some hybrid models
combine microscopic models for different degrees of freedom [144, 162] in an attempt to
choose the most appropriate description for the various stages of the collision. For example,
one approach [22, 146] employs a hydrodynamical model for the first part of the collision,
where the assumption of local thermal equilibrium seems most justified, and a microscopic

transport model for the later stages.

The theoretical and experimental study of relativistic heavy-ion collisions has yielded
a wealth of knowledge about the dynamics of a relativistic heavy-ion collision. When
viewing a relativistic heavy-ion collision at RHIC energies in the center of mass frame

(c.m.f.) the following time evolution picture emerges (see for example [113]):

In the initial stage two Lorentz-contracted nuclei head toward each other at nearly light
speed. These nuclei pass through each other within a time less than 0.1 fm/c, during which
on the order of 104 quarks, anti-quarks, and gluons interact with each other and create a
region of hot and dense matter. The system persists in this preequilibrium stage for about
0.1— 1.0 frn/c until local thermal equilibrium is established through collisions. In addition
to the deconfinement phase transition to a QGP that might take place in this stage, there

might be a concurrent restoration of chiral symmetry.

Once the equilibrium stage is reached, the system consists of either a liberated plasma
of quarks and gluons or excited hadronic matter, depending on whether there was a de-
confinement phase transition. In case a QGP was created, quarks and gluons have to be
confined into hadrons, once the system temperature falls below the critical value, which is
predicted by lattice QCD to be on the order of 150—200 MeV. The hadrons appearing at
this hadronization phase transition are formed significantly later in the collision than in the

scenario without QGP.

In either scenario, the system ends up in the hadronic stage, where hadrons decay and
interact with each other until the system is so dilute, i.e., the inter-particle spacing so large,
that the hadrons cease to interact. Chemical freeze-out is marked by the drop in the rates
of inelastic collisions which are required to maintain equilibrium in the abundances of
the various particle types. Kinetic (or thermal) freeze-out is marked by the drop in the
rate of elastic collisions necessary to maintain the thermal shapes of the local momentum
distributions. Chemical freeze-out is expected to occur simultaneously with hadronization
[127], whereas kinetic equilibrium should be maintained until the particles have their last

interactions.

The goal of this thesis is to study signatures that reveal the processes within relativistic

heavy-ion collisions and to provide means for modeling the dynamics of such collisions.

One proposed signature of the disoriented chiral condensate (DCC) is anomalous isospin
fluctuations, which would indicate a restoration and rebreaking of chiral symmetry. To
test the validity of this signature of the DCC, alternative sources for such fluctuations are
studied in Chapter 2 by considering the influence of total isospin conservation, quantum
symmetrization, and resonance decays. Exact and calculable expressions are derived for
multiplicity distributions and isospin fluctuations within a canonical ensemble, in which
additive quantum numbers and total isospin are strictly conserved. It is found that the com-
bined effect of total isospin conservation, quantum symmetrization, and resonance decays

cannot account for anomalous isospin fluctuations.

Chapter 3 of this thesis introduces Gromit, a modular framework for building micro-
scopic transport models for the study of relativistic heavy-ion collisions. Gromit was devel-
oped and written entirely within the RHIC Transport Theory Collaboration (R'ITC) [137].
The collaboration’s goal is to develop thoroughly tested and documented models that are
available to the community on an open—source basis. Model assumptions about degrees of
freedom and interactions are entirely modular and can be easily replaced to independently
test nearly every assumption made in a model. New features like mean fields, novel cross
sections, and dynamic linking to other models can be easily incorporated into the existing

framework.

The influence of numerical artifacts and arbitrary algorithmic choices in microscopic
transport models is investigated in Chapter 4. Observables like elliptic flow and particle
spectra are found to be sensitive to the sampling factor A, but fortunately, rather small values
of A. are sufficient for extracting reliable results. These findings indicate that transport
models based on the Boltzmann Equation are invalid in the cascade limit, A = l, and that

solutions converge in the Boltzmann limit, A —+ 0°.

The importance of viscosities in relativistic heavy-ion collision is explored in Chapter

5. The consequences of a nonzero interaction range are studied by varying the interaction
range in the scattering kernel, while leaving the mean free path unchanged. Finite inter-
action range is demonstrated to contribute viscous terms in a manner similar to the finite
mean free path, but with different dependencies with respect to density and temperature.
In particular, a comparison of heating derived from an analysis of the collision kernel with

the heating observed in a simplified simulation is presented.

The influence of resonance decays on pion correlation functions is studied in Chapter
6. Correlations between particles originating from a common source are generally intro-
duced through quantum statistics and contain information about the spacetime extent of
the source. The importance of quantum treatments for the extraction of correlation func-
tions from microscopic models is quantified by comparing to semi-classical forms for a
simple thermal model. It is found that quantum corrections become important when kine-
matics constrain the resonances to be off shell. An alteration of the methods for extracting
correlations from classical simulations is proposed to better account for quantum effects.
This modification can account for the quantum corrections by incorporating information

regarding the off-shell energy of the decaying resonance.

Finally, Gromit is used in Chapter 7 to explore the influence of hadronic interactions on
balance functions, a new observable proposed as a possible signal of late-stage hadroniza-
tion. Balance functions quantify on a statistical basis the relative location of charge-anti-
charge pairs, which are created together. It is found that secondary hadronic interactions
widen the WW“ balance function in relative rapidity, but not in qfiw. The wider balance
function in relative rapidity is argued to be merely a consequence of cooling due to hadronic
rescattering. It is also demonstrated that previous definitions of balance functions are ap-
plicable only for a small net charge. A new formulation for extracting balance functions in
the presence of net charges is introduced. Though not entirely efficient, the new balance

function is shown to remove the dramatic effects of the net electric charge of the protons.

The above findings on the effect of hadronic interactions on the balance function remain

valid when studied with the corrected balance function.

It should be noted that parts of Chapter 4 (namely Section 4.1 and Subsection 4.2.1)
and Chapters 5 and 6 have been previously published in a similar form [43, 45] and that
Chapter 2 has been submitted for publication [44].

Chapter 2

Isospin Fluctuations in Ensembles with
Exact Quantum Number and Isospin
Conservation

The temperatures and densities created in a relativistic heavy-ion collision, e.g. at RHIC,
are believed to be sufficiently high for the restoration of chiral symmetry. A brief descrip-
tion of this approximate symmetry in given in the next section. The subsequent rebreaking
of this symmetry as the collision region cools can then lead to the creation of a disori-
ented chiral condensate (DCC). This concept was first introduced to explain the observa-
tion of large fluctuations in the ratio of neutral to charged particles in cosmic ray events
[16, 17, 33, 133, 132]. These so-called Centauro events [105, 37], have sparked numer-
ous studies of isospin fluctuations during the last decade. It has been proposed that the
melting and subsequent recondensation of the chiral condensate could provide a dynam-
ical means for coherent pion emission where dozens of pions are emitted with the same
isospin. If N pions are confined to a single quantum state in addition to being in an isos-
inglet [92], the probability of finding no neutral pions takes a simple form in the limit of

large N = n+ + n- + no,
dN l

Zf-Zfi’ fEI’lO/N. (21)

IO

The same result can be obtained by considering a coherent state

lfi) = exph'i 71) IO), (22)

where the pion field operators are no : nz, rti 2 (it, i: iny)/\/2, and the direction of ii is
averaged over all directions. The source of the field '7'] has been proposed to be the chiral

condensate which might disorient itself in a quenching scenario.

Alternative sources for such isospin fluctuations have been studied. The dramatically
broad isospin distribution of Eq. (2.1) relies on the assumption that the emission proceeds
via a single quantum state. The inclusion of Bose-Einstein effects in the thermal emission
of pions from a nondegenerate array of states was shown to broaden the isospin distri-
bution [124] with respect to a random distribution, but not nearly as much as in Eq. (2.1).
Neglecting isospin conservation in [124] accounted for the reduced broadening of the peak.
A crude accounting for isospin conservation was suggested by considering the emission of
neutral pion pairs (2/3 n+1? and 1/3 27:0) [128], but came far short of considering the com-
plete ensemble of isoscalar states in a multi-level system. The emission of pairs through a
classical isoscalar field into nondegenerate single-particle levels, which may be considered

as an oriented chiral condensate, has been studied as well [150].

The effects of exact charge conservation in canonical ensembles were also studied in
other contexts with a projection method. Strangeness and baryon number conservation
were found to restrict strangeness productions in p'N collisions, particularly for small sys-
tems [56]. The confinement of the quark-gluon plasma to color-singlets was shown to lead
to a reduction in the number of internal degrees of freedom, which could lead to measurable

finite size effects in relativistic heavy ion collisions [64, 65, 63].

In the following, the sophistication of statistical treatments is extended by considering
the entire ensemble of isoscalar states available in a system with many single-particle levels.

Methods for calculating isospin distributions are presented, which include Bose-Einstein

11

symmetrization, the effects of resonances and the conservation of both total isospin I and its
projection M. Sample calculations are presented to illustrate the above-mentioned effects.
It is found that symmetrization effects are important for high quantum degeneracies, that
isospin conservation has little impact when the size of the domain exceeds a dozen pions,

and that resonances can strongly narrow the distribution.

2.1 Chiral Symmetry in Nuclear Physics

The notion of chiral symmetry was inspired by studies of nuclear beta decay. There, the
weak coupling constant for the hadronic vector current CV was found to coincide with the
leptonic analogue, and the coupling for the axial vector current CA was determined to differ
by only 25% from that for the leptonic counterparts. Chiral symmetry was proposed to
explain the conservation of both currents ' and can be well illustrated with the linear sigma
model which was proposed by Gell-Mann and Levy [74] as a model for nucleons before

QCD was known as the theory of the strong interaction.

The linear sigma model describes interactions between the nucleon field w as couplings
to two meson fields, the isovector pion field 1t and a hypothetical scalar sigma field 0. The

Lagrangian for the linear sigma model is given by
.- ._ _ l - l -
L = up at]; — gn(rwy51\un + wwo) + Eaundf‘n + Eduodf‘o, (2.3)

where T are the Pauli (iso)spin matrices and gn is the coupling constant for the interaction.
This Lagrangian was defined such that it is invariant under the following vector transfor-

mation

Av: v —+ e—iTG/zwz(l—iIG/2)w (2.4)

- +iTG/2

w ——> e \Dz(l+i‘CB/2)\]J. (2.5)

 

1The fact that the axial vector current is only partially conserved is ignored for the moment and introduced
later as explicit breaking of chiral symmetry.

12

Like any other symmetry, this symmetry is associated with a conserved charge or current

as stated by the Noether theorem. In this case the vector current associated with AV,
_ t“
V: 2 WHEWa (2-6)

is simply the isospin current. Another symmetry of the chiral Lagrangian is the invariance

under the axial vector transformation

AA : \u —+ ("SIB/2w z (1 -— iy510/2)\V (2.7)
q; ——> fibre/2o z (1 —— i‘y5TB/2)\I!. (2.8)
The associated axial current,
_ 1"
A2 = Writs—w, (2.9)

2

is not easily associated with a physical current because chiral symmetry is only an approx-
imate symmetry. The invariance of the Lagrangian in Eq. (2.3) under both transformations

Av and AA is referred to as chiral symmetry.

Although the Lagrangian is invariant under the transformations AV and AA, the ground
state is not. This is known as spontaneous symmetry breaking and is imposed through the
potential

V = E ((1:2 +02) —f,§)2, (2.10)

which is illustrated in Figure 2.1. The constant E determines the mass of the o meson.
The potential of chiral symmetry, which will become part of the Lagrangian, is invariant
under the transformations Av and AA. However, the ground state, so = fn, breaks chiral
symmetry. Excitations in the 0' directions require a finite amount of energy, thus giving the
o a finite mass, whereas excitations in the pion direction (orthogonal to both axis shown in
Figure 2.1) do not require energy and leaves the pions massless. The pion is, therefore, the

Goldstone boson associated with the spontaneous breaking of chiral symmetry.

13

 

 

 

In

Figure 2.1: Potential of the linear sigma model.

However, chiral symmetry is only an approximate symmetry because it is explicitly

broken by a term in the Lagrangian ignored so far,
5L = 80. (2.11)

Due to this small explicit breaking of chiral symmetry, pions acquire a finite mass, although
it is small compared to that of other hadrons. The pion mass determines the symmetry
breaking parameter 8. This small explicit breaking of chiral symmetry is also responsible

for the the nonconservation of the axial vector current.

At high temperature and/or density the potential of chiral symmetry is predicted to
change such that chiral symmetry is restored explicitly. Then, the ground state would be
located at 1: = O and 0' = 0. If the system cools off quickly, the potential of chiral symmetry
from Eq. (2.10) can be recovered with the ground state remaining in its previous location.
This state is, of course, instable and eventually chiral symmetry must be spontaneously
broken again. When this occurs all directions along the pion fields and the scalar 0' field are
equally probable and the chosen direction might not coincide with the usual ground state.

This disoriented state is referred to as disoriented chiral condensate (DCC) [133, 132]. As

14

chiral symmetry is also explicitly broken, the DCC is also instable and would eventually

decay to return the chiral condensate to its usual ground state.

This decay of the DCC might then lead to measurable signatures. One proposed signa-
ture of the disoriented chiral condensate (DCC) are anomalous isospin fluctuations, which
was discussed in the introduction to this chapter. To test the validity of this signature,
alternative sources for such fluctuations are studied in the following by considering the

influence of total isospin conservation, quantum symmetrization, and resonance decays.

2.2 Recursion Relations for Partition Functions

2.2.1 Nondegenerate Systems

When the number of available states is much larger than the number of particles, the prob-
ability for two or more particles occupying the same quantum state is small. In such a
nondegenerate system, quantum statistics can be neglected. The partition function for a
canonical ensemble of A particles conserving an additive quantum number or a vector of
such quantities Q can be written as a product of single-particle partition functions

21,9: 2 119—L. (2.12)

(End): :14) k: I
(ZVkCIk=Q)

where N is the number of particle types, vk is the occupation number of particle type k, qk is
the charge of one particle, and the particle number ak indicates how many times a particle
contributes to the main conserved quantity A. For example, if A is the number of pions
then ap = 2, since the p meson decays predominantly into two pions. The single-particle
partition function on = gk Ziexp(—e§k) / T) sums Boltzmann factors weighted by the the

spin degeneracy gk over all available single-particle levels i.

Summation over the immense number of partitions in Eq. (2.12) can be avoided by

15

rewriting the partition function as a recursion relation [42, 55, 130],

N akwk
2M): 2: A —ZA— —atQ- (11 (2.13)
k=l

If intermediate values of the partition function are stored, the computations required for
Eq. (2.13) scale linearly in both A and N, thus making it possible to quickly calculate the

canonical partition function numerically.

A partition function conserving total isospin as well as additive quantum numbers is
derived by adding a sum over all possible isospin configurations for a given partition, {vk},

and isospin weights §(I,M|{vk}) to Eq. (2.12),

N (0"!
£21m: 2 2&(1,M|{v,-}) [1%. (2.14)
(2v_,aj=A) {V1} 1‘:
<2V141=Q>

To convert this partition function into a recursion relation insert [1(2):] akv,< = 1 into Eq. (2.14),

N Vk— I v}
akcok (l) 03'
mm: 2 A 2 Elfin—{’17 §(I,M|{vj}). (2.15)
k=1 (EvjajzA) k #k I'M}
(XV/4}=Q>

 

The entire system can be broken into two subsystems, a single particle with isospin Ik and
projection rm and a remainder system with isospin I’ and projection M — mk, which are
coupled with the appropriate Clebsch-Gordan coefficients. All possible values for the total

isospin of the remainder system have to be summed over,

[(+1
2&(17Ml{vj})= 2 2€(I’.Ml{v3-})(Itmkzl’.M—n21-IIM>2. (2.16)
{V1} 1’=|M—mk|{V')}

With this modification, the summation indexes in Eq. (2.15) can be switched,

I__ Vj ajrfék
v _{v-—1 ,j:k (2.17)

16

and the partition function written as recursion relation

 

N akmk I+Ik / 2
$21,114 = 2 A 2 971-121.1214.an (Item/(:1 ,M—mkllM) . (2.18)
k=l l’:|l—lk|

Since the partition function is the trace of an isoscalar, i.e. e‘H/T

, it will not depend on
the isospin projection M. Hence, Eq. (2.18) can be further simplified by summing the RHS

over all isospin projections mk of an isospin multiplet,

 

I+lk
akmk
mp}: A 2 521-..,1'. (2.19)
1' l’=|I—lk|

where the sum over k’ includes isomultiplets, not individual particles species.

2.2.2 Degenerate Systems

Since several particles might occupy the same quantum state in a degenerate system, sym-
metrization of the wave function has to be accounted for. This discussion will be restricted
to studying Bose-Einstein particles. States with multiple particles have to be added to

Eq. (2.13), which only contains states that are occupied by zero or one particle,

co N
ak k'
ZA,M : 2 Z X it )ZA-nak,M-nmka (2-20)
n=lk=l
where the cycle diagram is

I- ~ _ ‘ k
cf. ’ = <a|e ”/7100: 231 expr—nel ’/T). (2.21)
I

Here, the state let) refers to an n-particle state of distinguishable particles and lot) is the
cyclic permutation of that state. The single-particle energy levels are 13$“ for particle type

k. A more rigorous derivation of Eq. (2.20) is given in [129].

Since the partition function for conserved total isospin QM)" is independent of the

isospin projection M, as mentioned above, a simple relation can be derived, ZAM = ZI>M QA, 1,

l7

which in turn leads to

9A,] = ZA,M:I — ZA,M:I+1- (2-22)

A second method for calculating the pion partition function constraining total isospin
is obtained by evaluating the cycle diagram cf,“ in Eq. (2.21) for its isospin content. These

new cycle diagrams are defined as

Cm E Z<B,n,ile‘H/T|B,n,i), (2.23)
B

where the sum over [3 represents a sum over all states with fixed particle number n and
isospin i. The particles are assumed to be distinguishable and B represents a cyclic permu-

tation of particles. The partition function for the pions in term of this new cycle diagram is

then
I A n l+i
£21,]:sz 2 Cut-$21-11,. (2.24)

where the new cycle diagrams C are yet to be determined. After obtaining a recursion

relation for these functions from Eq. (2.24) itself,

A-l n I+i
CA,I=AQA,I— 2 2 2 9,191-"). (2.25)
n21i=01’=|]—i|

these cycle diagrams are found to follow a simple pattern by considering a one-level system

where S2,“ is easily calculated,

C", i: n
—Cn, i: n _l
an -— C... i: o (2.26)

0, otherwise

Since resonances are more massive and have lower phase space occupations, the prob-
ability of creating several resonances in the same state can be neglected except in the limit
of extremely high densities. Therefore, resonances might be treated as independent, non-

degenerate subsystems, for which a partition function can be obtained through Eq. (2.18).

18

The partition functions of two subsystems, l and 2, can then be convoluted to obtain that
of the entire system,
A A’ I+I’ (r) (7)
9A,],M = Z 2 Z QA,’,,S2A"_A,,,H. (2.27)
A’=01’=01”=|l—l’|
If there are more than two subsystems Eq. (2.27) can be applied successively, i.e., partition

functions of any two subsystems are convoluted first to obtain a new partition function,

which is then convoluted with the partition function of another subsystem, and so on.

It should be pointed out that treating resonances as if they are in a different system is not
consistent with the indistinguishability of pions from resonances and direct pions. For nar-
row resonances, the pions could be uniquely identified by constructing the invariant masses
of the constituents. However, broad resonances cannot be identified confidently. The crite-
ria for resonances being narrow are identical to the criteria that they have sufficiently long
lifetime to decay outside the collision region. For all the calculations considered in this
study, it is assumed that the resonances are separable. This assumption is excellent for
pions from n mesons, good for pions from (n resonances and questionable for pions from p
decays. Including the effects of symmetrizin g the resonant and nonresonant pions remains

an open question.

2.2.3 Monte Carlo Algorithm for Particle Production

The partition functions with exact quantum number and isospin conservation obtained in
2.2, can be applied in a Monte Carlo algorithm for particle generation. Recently, there
has been much interest in modeling relativistic heavy ion collisions with hybrid models,
in which early, dense stages of the collision are described by a hydrodynamical model be-
fore switching to a hadronic cascade to simulate the freeze-out stage [22]. The change of
degrees of freedom at the interface between the two models from the energy-momentum

tensor to hadrons is generally modeled by a grand—canonical ensemble, which conserves

l9

charges only in the average over many events. However, event—by-event charge conserva-
tion is essential to calculating observables like fluctuations and balance functions, which

have been proposed as possible signal for the quark-gluon plasma [24, 146].

Energy and momentum conservation is generally enforced with the Cooper-Frye for-
malism [49, 50], whereas the object of this algorithm is to generate an ensemble of parti-
cles with given quantum numbers Q and total isospin I. The following algorithm is shown
to be consistent with partition functions derived in 2.2.1 for the nondegenerate case. The
algorithm can be easily adjusted to generate Bose-Einstein particles with only additive
quantum number conservation. An equivalent generation algorithm that conserves total
isospin while accounting for Bose-Einstein symmetrization is difficult to derive because of

negative contributions in Eq. (2.26) and is not attempted here.

An ensemble of particles with given quantum numbers Q and total isospin I can be
generated by following these steps: In the first step, the total number of hadrons A has to
be determined according to the probability for finding a certain hadron number in a grand-
canonical ensemble

9
WA 2 ——’1’i9——. (2.28)

2 S2A,1,Q
A20

Particles are then generated one by one in an iterative process until the number of hadrons
determined in the first step has been exhausted. In the second step, the particle type k to be
generated is obtained according to their relative weights

1+1,

= Q A 2 (Ikmk;1’,M—mkIIM)ZQA_ak‘,1,M_me_qk. (2.29)
A’I’M’Q l'=|M-Mkl

I am»,

 

Wk

The total isospin of the system I and the isospin of the particle 11- can couple to more than

one isospin of the remainder system 1’. Therefore, in the third step, that isospin I’ has to be

20

picked with the appropriate probability

(Ikmk;l’,M — mk |1M)252A_ak,,1,M_mk,Q_qk

l+lk
2 (Ikmk;I”9M-_ mkIIM)29A—ak,l”,M—mk,Q—qk
I”=lM—mkl

 

W‘Ir :

(2.30)

The last two steps are repeatedly iterated with the following substitutions

I —> 1’
M ——> M — mk
A —+ A — ak
until the number of hadrons acquired in the first step has been generated, i.e. A = 0. Con-
servation laws are ascertained at every step of this algorithm because the probabilities Eqs.

(2.28H2.30) would vanish for those configurations that violated the conservation laws,

hence, no such configuration will be encountered.

2.3 Recursion Relations for Multiplicity Distributions and
Isospin Fluctuations

2.3.1 Nondegenerate Systems

The multiplicity distribution can be calculated from a ratio of partition functions, where the
numerator includes an extra constraint,

ZA ,M.nj

 

PA,M(n,-) = (2.31)

ZA,M .
Here the numerator represents a canonical ensemble with the appropriate conservation laws
containing n j particles of type j. This additional constraint can be regarded as a “charge”
and, therefore, is added to the indices of the partition function. The partition function in

the numerator can be rewritten with the aid of Eq. (2. 13),

N
akook
ZA,M,nj : 2 A ZA—ak,M—mk,nj—dk‘ja (232)
k2]

 

21

where the feed-down factor dkdo indicates that a particle of type k decays into d,“- particles
of type j. This approach, however, will not work for a resonance that decays via more than
one decay channel. In such a case, a pseudo-particle is included for each decay branch with

its degeneracy g,- scaled by the corresponding branching ratio.

It will prove more convenient to write equations in terms of the product of partition

function and multiplicity distribution

WA,M("j) E ZA,MPA,M("j) (2.33)

instead of the multiplicity distribution itself. For a nondegenerate system conserving only

additive charges, the multiplicity distribution can be obtained from

akwk

N
WA,M("j) = Z A

kzl

 

WA—ak,M—mk (nj —dk,j)- (2.34)

The occupation number for particle type j including feed-downs from resonance decays
is determined by multiplying the occupation number of all particles, given in [130], by the

feed-down-factor and summing over all resonances,

I N
("12 : 2707,21. dk,j(1)kZA—ak,M—mk- (2.35)

The second moment of the distribution will be needed for calculating isospin fluctuations

as in Subsection 2.3.4,

N d -d .
k,j k’,’
(2.-12.): 2 ——J—

{Ska/1” meA~ak,M—mk + mkmk’ZA—ak—akl,M—mk—mk; } (2-36)
k,k’:l 2N”

The multiplicity distribution incorporating conserved total isospin can be derived in a

similar manner as was employed for Eq. (2.18),

akwk

 

N
WA,,,M(nj) = 2 2(Ikmk;l’, M — mk |1M)2WA_ak,,I,M_mk (n) — d“). (2.37)
k=l I’

22

2.3.2 Degenerate Systems

When considering only quantum numbers corresponding to additive charges, the multiplic-

ity distribution is obtained through Eqs. (2.20) and (2.31) and

0k

°° k
A 2C; )WA—Iaj,M——Inrj(nj_ldk,j)a (2.38)

121

N
WA.M("j) = 2
1:1

where C, is the cycle diagram defined in Eq. (2.21). Calculating the two-point function and

the four-point functions permits the derivation of the first two moments of the multiplicity

distribution. The 2-point function is

8. .
*L ZCXM—n8i/T)ZA—na,-,M—nm,- (239)
A,M n

(afaj) =

Summing over all particle types and states and multiplying by feed-down factors results in

an expression for the occupation numbers,

 

1 N °° (k)
-=— d.- CZ_ _,. 2.40
(I?!) Z/LM [(2:31 [1,112] 1 A lak,M Ink ( )
Similarly, the 4-point function
1 t _ II 11 1k 11
(“WW/<01) - ZA,M
X 2; €XP(—ni€i/T) exp(_"jgj/T)ZA—n,a,'—nja,.M—n,-nr,-—njmj (2-41)

711,11}

serves to obtain second moments of the multiplicity distribution

5. .,
' " _ —j‘] I I I
("1"!) — ZA,M EC”! ZA—(I+l )a,,M—(1+I )m, (2.42)
l
+ Z {ZSIJICIZA—lak~M*1mk + ECICI’ZA—lak—l'akr,M—lmk—I’mkr }
A,” I w

Neglecting the terms with Cg where l > 1, leads back to the nondegenerate result, Eq. (2.36).

Calculating the multiplicity distribution for degenerate systems with the constraint of

total isospin conservation becomes difficult because the analog of the cycle diagram, C5,“

23

in Eq. (2.38), needs to be analyzed for isospin i and charge 11}. Such a cycle diagram with

a particles and isospin projection m can be written as

Xa,i,m(’1j) E 2(d,a,i,m|eTH/T|B,a,nj,m) (B,a,nj,m|0t,a,i,m), (2.43)
(LB

where the sums over 0t and [3 correspond to sums over all states with fixed (a,i,m) and
(a,nj,m), respectively, and ("it represents a cyclic permutation of particles, which are as-

sumed to be distinguishable.

The multiplicity distribution can be calculated in terms of these cycle diagrams

1 .
WA7],M(nj) = K 2 an-‘m(723-)WA_a‘,r,M_m(nj—n3-)(I’,M—m;r,m|I,M)2. (2.44)

a,i,m,n’j,l'

The general cycle diagram accounting for all energy levels, xa,,-,,,,(n'j), can be derived in

terms of cycle diagrams for a single-level system, 7d ”,

I
Xa,i,nr,nk : x51,i),m,nk 285 CXP (—aEf/T)a (2-45)
3

where I indicates the single-particle energy levels with energy E5. Applying, Eq. (2.44) to

a single-level system leads to a recursion relation for x“)

(l (I)
XA,I,M(nj) -_— AWA,I,M(nj) (2-46)
I l) .
— 2 x33,"(nj)W/§_0J,‘M_m(nj—nj) (I',M—m;z,m|I,M)2,

a<A,i,m,nj,l’

where the distribution Wm

a im(nj) will be derived in the next subsection.
If other particles, like resonances, are to be included in the ensemble, one can either
insert a sum over species into Eq. (2.44); or calculate W separately for each species and

convolute them to find W for the entire system,

(I) <2) . 2
WA,,,M(nj) z 2 w ,~,,‘M,(n’j)w __A,J,,’M_M,(nj—n'j)(I'M',I",M—M’|IM) .
A' l’ M' n’- l”
1 9 7 j,

(2.47)

24

In these calculations, W was calculated separately for resonances neglecting symmetriza-

tion and then convoluted with W calculated for pions with proper symmetrization.

2.3.3 Single-level Multiplicity Distribution Conserving Total Isospin

By definition, the distribution W( I)

a,i,m

(n I), needed in Eq. (2.46), can be written as the product

of the partition function and the multiplicity distribution,

W“) :01“) l" (nj), (2.48)

a,i,m,nj a,i,mpa,i,m

where the partition function for a particles in a single level with energy E is

(r) _{exp(—aE/T) ,ifa+ieven

(o (2.49)

“31'" _ O ,ifa+ i odd.
(1)

Pa i,m(”j) is the probability of observing nj pions in a single-state system that contains a

total of a pions with total isospin i and projection m. This probability distribution has to
be calculated for only one type of pions because the pion occupation numbers are related

through
m = n+ — n_, (250)
a : no+n++n_. (2.51)

In the following, probability distributions will be derived for positive pions. The isospin

wave function of the system can be written in terms of eigenstates of the number operators

(a+m)/2
laaiam> : 2 aa,i,nr,n+an>ln+>ln—>a (252)

IL}. 2m

where no : a + m — 2n+ and n- = in — m. The coefficients in Eq. (2.52) are related to the
probability distribution by

l
p(a,i)’m(n+) : (aa,i,nr,n+)2- (2.53)

The isospin wave function |a,a,a) can only be constructed if all pions in the state are

positive, i.e., n+ = a, therefore

__ 1 ,If ".1. Z a
a“’“’”‘"* _ { 0 ,otherwise. (2'54)

25

Leaving the pion number a and isospin i = a fixed, the isospin lowering operator can be

applied 1.. : x/2(1t:1to + train.) to reach lower values of m,

 

I_|a,i,m) = \/i(i+ 1) — m(m—1)|a,i,m — 1). (2.55)
The LHS of Eq. (2.55) expands to

I- zaa,a,m.n+ la ‘1’" m —' 2n+)|n+>|n+ _ "1) :

 

 

[1+
(a+m)/2
V5 2 \/0 +m — 2n+\/n+ — m + 1a....n.,..la+ (m -1>— 2n+>ln+>ln+ — (m -1)>
"+2"!
(a+m)/2——l
+ fl 2 \/a+ (m — l) — 211+ \/n+ + laa,a,m,n++1
n+2m—l

x |a+ (m — l) — 2n+)|n+)|n+ — (m — 1)). (2.56)

When the coefficients on the LHS are match with those on the RHS of Eq. (2.55) a recursion

relation is obtained,

 

 

 

2
aa,a,m—l,n+ : \/I(I + I) { \/(a + m — 2n+)(n+ — m + l)aa,a,m,n+ (257)

— m(m — l)

 

+ \/(a+m — I — 2n+)(n+ +1)aa,a,m,n++l}-

So far, only the coefficients 01mm," + for a = i were found. With the help of the isoscalar

operator

U2 2 2min: — 11:71: (2.58)

that creates two pions without altering the isospin, higher values of a can be reached,

Ugla,i,m) : Na,,'la+2,i,m), (2.59)

where the normalization constant is

 

Na) = \/(a+2)(a+3) —i(i+ 1). (2.60)

26

Matching the coefficients on both sides of Eq. (2.59) to each other leads to

l
aa+2,i,m,n+ : N 2Vn+(n+_m)aa,i,m,n+—l

 

— \/(a+m — l — 2n+)(a+m — 2 — 2N+)aa,i,m,n+}- (2-61)

In case, a + i is odd, no combination of pions yields the isospin wave function la, i,m),

as is evident from Eq. (2.49). Therefore

610mm 2 o , ifa+iodd. (2.62)

Equations (2.54), (2.57), (2.61), and (2.62) completely determine the coefficients aa.i.m,n+ ,
(l)

a i,m

which in turn define the probability distributions p (n+) and, therefore, the single-level

partition functions Wm

a,i,m,n,--
If some states 11,12, - -- ,18 have degenerate energies, the amount of computations can
be reduced by including only one energy state I with a g-fold degeneracy. The partition
function for such a degenerate state is obtained by successively convoluting the partition

functions of states [1,12, - - - ,18 with the aid of Eq. (2.27).

2.3.4 Isospin Fluctuations

Given some system with quantum states 01, pion isospin fluctuations can be defined as
G2 = 2 (01 | (N+ + N_ — 2N0)2| or), (2.63)
(I

where N+, N. and No are the number operators of the respective pions. These isospin
fluctuations could be computed through multiplicity distributions or with the expressions
for densities and higher moments that were given above. However, when total isospin
should be conserved as well, multiplicity distribution calculations are slow and expressions
for densities and higher moments are difficult to derive. Instead, the isospin fluctuations

will be derived for a system in an isosinglet in terms of isospin projection states.

27

The operator in Eq. (2.63) is a product of two rank-2 spherical tensors components

(N1 +N_ — 2N0)2 = 6720720, (2.64)
where
720— _ 2— f (11+, [11,. ,- +11 Hn_ —2n3,,.n0,,-), (2.65a)
T211 2 2 71.5 (“3,1”¥,i— nlyinoj) , (2.65b)
.
T212 = 211L111). (2.650)

i

A product of spherical tensors can be decomposed into other spherical tensor components

720720 2 Z(20;20|JM)AJM. (2.66)
J,M

By the Wigner-Eckart theorem only Aoo contributes when contracted between isosinglet

states, and because Aoo is an isoscalar one can write

(GZ>_ 672%} (A00)—_ 65(20— M2) (MleolM). (2.67)
I: ._

Some algebra leads to

A00 2 2(2M;2, ——M|00)T2MT2,_M = Ago +7180“, (2.68)
M
where
3 3 1 2
A00: §N+ + —N_ + No+ —(N+ +N_ — 2N0) (2.69)
f2 6
and

AQM_ V5.2]; (211+in: jn+ jn_ ,- +n+ (n01n+ ,no,+n01nnojn_ ,) . (2.70)

The expectation of AgoM is nonzero when particles of different charges are in the same

quantum state, or when two differently charged pions are produced into two different states

28

with a quantum correlation arising from a resonance decay. The contribution to AgoM from

the degenerate nature of the pion states can be determined via Eq. (2.41),

r r
(nk,inkl,jnk,jnk’,i> : T'AL—Mz Cj(+)[IZA— —la—k l’,—ak;M lmk— I’mki (271)
M,ll’

This contribution can be ignored in the nondegenerate limit, where occupation numbers are

small.

Contributions to ASOM from the coherent correlation between pions from resonant de-
cays can be found by expressing the resonances in terms of pion creation operators. For
example, the p+ meson, which is a member of an isotriplet, can be considered as one pion

in an 5 wave and a second pion in a p state. Referring to these two states as i and j,

1 1 t
IN) = g (Mirror 111 ,43 ) I0) (272)
How the states i and j are chosen is irrelevant since they are summed over in Eq. (2.70),

but the coherent mixture of the two permutations, which is necessary for the p+ to be a

member of an isotriplet, results in a nonzero contribution to AQM,

(AQM)— _ — 2N p0 +Np+ +Np ). (2.73)

7(

The (1) and n mesons are isosinglets and can be treated accordingly. For instance,

(if—62 81,111.31 _ jno klO) (2.74)
finj, k
which adds another term to AQM,
<A80M) : _i (N11) +Nn)- (2.75)
\/§

2.4 Numerical Results and Comparisons

In heavy-ion reactions, and perhaps in pp reactions, pions reinteract with other pions in

their neighborhood, or domain, and might be expected to sample a large portion of the

29

available phase space. As isospin should be conserved in each domain, it seems reasonable
to explore distributions for a few dozen pions rather than creating an ensemble of a few
thousand pions, which could be treated as a grand canonical ensemble [71]. One of the
study’s goals is to understand how many pions are required for conservation constraints to

become irrelevant.

In the following, symmetrization and resonances are first ignored in order to focus on
the effects of conserving isospin, subsequently, the effects of symmetrization and reso-

nances are illustrated by considering a simple example.

2.4.1 Total Isospin Conservation

When quantum degeneracy and resonances are ignored, isospin distributions are indepen-
dent of energy levels or temperature. Therefore, the results presented in the following
are generic to any system where only pions are considered and the phase space occupation
numbers are small. A random distribution ignoring isospin conservation, i.e. a mixed-event
construction, will serve as a benchmark.

1 N N!
Prandom("0) : _ ——'—'-—'
3 n++n_+n0:Nn+.n_.no.

(2.76)
Unlike distributions that conserve isospin, this distribution allows both even and odd num-
bers of neutral pions and is, therefore, scaled by a factor of two to compare the width with
that of the other distributions. Secondly, when pion creation is constrained to isoscalar
pairs, as in [128], the distribution can be considered as a binomial distribution of pairs

where one third of the time the pair is comprised of two neutral pions and two thirds of the

time the pair is comprised of a positive and negative pion.

. . — 1 ”0/2 2 (IV—"0V2 (IV/2)!
ammo—(3) (3) (no/2>!(N/2—no/2)! (2'77)

This pairwise distribution is broader than the random distribution by a factor of \/2, as can

 

be seen in Figure 2.2.

30

 

 

1 r 1 T l l W
I // \ fl \
2 I; \
10 T\
/\O \
C '5 il
v
0' 10“ \\
\
I»
10-6 r J r r l J L r J l 1

Figure 2.2: The probability of observing no neutral pions is shown for a system of 12 pions.
Symmetrization and resonances are neglected. A random distribution, multiplied by 2 for
comparison, is represented by a solid line and is used as a benchmark. A distribution result-
ing from isoscalar pairs of pions (squares) is significantly broader, whereas a distribution
including all isosinglets (circles) approaches the random distribution.

Finally, the isospin distribution for nondegenerate particles is calculated with the meth-
ods of Subsection 2.3.1 with all 12-pion isosinglet states being considered. The constraint

of exact isospin conservation only modestly broadens the distribution relative to the random

distribution as shown in Figure 2.2.

These findings are underscored by comparing the isospin fluctuations as a function of
total pion number, as shown in Figure 2.3. When pion emission is constrained to isoscalar
pairs, fluctuations are twice as large as compared to the random case for all system sizes.
When considering all N-pion isosinglets fluctuations are larger by a factor which falls from

two to unity as N approaches infinity.

31

 

j I I I l I I I I I I I T I
2.0 —Q—a-Baaaeaaaana-a+n
_ \ _
.— \ _
is - \ -
NE ' 0 fi
0 1.5 :' \\ .:

\ \

N - _
o __ g _

 

1.0 -

 

\

- Geeeeeceo££%
I I I I I I I I I l I

O 10 20 30

A

Figure 2.3: The fluctuation G2 for nondegenerate systems as a function of system size for
an ensemble of isoscalar pairs (squares) and one including all isosinglet states (circles). The
fluctuations have been normalized by the fluctuation for a random system. The fluctuations
for the pairwise distribution are twice the fluctuations for a random distribution, whereas
including all isosinglet states relaxes the constraint and results in the same width as the
random distributions for large systems.

2.4.2 Including Symmetrization and Resonances

To illustrate the effects of Bose-Einstein symmetrization, total isospin conservation, and
resonance decays, an assumption must be made about the available single-particle energy
levels that are summed over in the single-particle partition function. Assuming the model

system is confined to a cube of volume V, the energy states are obtained

 

112
em) = \/F (n2 + m2 +12) + M, (2.78)

where R = V1/3, M is the mass of the particle, and n,m,l are chosen to be half-integers.
The choice of half integers, instead of the more usual integers, de-emphasizes zero-point
surface energy effects and seems more physical if the confinement to the volume does not

arise from an infinite potential well. This becomes important when systems are confined to

32

a small volume.

Isospin distributions were calculated for a system of 24 pions at two densities, 0.3 fm’3,
which is well above breakup densities for hadronic collisions, and 0.1 fm‘3. The tempera-
ture was chosen to be 125 MeV. The Bose-Einstein nature of pions was taken into account
when calculating multiplicities according to the formalism from Subsection 2.3.2. The mul-
tiplicity distributions were then convoluted with those of three resonant states, the isotriplet
p mesons and the isosinglet (u and n mesons, which were treated as nondegenerate systems

according to Subsection 2.3. 1.

As expected, the isospin distributions for symmetrized pions in an isosinglet are broader
than the random distribution when resonances are neglected, as shown in Figure 2.4. This
broadening is especially pronounced at high density. However, the inclusion of resonances
more than compensates for the symmetrization effects and results in distributions that are
narrower than the random distribution. Figure 2.5 displays fluctuations as a function of den-
sity for the 24-pion system. The dramatic broadening induced by symmetrization at high

density is counteracted by a remarkable narrowing when resonant states are considered.

The conditions for pions to prefer forming a resonance are similar to the conditions
for symmetrization to be important, i.e., a high phase space density. At low temperature,
where the mass penalty for resonant formation would play a larger role, resonant effects
would become relatively less important than symmetrization. However, for the temperature
of 125 MeV considered here, the resonant effects overwhelm the effects of symmetrization

at all densities.

The presented calculations with only three included meson resonance are largely schematic,
a more realistic calculation requires more resonances like strange mesons and baryonic res-
onances, which need to be added along with strangeness and baryon number conservation.

The formalism scales linearly with the number of particle species if the resonances were to

33

 

 

I I I I I I I I r I
-1 c : i‘
10 z :
”/V/ \\
1. \ 0
—3
10 #l n\\
r / Q
\\ \
_5 I: “ \
10 _ -3 \ 0
/\o p—O.1 fm \ \
C I I L I J I l I L L \ I
v I I kI‘ I I I T I I I
D_ _1 O ,__ = g
1O 121/ :3 is
\ \
/ Q

10 95 \151 \q

- I \
10 == p=o.3 fm-3 I
1

Figure 2.4: The probability for producing no neutral pions in a system of A = 24 pions
at T = 125 MeV is shown for ensembles restricted to an overall isosinglet and with sym-
metrization included. Calculations with resonances (squares) and without (circles) are dis-
played. Although symmetrization broadens the distribution relative to the random distribu-
tion (solid line), the inclusion of resonances results in a narrower distribution as compared
to the random distribution.

be included without additional conserved charges. Theoretically, any amount of quantum
numbers can be conserved as long as they commute with the isospin operator. Practically,
every new index, which has to be added to the partition function and multiplicity distri-
butions, to conserve another charge increments the number of loops by one. The increase
in run-time has little consequence for calculations of partition functions and direct calcu-
lations of isospin fluctuations, which are virtually instantaneous, whereas any additional

index would significantly slow down calculations of multiplicity distributions, which take

34

 

mixed
N
l
a

N ..
C) - 9/0’ -+
N\ 1 ,1 o/O/ _j
l— \B’ —
F El B‘G‘E‘B-EPEJ

I
l

 

 

 

0.0 0.1 0.2 0.3 0.4 0.5
density (fm’s)

Figure 2.5: Isospin fluctuations, or squared width of the isospin distributions, as a func-
tion of density scaled by the width of a random distribution. Calculations with resonances
(squares) and without (circles) are displayed. Calculations were performed for 24 pions
restricted to isosinglet states at a temperature of 125 MeV. At high density, the distribu-
tions are broadened by including symmetrization and narrowed by including resonances
produces.

on the order of one minute.

The aforementioned caveats are not expected to become major obstacles in the appli-
cation of the presented formalism. Although the particle multiplicities are high, in the
thousands, in possible physical applications like relativistic heavy ion collisions, the local
nature of charge conservation would limit the number of particles considered at any given
time. The system under consideration would have to be broken into domains, in which
the charges are conserved locally, then calculations proceed one domain at a time. Each

domain would have a relatively small number of particles with which to cope.

Also, the limitations of any comparison with experiment should be considered before

tackling the more numerically challenging problem of including strangeness and baryon

35

number. It is difficult to detect neutral pions since each neutral pion decays into two pho-
tons. Furthermore, the ability to identify neutrons and kaons, especially those kaons which
then decay into pions, should be considered. Given the inherent complexity of any such
measurement, it seems proper to stop short of performing more complicated calculations
without a commensurate consideration of the measurement. Nonetheless, several valuable

lessons can be gained from the schematic calculations presented here.

Finally, it should be noted that the isospin fluctuation G2 was calculated both from the
distributions themselves and from the methods described in Subsection 2.3.4. Although the
two sets of numerical calculations have little in common aside from the functions used to
generate the single-particle levels, the two sets of moments agreed with each other within

the numerical accuracy of the computer.

36

Chapter 3

The Microscopic Evolution Model
Gromit

The transient quark-gluon plasma (QGP) is not accessible by direct experiments due to its
short lifetime, which is on the order of 10’23 s, and its small size, which is on the order
of 10‘14 m. However, direct measurements may not be necessary to answer the most im-
portant question at hand. Did the QGP actually exist during some stage of the collision?
The existence of the QGP must be established from the information contained in the col-
lision debris recorded by the detectors. An observable that would indicate the presence of
a QGP during the collision is usually referred to as a signal of the QGP. Theoretical work
on this topic has been marred by the complexity of the underlying theory of QCD and by
the large number of produced particles. Typically, around 6000 particles are created in a

gold-on-gold collision at the Relativistic Heavy Ion Collider (RHIC).

After research on analytical theories failed to produce firm theoretical predictions for
a signal of the QGP, the focus shifted toward numerical modeling of relativistic heavy-ion
collisions. The application of such models has proven to be indispensable, particularly for
the design of experiments, the analysis of data, and the extraction of detailed predictions.
The nature of these models can be statistical, phenomenological, semi-classical, classical,

or inspired by field theories. Even a classical description of relativistic heavy-ion collisions

37

can be justified, given that the thermal wavelength of the system is small compared to the
size of the system. Unfortunately, the proposed QGP signals are generally model depen—
dent in the sense can that an observable proposed by one model to be signal of the QGP
can often be explained through other mechanisms by one or more other models. These
discrepancies stem from the incompatibility of basic assumptions in the different models,
which make it difficult, if not impossible, to compare one model to another. Hence, in ad-
dition to the search for model-independent signals much effort is spent on identifying the
model components that are responsible for the diverging predictions on observables. The
different components of currently available models are difficult to study in isolation, since
the implementations of the models are generally supplied “as is” with little room for ad-
justments. In addition, questionable components often cannot be switched off to examine

their influence.

In the following, a new modular framework for microscopic transport models based on
the BTE, Gromit, is introduced after a brief review of its theoretical foundation. Gromit
is developed and written within the RHIC Transport Theory Collaboration (RTTC) in the
hope of moving beyond the limitations of current microscopic models. It provides a gen-
eral framework that could be used for any degrees of freedom and any interaction; it is ex-
pandable, modular, transparent, and easy to debug and maintain. The reason for choosing
microscopic transport models was their capability of describing nonequilibrium dynamics,
which hydrodynamical models fail to describe due to their assumption of local statistical

equilibrium [85, 109].

3.1 Transport Equations

Transport theory has a long history dating back to the works of Boltzmann and Maxwell on

the kinetic theory of gases in the 1870’s. Over time, the original Boltzmann equation has

38

been modified to incorporate new physical insights like relativity, quantum statistics, and

the nuclear mean field.

The classical relativistic Boltzmann transport equation (BTE) describes the evolution
of the single-particle distribution f (x, p) in phase space with spacetime coordinate x and

momentum coordinate p.
ppaufctrp) : C(xap)v (31)

where C (x, p) is the Lorentz-invariant collision term

 

d3p d3pl ld3p I I I I I
C(41)): 2/p° m p")! p,olffiW (paprlp,pr)-ffrW(p,pilp.120]. (3.2)

with f1 :2 f(x,p1), f’ = f(x, p’), and fi = f(x,p'l). The transition probability W is defined

in term of the scattering matrix T as

W(p,pi IpCp’r) = 8”’(p + m — p’ - p'1)|<p',p’r lTlp,pr>|2- (3.3)

In Eq. (3.2) Boltzmann’s Stoflzahlansatz, or “assumption of molecular chaos” is used to
derive that the number of binary collisions is proportional to the products of phase space
densities. The collision term also assumes that interaction are local since all phase space

densities in Eq. (3.2) are evaluated at the same spacetime point x.

Uehling and Uhlenbeck modified the transition probability W in the B113 to account for

identical particle statistics [149, 148].

 

3p d3pl 3! ~
p”6uf(xp)= Z/dpfi ‘ ”"1 [mu:tn)(1:tf>W(p',palp,p1)—

p10 p10
ff1(1if')(1ifi)W(P MW, 191)] (3-4)

where the negative sign is used for fermions and the positive sign for bosons. The modifica-
tion in the Boltzmann-Uehling-Uhlenbeck (BUU) equation becomes important only at high
phase space densities, which are found for fermions in low-energy heavy ion collisions due

to the limited available phase space.

39

The BTE was further modified by Vlasov, who added a force term to the nonrelativistic

form of Eq. (3.1) to account for the nuclear interactions.

3 l
(5+V'Vr—;(VrU)'Vp) f(x,p)=C(X,p), (3'5)

where U an external potential, which is often approximated by a mean field. In the Vlasov
equation, interactions between particles are considered to arise from two distinct sources,

from the external potential (or mean field) and from cross sections.

3.1.1 Monte Carlo Solution

The BTE, Eq. (3.1), due to its nonlinear nature, cannot be solved exactly by either analytical
or numerical methods, except for grossly over-simplified situations. The common approach
is, therefore, to get an approximate solution of Eq. (3.1) from a Monte Carlo simulation.
Often, the quantity of primary interest is not the phase space density, f (x, p), itself but an
observable that is a function of the phase space density 0( f (x, p)). The following steps
are carried out to determine any observable, including f (x, p) itself. First, the phase space
density is sampled at randomly chosen points, such that the density of sampling points
corresponds to the phase space density. Then, the observable is calculated from this sam-
pled phase space density and, finally, averaged over a large number of samples to obtain an

approximation for the observable.

This procedure has rather strong analogies to the experimental determination of many
observables. Each sample in the Monte Carlo scheme is can be interpreted as one physical
event. The sampling points correspond to the physical particles. In experiments, like in
the Monte Carlo simulation, the observable is often calculated from each event and then
averaged over all events. For convenience, the analogous physical terms, particles, events,

and observable, will be used, when referring to the Monte Carlo technique.

However, not all observables can be determined in the fashion mentioned above. For

40

example, multi-particle correlations, are destroyed by the assumption of molecular chaos,
which introduced the one-particle phase space density into the BTE, and cannot be deter-
mined from Boltzmann-type models. Another class of observables, fluctuations, are lost in
the Monte Carlo approach, as observables are averaged and smoothed out. However, the

technique can be adjusted to simulate fluctuations [31].

The time evolution of the phase space density can be modeled by moving the particles
on trajectories consistent with the mean field, if one is present, between points of inter-
action. Collisions are determined to sample the transition probability from the transport
equation by prescribing that two particles interact if they get within a certain distance of

each other. This interaction range is determined from the cross section

0'
r = g. (3.6)

In an elastic interaction, particles merely change their momenta, whereas in inelastic colli-

sions particles from the incoming channel are replaced by other particles.

3.1.2 Scaling Properties of the Boltzmann Equation

The momentum change described in the previous subsection is generally executed instan-
taneously. In most cases information about a particle’s existence is transferred to another

with infinite speed as the particles are likely to be separated by a finite distance

d g \E (3.7)

This breaks special relativity and leads to a wide range of complications, fundamental
problems like the violation of causality [103] as well as practical problems with collision
prescriptions like those that will be discussed in Subsection 3.2.5. However, the impact of
these violations can be minimized by taking advantage of a scaling property of the Boltz-

mann Equation, which remains invariant under simultaneous transformation of transition

41

probability and phase space density [157, 32, 161]

w —. w’ = Y, (3.8)
A
f —> f’ = A f. (3.9)

The scaling factor, or sampling factor, 2» could be any positive real number. However the
sampling factor is usually constrained to be integer so that the the scaling of the phase space
density can then be implemented by multiplying the number of particles by it. The most
common way to implement the scaling of the transition probability is through reducing the
cross section by a factor A. This transformation leaves the mean free path, I = 1/no, and
the scattering rate constant, thus preserving the thermodynamic properties of the system.
The interaction range from Eq. (3.6) vanishes in the Boltzmann limit, A ——> 00, leading to
entirely local interactions. This restores Lorentz covariance and causality and removes the
problems that come along with their violation. Examples of these effects will be discussed

in Chapters 4.

Yet, there is a potential caveat to the preceding approach. Although the covariance
problems can be solved by choosing a large sampling factor A, it must be studied whether
ignoring the finite interaction range has affected the results. More realistically, particles
do interact at a finite range through mutually coupling to classical fields or through the
quantum exchange of particles. A variety of ideas have been discussed in the literature for
including such interactions in ways which preserve covariance and are consistent with the
uncertainty principle [39]. However, due to an assortment of technical challenges, all such
approaches remain in the development stage. The effect of finite-range interactions can be
incorporated into a Boltzmann treatment with the following prescription. The probability
of scattering, not the cross section, is reduced by a factor of A, i.e. not all particles passing
within a distance m of one another would scatter. The effects of finite-range interaction

will be investigated in detail in Chapter 5.

42

The use of a large sampling factor has an important advantage for the Monte Carlo
solution. As mentioned previously in (3.1.1), the sampling density in the Monte Carlo
calculation corresponds to the phase space density. Increasing the phase space density by
a factor 2» implies that the number of sampling points in the Monte Carlo method must be
multiplied by a factor A. This allows a denser and a more accurate sampling of the phase

space density f (x, p), hence, the Monte Carlo result is exact in the Boltzmann limit.

Nevertheless, an exact solution is not generally required for practical purposes, es-
pecially given the uncertainties in the experimental data to which model predictions are
compared. As a consequence a finite value for the sampling factor A can be sufficient for
extracting observables from Boltzmann-type transport models. The challenge, however,
remains to show that computed observables indeed converge in the Boltzmann limit and to

identify a lower bound for the sampling factor that still yields satisfactory results.

3.2 The Framework

Gromit has the potential for being employed for partonic as well as hadronic transport
models. One application of Gromit for a model with partonic degrees of freedom was
demonstrated in [1 15]. This thesis, however, will be limited to the discussion of a particular
hadronic transport model and will mention only those features of the framework that pertain

to hadronic degrees of freedom and their interactions.

In addition to its flexibility in the selection of degrees of freedom and interactions be-
tween those, Gromit offers the possibility of running simulations in two different coordi-
nate systems. Cartesian coordinates are generally sufficient for simulations at lower-energy
heavy ion collisions but their use becomes inefficient and numerically unstable at ultra-
relativistic energies such as RHIC energies. The reasons for this will be pointed out in

3.2.5.4. When using a Bjorken coordinate system, the time t and the coordinate along the

43

beam axis, z, are replaced by

I = VIZ — 22, and (3.10)

l t+z

=—1
n 20gt_z,

(3.11)
respectively. The coordinate r] coincides with the rapidity of an observer who has moved
from the origin at t = 0 to the spacetime point (1,2) at constant velocity. This observer will
be referred to as “co-moving”. The transverse coordinates, x and y, remain the same as in
the Cartesian coordinate system. Bjorken coordinates are very useful at ultra-relativistic

energies, particularly so if the system is boost-invariant along the beam axis. Details will

be presented in 3.2.5.4.

3.2.1 Hadrons

By default, Gromit contains all hadrons and their properties tabulated by the Particle Data
Group in the Review of Particle Physics [81]. However, Gromit is a general framework
for arbitrary degrees of freedom and is not limited to these particles or their properties.
It may be desirable to define hypothetical particles with simple properties, when studying
some schematic system; to deliberately alter certain particle properties; or to select only a
subset of known particles, when studying the influence of particular particles or properties.
Gromit facilitates these studies through a simple initialization file, in which the degrees of
freedom and their specific properties can be specified. An example of such a definition
with arbitrary particles is shown in Appendix A.l. To select a subset of known particles it

is sufficient to specify the names of multiplets as illustrated in Appendix A.2.

3.2.2 Initial Particle Distribution

Gromit offers the choice of several initial conditions for the transport model in order to

provide the greatest flexibility . In addition, the framework permits the introduction of user-

44

provided initial conditions as new modules. Among the currently available schemes for the
initial condition are the Heavy Ion Jet INteraction Generator (Hij ing), reading a file with

a list of initial particles, and boost-invariant initial distributions with simple geometries.

Boost-invariant initial distributions, fl, for a certain particle type i are usually gener-
ated by sampling their momenta according to the Jiittner distribution [96] for a canonical
ensemble

._ V
W) , (3.12)

l
f" 2 carrier 6"" < T
where pv is the momentum of the particle, nv is the velocity of the reference frame, in
which the thermodynamic temperature T is measured, and p is the chemical potential. The
particle’s spin projection is sampled uniformly over all possible values. In case the number
of particles N,- to be produced is unknown, they can be determined from a grand-canonical

ensemble.
N,- = fdx3 dp3 gif; = vififirefl/T [am m2K2(m/T)A,-(m), (3.13)

where V is the system volume, m and g,- are, respectively, the mass and the spin degeneracy
of the particle, and K2 is the second-order modified Bessel function of the second kind.
The spectral function, Ai(m), is equal to 5(m — mg) for stable particles and Eq. (3.14) for
resonances. There are several possible choices for the distribution in coordinate space,

including: uniform over an ellipse, Gaussian, and a binary collision profile.

The grand-canonical ensemble mentioned above does not conserve quantum number in
each event, although it does conserve quantum numbers when averaged over many events.
An initial condition based on the algorithm in 2.2.3 conserves quantum numbers as well as
total isospin and can be easily added to the framework. Such an initial condition would be
especially useful for studies of fluctuations and balance functions, when generating parti-

cles either at the interface to a hydrodynamical model or in a boost—invariant distribution.

45

3.2.3 Resonances

When a resonance is created, its mass is Monte Carlo sampled according to a Breit—Wigner

distribution

_ l F(m) .
A(m) — 1V (mR —m)2+I‘(m)2/4e(m —mmm), (3.14)

 

where N is the normalization constant, which in the case of mass independent width equals
2n, mR is the pole mass of the resonance according to [81], and mm," is the threshold mass

for the resonance production.

Gromit contains a loadable table of all known decay branches for resonances from the
Review of Particle Properties [81] as long as the final state contains two or three particles.
However, decays will be generated only for those particles that were explicitly included in
the simulation. Currently, decays into three particles violates detailed balance (see 3.2.4.6)
in this framework as the inverse process 3 —> 1 has not yet been included. Gromit offers two
options for the partial decay width, but it was ensured that additional functions for the decay
width could be easily incorporated into the model. The simplest choice for the partial width
is a constant, where kinematically allowed, and zero otherwise. The other option for the
partial width implements a parameterization for a mass and momentum dependent width

suggested by [99, 156].

1+1
(1) _ fl “5) 2'“ (W)
Y (Samnmlamil) — \/3 (k(n12)) k2(S) +62 3 (315)

r
where the parameters 1 and [3 denote angular momentum and a cut-off value, respectively.

 

The width is set to zero, where kinematically prohibited.
The width has to be related to the lifetime of a resonance 1: in a manner consistent with

the ergodic theorem [54]. The lifetime for a usual Breit-Wigner-type resonance is

r = (3.16)

1
F

46

At the end of its lifetime, the resonance is decayed into the final-state particles. A wide-
spread practice, when modeling resonance decays in microscopic models, is to place the
decay products at the spacetime point where the resonance was located. This is also adopted
in Gromit. This procedure, however, moves particles closer together, on average over the
entire evolution of the system, thus contributing a negative pressure. An alternative to
this practice could be to remember the spacetime separation of the particles forming the
resonance and to assign this separation to the decay products when the resonance decays.
The effects of the different decay product placement on the pressure should be studied in

the future.

The assignment of outgoing momenta is straight forward, when there are only two
particles in the final state. In the rest frame of the resonance, each daughter particle is
assigned a momentum according to Eq. (3.18), such that the total momentum in the center
of mass frame of the resonance vanishes and the direction of the relative momentum is
distributed to match s-wave scattering. Other implementations for the angular distribution
of the relative momentum are readily accepted by Gromit. The case of three particles in the

decay final state will be treated next.
3.2.3.1 3-Body Final State

The case of three particles in the decay final state is more complicated. According to [81]

the 3-body phase space is

1 1 , .,
r: WW |fM|2|p1| |p3|dm|2dQldS23, (3.17)

where M is the mass of the decaying particle, (|p’{|,$2’{) is the momentum of final state
particle l in the rest frame of l and 2, and $23 is the angle of final state particle 3 in the

rest frame of the decaying particle. The momenta Ipil and |p3| can be expressed using the

general form for the particle momenta in the center-of-mass frame of two particles with

47

momenta p1 and p2 and masses m1 and m2

 

|p1| = |p2| = k(\/s,ml,mz) : fifls— (ml +mz)2] [s— (ml —m2)2], (3.18)

where s = (p1+p2)2. Then |pi| 2 k(\/§,m1,mz) and |p3| z k(M,\/s,m3). The matrix
element M is generally assumed to be independent of mu, $2}, and 523. Therefore, the
3-body phase space can be sampled by sampling mlz, 82}, and S23 according to the the
bounded weight Ipil |p3|. This can be incorporated into a keep-or-reject algorithm that
compares |p§| |p3| for a given set of sampled parameters, mlz, S2}, and £23, to an upper
bound, which is ideally an absolute maximum. As this absolute maximum is difficult to
find, an upper bound of lp’ilmax lpglmax, is used, where |p’;|max = k(M —m3,m1,m2) and

lp3|max : k(Mrml + ”12,"!3).
3.2.4 Cross Sections

In a Monte Carlo simulation of the BTE, cross sections are needed in three distinct ways.
The total cross section om. is used to decide whether two particles are going to interact,
which they do if their relative distance is less than M. Once it has been established
that the interaction in going to happen, the partial cross sections 0’,- serve to choose the
specific reaction channel, i.e., what kind of reaction will occur and what the final state
particles will be. The probability of a certain reaction channel equals (Si/0.0.. An angular
distribution, or the differential cross section, is then needed to assign the directions of the

momenta in the final state.

Several functions have been implemented in Gromit for each of these, total, partial, and
differential cross sections. Among the options presently offered for the total and partial
cross sections are, for example, a constant cross section and a tabulated cross section that
reads in a file with cross section data. Files of experimentally measured cross sections are

supplied for a variety of hadron interactions and can be used at the user’s discretion.

48

3.2.4.1 Total Cross Section

Unless otherwise determined, the total cross section is obtained by addition of all partial

cross sections,
om. 2 20,-. (3.19)
1'

When some total cross section is very well known, it might be desirable, instead, to use tab-
ulated values, a parameterization, or a combination of these two. Gromit supplies loadable
files with experimentally measured total cross sections for a variety of hadron interactions

as well as a customizable cross section function for any user-specified parameterization.

In case a specific total cross section is requested, discrepancies can arise between that
cross section and the sum of partial cross sections. The default choice in Gromit is to ignore
this discrepancy as the two values are used for different purposes. The specified total cross
is only used for the determination of an interaction, while the sum of partial cross sections

is only used for the probability of finding a certain reaction channel.

A different scheme can be used to remove the discrepancy entirely: If A0 = 0,0, ~
20,- > 0, a default inelastic channel with partial cross section A0 is added. On the other
hand, if A0 < 0, the uncertain partial cross sections are rescaled to match the total cross
section. Most partial cross section earn the designation of “uncertain” since experimental
data are available for only a few reactions and theoretical predictions lack sophistication.
The drawback of this rescaling of partial cross sections is its inconsistency with detailed
balance (which is described in 3.2.4.6). For instance, in reactions involving well-measured
cross sections for the forward reaction but unknown cross sections for the reverse process,
e.g. np —> X Y, rescaling would take place only for the forward, but not the the reverse,

reaction.

49

3.2.4.2 Additive Quark Model

Nucleon-nucleon experiments in the 1970’s and 1980’s demonstrated that, except for the
effects of the usual conservation laws, particle production in high energy jets is largely
independent of the incoming hadron [116, 57, 20]. This concept is usually referred to as
hadron universality and it suggests that interactions between hadrons are mainly governed
by similar mechanisms. Hadron universality provides, therefore, the basis for the addi-
tive quark model that parameterizes unknown hadron-hadron cross sections by the quark
content of the hadrons [80]. In this model the total and elastic cross section are given,

respectively, by

 

 

2 "11+th S] 82
:4 —- 1- .4 — .4 .2
no 1 . 3-...)(1 . 3-..) <3 .
and
0., = 0.0390342, (3.21)

where the cross sections are given in units of 10‘3 barns, s,- is the number of strange quarks
in particle i, and m,- : 1, if particle i is a meson, m,- = 0 otherwise. This approximation of

constant cross sections is good only for very high energies.
3.2.4.3 Resonance Formation

In this model, a resonant scattering process, A] +A2 —R) 81 + 82, is modeled as two distinct
processes, a resonance formation process, A; +A2 ——> R, followed by a resonance decay,
R —+ 31 + 82. The resonance formation is assigned a Breit-Wigner cross section, that is
consistent with the measured resonant elastic and inelastic cross sections, when folded with
the decay branching ratio. The spin-averaged Breit-Wigner cross section for a collision of

two particles with spin S 1 and $2 to produce a spin-J resonance is

1t rR—MBFR 21+ 1
k(t/Emzmmaflz (MR — t/EV + 1“}./4(251+1)(252+1)’

 

CAB—4R“) : [ (3-22)

50

where MR is the pole-mass of the resonance, mA and mg are the respective masses of the
incoming particle, k(\/§,mA,mB) is the momentum of one particle in the center of mass
frame given in Eq. (3.18), FR is the total width, and rte—ma is the partial width of the

resonance decay.

The cross section given in Eq. (3.22) might diverge at the threshold, as the momentum
k in the denominator vanishes, unless the mass-dependent width in the numerator goes
to zero faster than the momentum. In case a constant width is used, the divergence in
the cross section has to be removed. One possible scheme is to introduce a cut-off in
the momentum kw. = 10 MeV, below which the cross section vanishes. This scheme is
presently implemented in Gromit but may be changed in the future. Another alternative is
to replace the mass dependent momentum in the pair rest frame k( \/§, mA , mg) by a constant

momentum k = k(MR, mmmg) to eliminate of the singularity.

Equation (3.22) is only valid for stable incoming particles because of the mass depen-
dence in the momentum k(\/s, mA,mB). This can be remedied by replacing the momentum

in the numerator by an averaged value, as suggested in [23],

fi—M'z’m' \/§—m]
(k(fi» = /dm1A|(m]) [ding A2(m1) k(\/s,m|,m2). (3.23)
Mrlmn Mani"

3.2.4.4 N N -excitations

Nucleon-nucleon excitations (NN-excitations) are much more important at AGS and SP8
energies than at RHIC energies because nucleons make up less than 10% of the particles
at the higher energies. NN-excitations are included in the model to extend its range of va-
lidity, however, it is anticipated that the main focus in the research of relativistic heavy-ion

collisions is increasingly shifting to RHIC energies. The following approximate parame-

51

terization is used for the various NN —+ X Y cross sections [23]

k 1
012434 3 (253 +1)(254 + ”if; IM(51"13,M4)|2 (324)

where km and k34 are the momenta of the particle pair {1,2} and {3,4}, respectively, in
the center of mass frame given in Eq. (3.18). The matrix element M(s, m3,m4) depends on

the particle types in the outgoing channel.

For NN —+ NA the matrix element is given by

"1.211%
2 2 2 2’

 

|914|,2(,,,,__,NA = 40000 (3.25)

where mA 2 1.232 GeV, FA 2 0.12GeV. The corresponding matrix element for NN —+ AA
is
lMl12\/N—>AA : 2.8, (3.26)

and for the other NN-excitation the matrix element is

A
M ., , 2 = , 3.27
I (9 m3 m4)| (m3 —m4)2(m3+m4)12 ( )

 

where mN. = 1.44 GeV, mA. 2 1.6GeV and the values for the parameter A are given in

 

 

 

Table 3.1. '
[Excitation [ A [GeV4] ]
NN —> NN“ 6.3
NN —> NA“ 12
NN —> AN" 3.5
NN —+ AA“ 3.5

 

 

 

 

Table 3.1: Parameter A for NN-excitation matrix elements.

3.2.4.5 Baryon-Anti-Baryon Annihilation

Rapp and Shuryak [134, 135] found that in relativistic heavy-ion collisions the proton—

anti-proton annihilation with several pions in the final state, pp —> fin, where ii = 5 —- 6, is

 

1Whereas the other A parameters are determined to best fit the experimentally available data, no such data
are available for NN —+ AA‘. Its parameter is, therefore, taken to be the same as in N N —+ AN ‘

52

substantially counterbalanced by its inverse reaction. This effect is believed to be relevant
at RHIC energies, albeit less than that at SPS energies. It would, therefore, violate time
reversal invariance, if BB-annihilation were included in a microscopic simulation without

accounting for the inverse many-body processes 1m —+ BB with n 2 2.

However, including such many-body collisions explicitly into a microscopic model is
presently difficult. Besides the practical problems involved with numerical simulations, it
remains unclear how many-body collisions can be included in the theoretical framework
of Boltzmann-type models because the Boltzmann transport equation and its modifications
only describe binary collisions. One possibility of including many-body collisions implic-
itly in future models is to subtract the cross section of the many-body process from the
cross section of the inverse binary collision. Of course, this would be possible only, if the

population of the incoming channel was overpopulated relative to the outgoing channel.

3.2.4.6 Detailed Balance

Detailed balance is a consequence of time reversal invariance of the scattering matrix and
relates the cross section of a reaction to that of its inverse. This relation can be derived start-
ing from a time-dependent Hamiltonian H (I) that is essentially time-independent (given by

HQ) with a small time-dependent perturbation H’ (t)
Hm=m+mm. aw)
Fermi’s golden rule then states the transition probability per unit time as

2 , ,
W11: §1<411H1¢o>1~p<511 (3.29)

where the (11’s are eigenfunctions of the time-independent Hamiltonian Ho and the density

of states p(Ef) is given by

—. 3.30
.115, ( >

v
pw0=mmfifi

53

If the matrix element 9%; = (¢f|H’|¢0) is time reversal invariant, i.e. leilz = lelz, then
Wfi Wif

= —. (3.31)
0(Ef) NE)

The transition probability is related to the cross section by Wfi = o fivi/ V, where ofiv; is

 

the volume “probed” by the particle per unit time. Substituting the transition probabilities

in Eq. (3.31) by the cross section leads to

ofiv;p(E,-) = O'ifop(Ef). (3.32)

Finally, the cross section for the inverse process can be written in terms of the cross section
of the forward process
2
8' P-
0,, = 442-0]... (3.33)
8f Pf

Eq. (3.33) is strictly only valid for stable particles. In case there are resonances in the
incoming or outgoing channels, the corresponding momentum in Eq. (3.33) has to be re-

placed by the average momentum given in Eq. (3.23).

3.2.4.7 Differential Cross Section

At this point, Gromit provides two different expression for the differential cross section but
more functions can be added easily, if desired. The differential cross section are given in
the center of mass frame of the binary collision. The choices for the function include the

s-wave distribution
do
d (cos 0)

and the forward-peaked screened Rutherford scattering, or “gluonic”, distribution [163]

 

: cons! (3.34)

do 91: or?
-——— = —————————-. 3.35
dt 2 (t—uz)2 ( )

Here, I 2: (p1 — p3)2 is one of the usual Mandelstam variables, or, = r/ 2 / 9, and the screen-

ing mass [.1 is chosen to match the cross section, such that

91:01?
0 = ., .
21i-

 

(3.36)

54

3.2.5 Collision Handling

It has become apparent that arbitrary algorithmic choices in transport models constitute
themselves major mode] assumptions that alter the extracted results. It is, therefore, imper-
ative that such algorithmic choices are thoroughly documented and studied. Several such

studies will be discussed in Chapter 4.

3.2.5.1 Invariant Impact Parameter

The impact parameter between two particles, 1 and 2, is defined as the distance of closest

approach in the pair’s center of mass frame,
b=min(|x2(t)—xl(t)|). (3.37)

Boosting into the pair’s center of mass frame can be avoided by writing this impact pa-
rameter only in terms of Lorentz invariant scalars, i.e. in a manifestly invariant form. The
following quantities will prove useful in the calculation of this invariant impact parameter
for two particles at spacetime positions x1 and x2 with 4-momenta p1 and p2 in an arbitrary

reference frame,

 

x = x2 — x1, (3.38)
P = 121 +122, (3.39)
q = p2 — p1 — Egg-’11P. (3.40)
b=x— :étP—gééfq. (3.41)

The impact parameter vector b is defined in Eq. (3.41) such that it is perpendicular to the
relative momentum, i.e., b - q : 0. With the above definitions, the square of the invariant

impact parameter is given by

 

 

P', 2 . 2
b2:_b#b,,-_—_x2+( Pg) +(qqf) , (3.42)

55

or expressed in terms of the original momenta,

b2 : —x2 + P%(P1‘X)2 +Pi(P2 ~X)2 - 2(P1-P2)(P1'X)(P2'X)
PiPi - (P1 P2)2

 

(3.43)
This expression for the invariant impact parameter can be substantially simplified if the
particles are massless particles,

b2 Z —x2 + 2(P1 ~x)(p2-x).
(pr-p2)

 

(3.44)

3.2.5.2 Collision Times and Ordering

The points on the particles’ trajectories, A”! and 21’ , where they reach the point of closest

approach in the center of mass frame are given by

p§(p1-X)-(p1-p2)(p2-x)

 

I
X] ZXI‘I‘ p], (3.45)
PIP%’(P1'P2)2
2 . _ . .
x12 :x2_ P1(P2 36) (P1 P2)(P1 X) (3.46)

 

Pin - (P1 °P2)2
where the same notation is used as in the previous subsection. The time components co-
incide, t; = 15, in the pair’s center of mass frame. However, since the two particles are
generally separated by a finite distance, these two times are not necessarily simultaneous in
another reference frame, which must be adopted in the numerical simulation for ordering

the interactions in time.

It is impossible to honor the individual collision times for each particle as this would
lead to acausalities in the ordering of collisions. For example, consider three particles that
are going to interact pairwise and inelastically. The time in particle i’s frame, at which
particle i collides with particle j, is denoted by tij. A possible situation that could arise for
the relationship between the various collision times is: m, < tlz, I32 < t3], and m < t23.

In order to determine which particle pair should be collided first, the following constraints

56

have to be considered. Particles l and 3 should be collided before 1 and 2 as the former
inelastic collision removes particle 1, thus rendering the latter collision invalid. However,
before 1 and 3 can be collide, particle 3 has to be collided with particle 2. Finally, before
colliding with particle 3, particle 2 has to interact with particle 1 first. In this circular

deadlock, it impossible to decide which pair of particles to collide first.

Therefore, a strategy needs to be devised for assigning a unique time to a collision for
the purpose of time-ordering the pairwise interactions. In principle, any function that takes
those two individual times and computes a unique collision time for the pair is permissible.

Obvious choices include

t t
521 = #1: 2, (3.47)
67;,“ = min(t1 ,12), (3.48)
531‘“ = maX(Ii,t2). (3.49)

Because this is an arbitrary algorithmic choice, the results of the model must be indepen-
dent of the selected prescription. This will be demonstrated for several Boltzmann-type

models in Chapters 4 and 5.

When using Bjorken coordinates, there is an additional choice for obtaining an unique

collision time for collision ordering. Either, the individual t’s are combined to a collision 1:

Tacgfth : %_E, (350)
2 2
.. 7: +r
1290?: ‘2 2 (3.51)

or the t’s and 2’s are combined first and then used to calculate a collision t

2 2
‘Féi,|=\/(n:t2) _(Z1:22) . (3.52)

 

 

 

57

3.2.5.3 Collision Prescription

Although for the purpose of collision ordering the two individual times, In and t2] , have to
be combined to a unique collision time, there is still the option of propagating the particles
to a time other than the collision ordering time when the interaction is processed. For
example, particle 1 could be propagated to I12 and particle 2 to :21. This might move a
particle back in time, perhaps even to a period before they were created. The evidently
illogical situation could be avoided by simply rejecting any such collision. However, there
is another serious problem with moving the particles to their individual times 1112 and :21.
The particles do not collide at the distance of closest approach in this scheme, which means

that the actual interaction range can be greater than bmax = ‘ /o/7:.

To demonstrate this problem, the particle separation in the center of mass frame is
calculated,

rsepz x2(t§0')—x.(zf°') , (3.53)

 

where If”; are the times and xlyz the positions in the center of mass frame, at which the

particles receive a momentum kick. The particle separation rZep needs to be compared to
d2, the squared distance of closest approach from Eq. (3.42), which is always bounded by
d2 < o/n. An analytic solution is presented here for massless particles, while the case of

massive particles is treated numerically.

The first prescription to be examined assumes that massless particles are moved to their

individual collision times, In and 121 , then

-12: 2(1—v1-v2)+x-(v1—v2) 2
56p (l-Vr-Vz)

 

, (3.54)

where x = x2 — x. and v1 and V; are the particles velocities. This means that the particle
separation at the time of interaction is larger than the impact parameter. What is even worse,

there is no upper bound to the particle separation. In case, the prescription propagates both

58

particles in the collision to a common time tcol = (112 + 121 ) / 2, it is found that

b2_r2 _ [X'(V1+V2)*t(1+Vr-V2)12.

scp“ 2(1—V1-V2)

 

(3.55)

The denominator is always nonnegative as v1 -v2 < 1. Therefore, the particle separation is

O'
rsep g b 3 g. (3.56)

Numerical studies with massive particles confirm the analytical findings for massless

bounded

particles. They also show that the average-1: prescription yields no collisions where the
particle separation exceeds «(S/7:. As a result of this discussion, particles are generally

propagated to the common collision ordering time in the Gromit framework.
3.2.5.4 Collision Mesh in Cartesian and Bjorken Geometries

To find all collisions between every particle pair in a system with N particles requires
an 0(N2) algorithm. At initialization the collisions have to be determined this way and
stored. However, after a collision has taken place, particles pairs, where neither particle
was involved in the collision, do not need to be checked again for collisions. New colli-
sions could only result from the changes to the reacting particles between those and the
other particles. The scaling of this search is 0(N) and cannot be fundamentally improved.
However, the number of checks can be significantly reduced because a particle can only

collide with other particles within its range of interaction.

This interaction range cannot correspond to a specific cross section because the poten—
tial reaction partners are unknown until after an interaction range was determined and the
adjacent cells searched for collision partners. Instead, the interaction range has to be deter-
mined form a overall maximum cross section. There is, however, no analytical or numerical

way to determine this maximum cross section in a general framework, where any function

59

 

could be used for the cross section. Therefore, the maximum total cross section are deter-
mined through a Monte Carlo scheme. They are initially estimated for each particle and
then updated as larger cross sections are found. The information about the maximum cross
sections are stored across events and yield satisfactory upper bounds for the cross sections
after a few simulations. Since cross sections for the different particles vary greatly, the
collision finding algorithm can be further sped up by recording the maximum individually

for each particle type.

To optimize the search for collisions, particle are placed into cells within a mesh. The
particle-cell assignment have to be updated continuously, such that it is accurate at any
point in time. This requires additional computational overhead but is generally justified by
the speed—up in the collision finding algorithm. After an interaction has taken place, the
modified particles need to be checked for collisions only with other particles that are in its
own cell or in adjacent cells within its maximum interaction range. The optimal size of cells
sensitively depends on the density and is determined by the trade-off between increasing the
number of updates of particle—cell assignment and decreasing the number of neighboring
cells searched after a collision. Often times, yet, it may not be feasible to use the optimal
cell size since the number of computational cells is limited by the available memory. This
restriction becomes serious particular when a high sampling factor A is chosen along with

a reduction in the interaction range.

Gromit can be run with meshes in either Cartesian coordinates (t,x, y, z) or Bjorken
coordinates (1:,x, y,r]). The latter is required if one is to implement periodic boundary
conditions in rapidity. Bjorken coordinates are also recommended if the system covers a
large spread in rapidity, e.g. at RHIC or LHC energies. In these case, the Bjorken mesh
naturally expands along with the system and, hence, maintains a more even distribution of

Particles across the mesh.

60

A collision finding mesh as described above can be useful only if the particle separa-
tions are truly bounded by the interaction range. It was shown in 3.2.5.3 that the particle
separation are bounded only if certain collision prescriptions were used. Yet another com-
plication might arise from a collision prescription that allows collisions in a particle’s past.
For instance, consider a particle for which a collision time is determined in its past. At
the time when the collision is inspected the particle will have passed the point of closest
approach and might be a great distance away from it. Hence, this collision would not be
registered when using a collision finding mesh where only neighboring cells are searched
for potential collision partners. Consequently, use of a collision mesh is inconsistent with
particle propagation backward in time. Since not all collision prescriptions are consistent
with the use of a collision finding mesh, care is required when choosing the collision pre-

scription.

Another benefit of using a collision mesh is that it facilitates implementation of periodic
boundary conditions. When using Bjorken coordinates the following transformation are

required to move a particle at one end by A11 to the other end,

1’ = ’t, (3.57)
11’: n in Ar). (3.58)
m; = "1,, (3.59)
y' = y i AT], (3.60)

where An is the system’s extend in the n-direction, the sign depends on the direction of
motion, the transverse mass m, 2 ME? — p2, and y is the rapidity. This transformation is

chosen to preserve the energy in the comoving frame, m; (n’ - y’) = m, (n — y).

61

3.2.5.5 Run-time Estimates

The use of a collision mesh as described above greatly reduces the run-time required for
the simulation of one event. This facilitates faster testing of new ideas but even more im-
portantly, it improves the quality of the results obtained from the Monte Carlo simulations
in two ways. First, more samples can be computed to increase the reliability of the Monte
Carlo solution and to reduce the statistical error bars associated with the results. Second,
numerical calculations with higher values of the sampling factor A become feasible, which
are required to extract correct predictions from the simulations. The latter point has been

mentioned in Subsection 3.1.2 and will be studied in more detail in Chapter 4.

The run-time required for one event depends greatly on the details of the model imple-
mented within the Gromit framework and can range from fractions of a second to several
minutes. On a system with a Pentium 4 processor rated at 2.2 GHz, the simulation of one
gold-on-gold event at an energy of M = 200 GeV, as described in Chapter 7, requires
about 40 seconds. However, more important than the absolute run-time is the scaling of the
run-time as the parameters of the models change. For typical simulations, the run-time of

models based on the Gromit framework are found to scale as
0(anV), (3.61)

where n is the density and V the volume of the system. This scaling is approximate but
deviations are relatively small. The linear scaling with the sampling factor is only valid if
the cross sections are scaled by l/A, which allows cell sizes to be adjusted to match the
smaller interaction range. This, however, comes at the expense of using more cells in the

collision finding mesh and, accordingly, more computer memory.

As a consequence of the linear scaling with the sampling factor, most calculations with
larger sampling factors do not require more run-time if the statistical errors are to be kept

constant. Although the run-time for one event increases by a factor A, the event also in-

62

volves a factor A more particles. Therefore, in a given amount of time, the same number
of particles can be processed, keeping the statistical errors constant independent of the

sampling factor.
3.2.5.6 Secondary Scattering

Particles that have just scattered attractively, as displayed in Figure 3.1, are very likely to

attractive repulsive

Figure 3.1: Illustration of attractive vs. repulsive scattering.

scatter again within a very short time due to their proximity. If attractive, this second scat-
tering, might lead to a third, and so on, until one scattering is repulsive. In fact, if the density
is high in some region, particles within this region scatter among each other so many times
that it effectively halts the simulation. More seriously, this kind of attractive scattering has
been shown to modify the thermodynamic properties of the system by effectively reducing
the pressure [54]. The cause for this unphysical behavior of the secondary collisions are
the short time separations between collisions. This, in turn, is a consequence of the instan-
taneous collisions that were introduced in 3.1.1. Instantaneous collisions, though, are not

found in nature where finite time scales are associated with interactions.

Such interaction time scales can be mimicked in the Monte Carlo simulation by re-
quiring that particles involved in an interaction do not interact within a certain dead time
tdcad. This time would allow the particles to separate from each before colliding again. This

prescription, however, reduces the overall collision rate.

An alternative to blocking all interactions within a certain dead time, is to restrict par-

63

ticles involved in an interaction from colliding with each other until either of them has
interacted with a third particle (“weak blocking”) or until they have both interacted with
other particles (“strong blocking”). Gromit implements both these possibilities and offer

the choice between weak and strong blocking.

3.2.5.7 Conservation Laws

Conservation of energy and momentum is relevant for Boltzmann-type models on two
scales, the macroscopic and microscopic scale. On the macroscopic scale, the fundamental
conservation of energy and momentum throughout the evolution of the collision crucially
affects the thermodynamic properties of the system like temperature and pressure. Since
the goal of modeling relativistic heavy-ion collisionsis to discern bulk properties of the col-
lision, it is important for any model to conserve energy and momentum on the macroscopic

scale.

This, though, does not necessarily imply that energy and momentum have to be con-
served microscopically, i.e., in each interaction between the individual particles in the sim-
ulation. In fact, some Boltzmann-type models choose to model the average effect of col-
lisions and, therefore, only conserve energy and momentum on average. In the Gromit
framework, energy and momentum are exactly conserved in each microscopic interaction
as well as macroscopically across the entire evolution of the collision. This is numerically
ascertained for every individual interaction in the Gromit framework and has been verified

for the macroscopic evolution of the model.

Preservation of the reaction plane has previously been implicated to affect calculations
of another important observable, elliptic flow discussed in Section 4.2, in the cascade limit
(I. = 1) [97]. The Gromit framework facilitates the investigation of this possible effect
by providing the possibility of preserving the reaction plane in each binary collision. The

stricter conservation of angular momentum in each microscopic collision cannot, even for

64

nonrelativistic cases, be implemented consistently with scattering at the distance of closest

approach.

Furthermore, it should be noted that due to the finite range of interactions in micro-
scopic simulations, charges are generally not locally conserved. In an inelastic collision,
the incoming charges are removed when they reach their point of closest approach. These
charges are then replaced by the outgoing particles that have the same net charge but dif-
ferent individual charges. In addition, the outgoing charges might not emerge from the
same points, at which the incoming charges were removed. As a consequence of moving
the charges across finite distances charges are not locally conserved, although charges are

conserved overall.

65

Chapter 4

Results and Tests of Microscopic Models

The search for the quark-gluon plasma with relativistic heavy-ion collisions is mainly ham-
pered by the lack of firm theoretical predictions of QGP signals. Theoretical predictions
are frequently made for the early stages of the collision where novel phenomena are ex-
pected to take place. The predicted effects are then followed through the collision dynam-
ics with transport models to arrive at a prediction for the final state of the collision. Besides
disagreements between different models, there are considerable uncertainties about the dy-
namics contained in the models. It is, therefore, imperative to thoroughly study the contri-

butions of model assumptions, arbitrary or not, and their effects on extracted observables.

In microscopic models based on the Boltzmann equation, the most important uncon-
strained model parameter is the sampling factor I discussed in Subsection 3.1.2. It deter-
mines how densely the phase space distribution is sampled. The exact solution of the Boltz-
mann equation is obtained with Monte Carlo methods only in the Boltzmann limit, A ——> 00.
Nevertheless, calculations have to be performed with a finite sampling factor 3., which in-
troduces instantaneous interaction across finite distances. This leads to issues involving

Lorentz covariance [103, 163], which have to be resolved by arbitrary prescriptions.

A number of these arbitrary prescriptions were mentioned in this thesis. For instance,

a prescription is needed for ordering collisions in time. Each collision can be assigned a

66

point in space time. However, since many points would have space-like separations, time
ordering is frame dependent. Another aspect of arbitrariness derives from the collision pre-
scription for altering the trajectories of the two particles. Since the particles are separated
by a finite distance the notion of simultaneously changing the momenta is inherently frame
dependent. This finite distance is dependent on the scheme chosen to implement the scal-
ing of the transition probability as the sampling factor 3. changes. In one scheme, the cross
section is scaled by l/A to account for the change in the transition probability, whereas
in the other scheme, the appropriate number of collisions is rejected while leaving the in-
teraction range constant. Preservation of the reaction plane in each scattering is another
algorithmic choice that has been shown to affect elliptic flow in the cascade limit 0» = l)
[97]. Finally, even in the nonrelativistic limit, scattering at the point of closest approach

violates conservation of angular momentum.

Furthermore, there are other issues with scattering algorithms that are mentioned but
will not be addressed here. For example, particles might scatter with an impact-parameter-
dependent probability that has a more complicated form than a simple step function. Also,
scattering need not occur instantly. Delaying the particles a certain time, before they
emerge with their asymptotic momenta, can effectively alter the pressure [54]. Such time
delays must be coordinated with in-medium modifications and knowledge of the phase
shifts in order to be ergodically consistent with thermodynamic properties at the level of

the second virial coefficient [54].

Despite this multitude of arbitrary algorithmic choices, results extracted from Boltzmann-
type models should converge in the Boltzmann limit. In the following, this will be ascer-
tained for particle spectra and elliptic flow by studying their sensitivities to some of the
above mentioned prescriptions in a schematic model. Calculations of elliptic flow in an-
other simplified model with only elastic collisions are then compared to experimental data

to verify the dynamics contained in the Gromit framework. The same model is also utilized

67

to study the effect of elliptic flow in the initial condition on the elliptic flow in final state.
Finally, a more realistic model with mesonic and baryonic resonances as well as inelastic

collisions is applied to explore the regeneration of p mesons in the collision dynamics.

4.1 Particle Spectra

To illustrate the sensitivity of results to algorithmic choices, we compare results from four
numerical approximations to the Boltzmann equation which are realized in four separate
models, each arising from a different author. All the algorithms discussed here scatter
particles instantaneously at the point of closest approach, and all the models assume a
simple s—wave form to the cross section. Each algorithm is executed using the same initial
conditions, and each is executed both in the cascade limit (A = l) and with a high sampling
factor, A = 16. Due to the arbitrariness inherent to algorithms with finite-range interactions
it should not be surprising for results to vary between codes with A = 1, but it is expected

that results from the four approaches should converge in the large 7t limit.

The four models differ in their definition of “simultaneous”, in their time-ordering pre-
scriptions and in whether they preserve the scattering plane in two-particle scattering. The

four algorithms compared here represent the following choices:

1. ZPC: A parton cascade code authored by B. Zhang [159]. Particles scatter at the
point of closest approach, simultaneous in a fixed reference frame. This time is
determined by first finding the spacetime points if, and x” , where the two particles
would scatter should they alter their trajectories simultaneously in the two-particle
center-of-mass frame, as given in Eqs. (3.45) and (3.46). In a designated laboratory
frame, a time is determined by averaging ti and t5. The trajectories of both particles
are simultaneously altered at this average time. The reaction plane is preserved in

each two-body scattering, and scattering between the same partners is allowed after

68

one particle has scattered with a third particle.

. Gromit-t: A generic scattering engine developed for the RHIC Transport Theory

Collaboration (RTTC) [137]. This should be identical to ZPC.

. MPC(0.4.0): A Boltzmann description authored by D. Molnar [109]. The principal
difference with the codes above is that the reaction plane is not preserved in the two-
body scattering. Another difference is that scattering between collision partners is
not allowed until both particles have scattered with other particles. It should be noted
that MPC has switches which make it possible to reproduce the choices described in

ZPC and Gromit-t.

. Gromit-t: Another version of cascade/Boltzmann engine developed by the RTTC
collaboration with collisions ordered by t to be appropriate for ultra-relativistic col-
lisions [137]. The spacetime points at which the two particles would scatter had the
particles scattered simultaneously in the two-particle rest frame, x’l and A” , are used
to generate two proper times, 1”] and t’z. The average 1' _-_—. (1"l + 1’2) / 2 is then used for
ordering. Both trajectories are altered when their proper time equals 1’. The reaction
plane is preserved in the scattering and scattering between the same pair of particles

is allowed after one of the particles has scattered with a third particle.

The initial condition is chosen to be a boost-invariant particle distribution as described

in Subsection 3.2.2. The initial particles are confined to lnl g 2 and a circle of radius 5 fm

in the transverse spatial coordinates. The simulation involved 2400 pions and 240 protons

at a temperature of T = 180 MeV and T = 165 MeV, respectively, to roughly achieve

consistency with spectra resulting from heavy ion collisions. The initial particle densities

are 7.64 pions and 0.764 protons per fm3. To be as simple as possible, constant cross

sections are imposed independent of the species involved.

69

 

I I I I I I I I I I I I I I

 

 

 

 

 

10
i 1 pions
<1) 10 o a=0mb oo
8 A a=10mb
A 10-1 U 0:20mb
.3; o o=40mb
E
.0“
E 101
l:
\N/ o protons
\ 10
Z
N -1
'0 10
, I 1 I I J L l l I I L I
0.0 0.5 1.0 1.5
mt—m (GeV)

Figure 4.1: Transverse Mass Spectra from Gromit-t are displayed for four cross sections.

The resulting transverse mass, m, 2 W, spectra at mid-rapidity as calculated
with Gromit-t are shown in Figure 4.1 for four different cross sections, 0, 10, 20 and 40 mb.
These calculations are performed with a sampling factor A 2 32. For larger cross sections,
the pion spectra are cooler while the proton spectra are hardly altered. The difference would

be more if the pion and proton distributions were initialized with identical temperatures.

Figure 4.2 displays spectra from Gromit-t with 0 =40 mb for four different sampling
factors, A 2 1, 2, 8 and 32. The sensitivity to A is remarkable. For small A it appears that
both the pion and proton spectra become hotter. Results converge for A ——> 00 as the A 2 8
and A 2 32 results are already barely distinguishable. In all the models investigated here,
the final-state transverse energies are higher for small sampling factors. This reduction

in cooling is attributed to viscous effects arising from finite-range interactions, as will be

70

 

IIIIIITIIIIIII

 

-2
dZN/(Z‘nmtdmtdy) (GeV )

 

 

 

Figure 4.2: Spectra for pions and protons resulting from Gromit-t with 40 mb cross sections
for four sampling factors, A 2 1,2, 8,32.

discussed in Chapter 5.

In the upper panel of Figure 4.3 the resulting spectra from the four algorithms are dis-
played in the cascade limit. The slopes vary by approximately 20 MeV. Since the hadrons
are initially separated by approximately 1/3 fm, while the interaction range is 1.1 fm, it is
not surprising that the A 2 1 result is sensitive to the scattering algorithm. The same results
are displayed in the lower panel of Figure 4.3 but with a high sampling factor, A 2 16. The

various algorithms have converged.

For each algorithm the spectra fall significantly more steeply when produced with a
higher sampling factor A. This is especially true for the proton spectra, which suggests that

radial flow is stronger in the cascade limit. Since the average transverse energy for pions

71

 

 

dZN/(antdmtdy) (GeV‘Z)

 

 

Figure 4.3: Spectra for pions and protons resulting from four models run with 0' 2 40 mb:
(ZPC — squares, MPC — diamonds, Gromit-t — triangles, Gromit-t — circles). In the upper
panel calculations are performed with a sampling factor of 1, while in the lower panel each
model used a sampling factor A 2 16. Due to different scattering algorithms, the models
generated different results for A 2 1, while generating identical results for larger A.

is also larger in the cascade limit, one can infer that longitudinal cooling is significantly

suppressed through particles interacting at a finite range.

4.2 Elliptic Flow

The elliptical shape resulted in elliptic flow which is parameterized by the observable v2.

v2 = (608901-900), (4.1)

where 90 specifies the direction of the reaction plane which contains the short axis of the

ellipse. Elliptic flow has been proposed as a means for determining the equation of state of

72

the matter at early times of a relativistic heavy-ion collision [117, 53, 121, 123, 142, 160, 4]
and is, therefore, one of the most important observables in the research of relativistic heavy-

ion collisions.

Recently, the inability of traditional microscopic models to reproduce the experimen-
tal elliptic-flow results from gold on gold collisions at m 2 130 GeV [4] has stirred
much recent theoretical research activity [145, 108, 18, 158, 147]. Some authors matched
the experimentally measured elliptic flow with a Boltzmann-type model by increasing the
model cross sections beyond physically conceivable values, which effectively mimics a
hydrodynamical description [110, 111]. Other authors could produce an agreement to ex-
perimentally data by using a hydrodynamical model for the early stages of the collision and

a microscopic transport model for the later stages [102, 93, 101].

4.2.1 Sensitivity to the Sampling Factor A

Given the interest in elliptic flow, it is important to obtain an accurate calculation of elliptic
flow from microscopic simulation. Of foremost importance is the study of sensitivities of
elliptic flow to the sampling factor. For this study, the same model is adopted as in the
previous section. The initial condition differs from that used previously in that the particles
here are confined to an ellipse in the transverse direction with a major semi-axis of 5 fm
and minor semi-axis of 2.5 fm to crudely account for the initial spatial anisotropy in a
relativistic heavy-ion collision. In addition, cyclic boundary conditions in n are employed

to simulate a truly boost-invariant system.

Both schemes for incorporating the sampling factor into a Monte Carlo simulation,
which were introduced in Subsection 3.1.2, will be used when extracting the elliptic flow.
In one scheme, the increase in the sampling factor by A is countered by scaling the cross
sections by 1/ A to maintain a constant collision rate, which is proportional to the product of

density and cross section. Extracted elliptic flows are illustrated in Figure 4.4 for two cross

73

sections, 10 mb and 40 mb, and for four sampling factors. The results for 40 mb confirm

both expectations about the effects of the sampling factor. First, observables, elliptic flow

=10mb o=40mb

 

 

 

 

 

 

0.6 TIII'IIII‘IIIWIIIIIIWI'II
BBBA=1
' 30-0124 “f.
’3 0.4 —A-A-AA=16 —1,— ,.;@==§’°'-S‘T
g ooox=32 a,» 3 0 flag.
- " .. t .31 “O _.' .. v =4
38. 0 2 «”5353 GIL- 3’9 F" _
> i- (.6. £33.. . - a
_ .‘ .' [‘1'
fl - ..
00 11::[1111114é11111111111
L -1. §._
w i 22-? '
g 0 4 - -— 9,.-Q- :2}-
46 - -§ $912.. :4- {2:115- ,.0"'§ a“;
a 553;... 2. £7" “E"
v 0 2 _ 69' Briti— u: .' 'B .3
N , : .-
> _ ,3. _ 5.0 .13 3
j H.
0.0 l l I l I l 1 1 III I 8.1 I I II I I LLLI
0.0 0.5 1.0 0.0 0.5 1.0
mt-m (GeV)

Figure 4.4: Elliptic flow is shown as a function of transverse energy for both pions and pro-
tons for several sampling factors. In this case the cross sections are scaled by a factor 1 / A
to account for the over-sampling. The reduced flow resulting for calculations with small
A (larger cross sections) is attributed to viscous effects arising from the finite interaction
range.

in this case, are sensitive to the sampling factor. This effect can be seen from Figure 3.2.2,
where larger sampling factors lead to pronounced increases in elliptic flow. Second, the
observables converge for large sampling factors. In this case, the results do not change
significantly from A 2 16 to A 2 32, so that a sampling factor of A 2 32 already approaches
the Boltzmann limit. The minimal change in the elliptic flow for 10 mb as a function of the
sampling factor suggest that the source for the sensitivity to the sampling factor are linked

to the range of the interaction. The interaction range decreases as the cross sections are

74

scaled by l / A in this scheme. Small sampling factors A correspond to the largest interaction
range, which may lead to a reduction in elliptic flow by homogenizing the collision region

across larger distances.

This hypothesis is further examined by comparing the above results to the elliptic flow
extracted with the second scheme for modeling the sampling factor. In this scheme, the
cross section, and thus the interaction range, is left constant as A is varied. Instead, only
a fraction of 1 /A of all possible scatterings are executed to maintain a constant reaction
rate. The elliptic flow extracted from this simulation proves to be insensitivity to A as

demonstrated in Figure 4.5. This result shows that the finite range of the interactions is

0=1Omb U24Omb

 

 

 

 

 

 

 

0.6 I T I I l I 1 I I l I I I I I I I I I r Fl I I
DBBA=1
" o-o-ox=4

”a? 0.4 ~AA-Ax=16 —— —
g GOOA=32 1.3:“: ”'33:

o 5’ 45"{3g‘ . .gi‘fsg’ .
1110.2 -— ,9 —— ,9 -
_ ,P’ -. ,P’ .

1’ 9'
oo 9H:}HH%H?%H}HH%%1
A P T _
(f)

g 0.4 — , -P 2
o . $1
45 ‘ 3:6“. '. ‘1' i!‘

a 5"“- ' g3
\{V O 2 — g? """‘ ’,,~_ E] —
> _ f __ #6348. _
0'0 2‘] I I I I I I I I I I ‘_l LI I I I i L I JJ

0.0 0.5 1.0 0.0 0.5 1.0

mt—m (GeV)

Figure 4.5: Elliptic flow is shown as a function of transverse energy for both pions and
protons for several sampling factors, A. The cross sections is not varied as A is varied, but
instead only a fraction 1 / A of the particles are scattered. The insensitivity to A demonstrates
that the Boltzmann limit is effectively realized, even for A 2 l.

75

indeed responsible for the sensitivity of the elliptic flow to the sampling factor, when the
cross sections are scaled with l / A. The explicit incorporation of finite-range interactions
into a Boltzmann treatment and their link to viscous effects will be pursued further in

Chapter 5.

4.2.2 Comparison to Experimental Data

Gromit is able to reproduce elliptic flow measured by the STAR collaboration [4] with a
simplified model. The goal is not to provide a physically accurate description of a rela-
tivistic heavy-ion collision, but to examine whether the model assumptions in Gromit are
consistent with physical principles. The simplified model includes only gluons that interact
elastically with different energy-independent cross sections. Each gluon is converted to a
pion via parton—hadron-duality with n+, n‘, and no being equally probable. The initial
condition is obtained by a Glauber—type model [77, 101], in which the particle distribution
is taken to be boost-invariant along the beam axis with the following transverse profile.

dN(X‘L)
dnafxzi

 

= oTA(xL+b/2)TA(xl —b/2), (4.2)

where o is chosen to fit the required multiplicities and the thickness functions are defined

7;I(XL) = /dZ PA (t/Zzi-Xi) (4-3)

The nuclear density distribution pA(r) is modeled by a Woods-Saxon (or two-parameter

as

Fermi) form

___ Do
l+exp{(r-a)/b}’

where for any nucleus po 2 0.16fm“3 and a 2 (1.12/11/3 — 0.86A"/3) fm and for the gold

 

PA (r) (4-4)

nucleus A 2 197 and b 2 0.54 fm .

This theoretical initial condition uses the impact parameter, whereas the STAR collab-

oration parameterizes the centrality of the collision by the particle multiplicity in the event,

76

The transverse profile in 4.2 has to be integrated over xi to obtain a relation between mul-

 

I
1

03

05 - r
04
Q3
02
DJ

0 I 1 1 1 1 1
0 2 4 6 8 10 12 14

b [fm]

NINMIX

I I I
I I

I
I

 

 

 

Figure 4.6: Relative multiplicity distribution as a function of impact parameter according
to a Glauber-type model.

tiplicity and impact parameter, which is shown in Figure 4.6.

Monte Carlo simulations are performed with different values for the energy-independent
cross section as well as different distributions for the differential cross section. A selection
of these calculations are shown in Figure 4.7 with cross sections of 8 mb, 10 mb, and 40 mb
and either s-wave or screened Rutherford differential cross section. Despite the simplicity
of the employed model, the elliptic flow from the simulation with 8 mb, s-wave cross sec-
tion describes the experimental data reasonably well. However, the purpose of this study
is not to fit experimental data but rather to gain insight into the applicability of the Gromit

framework introduced in the previous chapter.

It should be noted that the experimental values for the elliptic flow shown in Figure 4.7
are in the process of being revised [5]. This new analysis became necessary after novel
theoretical insight by Borghini, Dinh, and Ollitrault indicated that nonflow effects con-

tribute substantially to the elliptic flow determined experimentally and have to be removed

77

 

 

 

 

0'12 soGev Au+Au, STAR -e—.
9.. Omb, gluonic --a---
""" 9". mb, s-wave -A~
0.1 r Om b, s-wave ---¢--- _
? "‘9,
0.08 - 9 t.‘ .

>“‘ 0.06 e,
:‘
0.04 5

0.02 F

 

 

I I I

L I

o 1 L I 1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
NINmax

 

Figure 4.7: Comparison of elliptic flow from a Gromit model to experimental data from
Au+Au at 130 AGeV measured by the STAR collaboration for some choices of cross
sections. The experimental results are shown without systematic errors. The calculation

were performed with a sampling factor of A 2 32.

[36, 35, 34]. The new results have not been published yet, but the improved analysis are

expected to yield slightly smaller values for the elliptic flow for peripheral collisions. The

changes would have little impact on the discussion presented here.

4.2.3 Elliptic Flow in the Initial State

So far, theoretical studies on elliptic flow generally assumed for simplicity that the mo-
mentum distribution in the initial state was isotropic, or in other words, that there was no
elliptic flow in the initial condition. This, however, need not be the case in a relativistic

heavy-ion collision and will be examined with a Glauber-type model as explained previ-

78

ously in Subsection 4.2.2.

In a Glauber-type model, the nucleus-nucleus collision is assumed to be a superposition
of many nucleon-nucleon collisions. A certain angular momentum is associated with each
such collision between a pair of nucleons. The sum of the angular momenta would vanish

if the distance vector between the two initial nucleons at positions x1 and x2,
Ar 2 x1 — x2, (4.5)

were distributed isotropically. This condition can be inspected by examining the orientation
of the distance vector, which can be quantified by the angle between the distance vector and

the x-axis, Z (Ar,ex). The average over the cosine of twice the orientation angle
(cos(2 Z(Ar,ex))), (4.6)

would vanish if the distance vector was distributed isotropically. Figure 4.8 displays that

average for the adopted model as a function of impact parameter. As expected, the distri-

0.25 I j CF I fl

 

0.2

0.15

<cos[2 ((Ar, ex)]>
p

0.05

 

 

 

 

b [rm]

Figure 4.8: Anisotropy in the distribution of the distance vector in the initial condition
according to a Glauber-type model of a relativistic heavy-ion collision.

bution of the distance vector is isotropic for central collisions, however, that distribution

79

becomes more anisotropic as the impact parameter increases. The anisotropy in the distri-
bution of the distance vector yields a net angular momentum that might lead to momentum

anisotropy in the initial state of relativistic heavy-ion collision.

In addition, the elementary nucleon-nucleon collision employed in the model for the
initial state might contain elliptic flow. This aspect has neither been addressed theoretically
nor experimentally up to this point. Yet, rather than seeking to establish the existence or
source of elliptic flow in the initial state of a relativistic heavy-ion collision, this study will
investigate whether and how the elliptic flow in the initial state affects the elliptic flow in

the final state.

 

0.12 . . . . . J . . 4
nitial v2: 0 —e—
nitlal v1: 0.3 ----a-

 

 

o
"'
a
"

 

   

0.08 ~

0.06 H.

.3. -

 

final v2

0.04 l

0.02 ~

 

 

   

                 

I I

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Nleax

L L I l J I I

-0.02

Figure 4.9: Elliptic flow in the final state as a function of centrality with and without
elliptic flow in the initial state. The initial condition and dynamics of these simulations are
described in the text in Subsection 4.2.2. The sampling factor is A 2 32, the cross section
40 mb.

The elliptic flow is calculated from a simulation similar to the one specified in Subsec-
tion 4.2.2 with and without initial elliptic flow. The calculations are carried out for a range
of values for cross sections and initial elliptic flow and show a consistent picture. For sim-

plicity only the calculation with those parameters will be presented that promise the most

80

I I ‘2.

visible effect. The results demonstrate that even with exaggerated, large values, 40 mb for
the cross sections and 30% for the initial elliptic flow, the elliptic flow in the final state is

not affected significantly, as shown in Figure 4.9.

4.3 Regeneration of p Mesons

So far in this thesis, the Gromit framework has been employed only for simple models
that included one or two particle types and elastic collisions with energy-independent cross
sections. To move further toward a more realistic model, the particles that are most com-
monly found in a relativistic heavy-ion collision are included. The following particles and
resonances are included: the meson pseudo—scalar and vector octets and singlets as well as
the baryon octet and decuplet along with their anti-particles. These particles are allowed
to scatter elastically with cross sections according to the additive quark model (AQM). For
the included resonances the known decays are incorporated, as are the corresponding in-
verse resonance formations. This model is then used to gain insight into the regeneration
of p mesons in the collision dynamics. This insight can be used in the search for the chiral

restoration phase transition.

A possible signature of chiral restoration and rebreaking are anomalous isospin fluctu-
ations, which were mentioned in Chapter 2. Another indicator for the restoration of chiral
symmetry is derived from a predicted consequence of this phase transition. Resonance
mass distributions are believed to broaden and shift toward lower invariant mass. The p
meson is a candidate for a probing this effect. It has a lifetime on the order of l fm/c
and is, therefore, sensitive to the rapidly changing environment in a relativistic heavy-ion

collision, which is though to last on the order of 10 fm/c at RHIC energies.

The p meson has a dominant decay branch into two pions. These are not able to es-

cape the collision region without further collisions, when emitted during an early stage

81

when hadron densities are high. Subsequent collisions erase the information that the pions
stored about the decayed p, which suppresses the experimental signatures of the p mesons.
Nonetheless, at later stages of the collision, when densities are significantly lower, the pi-
ons from p meson decays are able to exit the collision region without further interactions.

These pions can then be used to deduce the properties of p mesons at lower densities.

There are other p decay channels into electron and muon pairs with respective branch-
ing ratios of 4.5 x 10‘5 and 4.6 x 10‘5. Electrons and muons have a high probability
of leaving the early, high-density phase unscathed, as their cross section with hadrons is
relatively small. It is, therefore, possible to experimentally determine properties of the p
meson at high densities from decay leptons. The small numbers of leptons produced in p
mesons decays complicate their use in the experimental extraction of p signatures. This
task, however, may be more feasible as there are indications that the branching ratios of
the p meson decay into leptons are enhanced under the conditions found in the core of the
collision region. Consequently, the properties of the p meson under different conditions
like temperature and density can be inferred from analyses of the different decay channels.
Consequently, the properties of the p meson under different conditions like temperature

and density can be inferred from analyses of the different decay channels.

The short lifetime of the p meson that makes it possible for it to probe a small range
of temperature in the collision region might also pose a problem. The p mesons from
the early stages of the collision decay exponentially, leaving few surviving resonances to
probe the later stages of the collision. This depletion might, however, be compensated
for by the formation of new p mesons in the Zn —> p reaction. This process is referred to
as regeneration. The above-described model is employed to investigate this regeneration
quantitatively. The simulations will be performed with collisions included and blocked to

illustrate the effect of the p meson regeneration.

82

The parameters of the model are chosen to provide realistic conditions. In the initial
condition, the particles are distributed according to the boost-invariant J tittner distribution,
given in Eq. (3.12), with initial local temperature Tlocal 2 278 MeV at to 2 l fm/c. The ini-
tial temperature is determined such that in a grand-canonical ensemble dN/dyIFO 2 750,
the value measured at RHIC. In coordinate space the particles are uniformly distributed
over a cylinder with radius 5 fm and |n| < 1. To make maximum use of boost-invariance,
the mesh in Bjorken coordinates are employed with periodic boundary conditions in n. The

simulations are performed with a sampling factor of A 2 10.

In the absence of interactions, the initial p mesons follow the usual exponential decay

law, as shown in Figure 4.10. The decay rate tail off after 1.‘ z 5 fm/c and is insufficient to

 

 

 

 

 

 

L l r . .1

50 collusrons ‘

50 _ \ ---- no collisions A

$540 —- a
|~

Q 30 — a
Z

'0 20 - _

10 - '-

0 1 ‘“2~k_-__4 _____
O 5 1O 15 20

Tdecay [fm/c]

Figure 4.10: Distribution of 1’s at which the p’s decayed with and without collisions.

sample the entire evolution of the collision. On the other hand, inclusion of the 21: ——> p
reaction among other interactions leads to creation of p’s throughout the evolution of the
collision. The decay of these collision-produced p mesons take place at a later time and
alter the exponential decay rate, demonstrated in Figure 4.10. The decay rate remains

sufficiently high even at times well past 10 fm/c and supplies sufficient probes for the latest

83

stages of the collision. This shows conclusively that the p meson can be used as a probe for
all stages of the collision. Experimental work is under way to extract in-medium resonance

properties.

The model can also be used to gain insight about how the regeneration of p mesons
might help in the search for the QGP. If no QGP was formed in the early, hot and dense
phase of the relativistic heavy-ion collision, the collision region would mainly consist of
hadrons. In this case, more p mesons would be regenerated than in a cocktail of overlap-
ping nucleon-nucleon collisions due to the larger number of collisions between hadrons.
As mentioned before, the number and properties of p mesons in the early stage of the colli-
sion can be experimentally determined from the leptonic decay branch of the p meson. To
investigate whether indeed more p mesons are regenerated through hadronic rescattering in
a relativistic heavy-ion collision, the same model as above is used to calculate the mass dis—
tribution of all p mesons that were created throughout the event. The calculation without

collisions shown in Figure 4.11 represents the result expected for a cocktail of nucleon-

 

 

 

 

 

2000 l I l l .1 l
--\Ndh coL
1500 _ ---no coL _
,8,
\1000 — -
2:
“O
500 - 2
O ——k”i 1 \4__ 1
0.2 0.4 0.6 0.8 1.0 1.2 1.4

rn [GeV]

Figure 4.11: Distribution of p masses integrated over the entire collision history with and
without collisions.

84

nucleon collisions. The distribution is entirely determined by the Breit-Wigner form of the
spectral function for the resonance mass given in Eq. (3.14). As expected, the mass distri-
bution receives significant contributions from 21! ——> p formation processes, when collisions
are included. In the Gromit framework, the resonance formation cross sections explicitly
contain a Breit-Wigner distribution to reproduce the resonance production correctly, which
is confirmed by the mass distribution displayed in Figure 4.11. The resonance formation
cross sections, however, do not account for in-medium modifications to the resonances and,

therefore, lack the shifts and broadening effects referred to earlier in this discussion.

Evidently, inclusion of collisions in the system affects other particles as well. For in-
stance, light particles in the final state are colder when there are collisions, as can be in-

ferred from the pion spectra shown in Figure 4.12. The temperature can be extracted from

 

   
     
 

 

 

 

I I I I I I r I I
e —— 1r col 1
103 :~ ——-- 1r no col _:
E ---~- 9 p" col 3
I —--pp' no col;
/'\ " ‘
>\ 2
'0w 10 5‘ :-
‘0“ I \ I
E _\.. _
t: 1 “\
N 10 E- \i‘.
v : \Q.
\ I “\N— ..
NZ _ “-..\,'__ .
_ '..\ -
‘0 o .k§
10 E‘ '..\. T1
: "\;
: I
I I I I I I I I I
0.0 0.5 1.0
mt—m (GeV)

Figure 4.12: Transverse mass distribution for pions and combined proton and anti-protons
with and without collisions.

85

the slopes of the particle spectra as they are proportional to the inverse temperature. Figure
4.12 also indicates that the temperatures of heavier particles like protons are not affected
as much by the inclusion of binary collisions. This is a result of the increase in the average
transverse energy of protons relative to that of the pions due to the transverse expansion.
This increase compensates for the decrease due to the work performed on the longitudinal

expansion, which affects all particles.

86

Chapter 5

Viscosities from Finite-Range
Interactions

A principal goal of relativistic heavy ion collisions is to experimentally discern bulk prop-
erties of the excited vacuum. To accomplish this aim it is imperative that the implications
of the finite size and lifetime of the global reaction are understood. Two microscopic length
scales govern the importance of finite-size effects, the mean free path and the range of in-
teraction. Microscopic models based on the Boltzmann equation, easily incorporate the
effects of a finite mean free path. Such effects can be linked to viscous terms in analogous

hydrodynamic descriptions.

In intermediate-energy collisions, where excitation energies are tens of MeV per nu-
cleon, the role of the finite range to the strong interaction has been studied in its relation to
the surface energy of nuclear matter. In such Boltzmann descriptions the binding energy of
nuclear matter is introduced via the mean field, with the coarseness of the mean field mesh
being adjusted so that the effect of the interaction range is effectively tuned to reproduce
the surface energy of nuclear matter. In nonrelativistic molecular dynamics, the effects of
hard-sphere interactions have also been investigated. In this case the size of the spheres
represents a length scale which can strongly affect bulk properties of the matter at high

density. However, in the context of a Boltzmann description, where n-body correlations are

87

explicitly neglected, the effects of a finite range inherent to the scattering kernel have not

been analyzed for their impact on final-state observables.

It is the goal of this study to ascertain the importance of this second length scale in
high-energy collisions. By varying the interaction range in the scattering kernel, while
leaving the mean free path unchanged, the consequences of a nonzero interaction range are
studied. It is demonstrated that a finite interaction range contributes to viscous terms in a
manner similar to the finite mean free path, but with different dependencies with respect to
density and temperature. The finite range of the interaction is found to affect the evolution

of heavy-ion reactions and to alter final-state observables, especially elliptic flow.

After formally reviewing the collision kernel in the Boltzmann equation, it is shown
how viscosities arise from the interaction range. the finite range. Viscous heating is then
explained by analyzing the collision kernel. In particular, a comparison is presented be-
tween the heating derived from analysis of the collision kernel and the heating observed in

a simple simulation.

5.1 Connecting Viscosities to Finite-Range Interactions

5.1.1 The Role of the Collision Kernel in Boltzmann Descriptions

The Boltzmann equation can be expressed,

a
affirm!) +vp - Vf(p.r,t) + F(r,t) -fo(p,r,t) = / d3qd3q’d3p’d3r'dz'
X {f (q,r,t)f(q’,r’,t’) 9C(r— r’J - t’;q,q’;p,p’)

-f(p, r,t)f(p’, r’,t’) 7dr — r’,t — t’;p,p’;q,q')}, (5.1)

where f is the phase space density and F is the force dp/dt felt by a particle at position
r. The collision kernel 7((r — r’ ,t -— t’ ;q,q’;p,p’) describes the differential probability for

scattering a pair of particles separated in spacetime by x — x’ with initial momenta q and

88

q’ into final states with momenta p and p’. The range of the collision kernel in coordinate

space is the subject of this study.

Integrating over the collision kernel should yield the cross section,

1 do

[d3r’dr’x(r-— r’.t —t’:q.q’;p,p’) = 22—70333—37vaka —q—q’), (5.2)
re

where pm is the relative momentum of the outgoing pair (p — p’ ) / 2. By inspection of
Eq. (5.2), it can be seen that the coordinate-space dependence of 17C appears rather arbitrary
as long as it integrates to the free-space cross section. Indeed, results at low density, where
particles interact only pairwise, are unaffected by the form of K as long as the range is
much less than the mean free path and much less than the characteristic dimensions of the

reaction volume.

Any nonzero extent of the collision kernel leads to problems with super-luminar trans-
port. However, these problems are easily defeated by restricting the kernel to being local,

i.e.,

 

I I I I I I do I
7((r—r,t—t ;q,qum): 53(r-r)5(t-t )-—-—-vre153(p+p -q-q’). (5.3)

l
(2703 d3prel

The Boltzmann equation can now be written in a manifestly covariant form.

1 d3q’d q (5 4)
(21:)3 qu Eq '

 

 

a
(#58qu Pug) f(p,r,t) 2

 

 

 

 

do do
- ’2—m4 ~ .r,t ’,r,t — - ’2—m4
{\/<q q) d3pmf<q , )f(q ) 7(1) p)

d3fIrel

f(p,r,r)f(p’,r,t)}

Here, up is the four—velocity of a particle with momentum p, F “ is the force dpu/dt, and

pm. is the relative momentum of the outgoing particles in the center of mass.

5.1.2 Effective Viscosities from Finite-Range Interactions

The interaction over a finite range contributes to the shear viscosity, 1], the bulk viscosity,

C, and the thermal conductivity, x. The range of the interaction is related to all three

89

coefficients. In order to make this connection, two particles are considered which scatter
from one another, separated by a distance r 2 r2 — r1. Combining this finite separation
with the velocity gradient, it can be seen that the first particle interacts with particles which
have a higher average energy. By evaluating the rate at which energy is transferred to the
first particle from colliding with more energetic particles, an expression is found for the
rate at which heat is deposited to the region defined by n. By comparing to analogous
expressions from hydrodynamics, expressions can be derived for all three viscous terms in

terms of the interaction range, r 2 In — r2], the density n, and the collision rate F.

Choosing a reference frame such that the velocity of bulk matter at the location of the

first particle is zero, the collective velocity at r2 is

8v,-

v,- 2Aijrj, A

For an elastic collision where two particles of identical mass simultaneously change
their momenta, the radial components of the momenta must be interchanged by the collision
if energy, linear momentum and angular momentum are to be conserved. Physically, this
corresponds to the scattering from the interior or exterior surface of a hard sphere. The

average energy change of the first particle is then
(512,): <E2.r> —<E1,.> = —<(v - 02>. (5.6)

The mass term m is not to be taken literally as the mass of the particles, since the averaging
may include factors of the velocity to account for the flux or it may have a complicated
form to accommodate a desired differential cross section. For relativistic motion, the mass
might incorporate the lateral motion of the particles. Writing (513;) in terms of A,

(551) ‘ £<(riAijrj)2>~ (5.7)

— 2r2
The average over the directions of r can be performed using the identity,

4
r
AijAkl (firjrkrl) = AijAkIE (5ij5k1+ 8ik8jl + 5115jk) - (5.8)

90

Then, (5E1) can be expressed in terms of A and r.

,2
(SE): "32.—0 ((TrA)2 +— ;(2 (AU-Mfg) (5.9)

U

:2: [(V v)2+ %z(§%+%:—:)z].

The rate at which the entropy increases due to these interactions is given by the density

multiplied by the rate at which collisions deposit energy nonlocally,

as = ETE<8E1> (5.10)

_nl"mr2 av, +812]
_ 307" [(V ”2+ 2Z(ar—,+ an)?

Here, F is the collision rate experienced by a single particle and n is the density.

 

It is notable that only the symmetric part of A contributes to (SE). This owes itself to
conservation of angular momentum which forbids rotational motion from being transferred
between particles. In fact, if an expression for (SE) had been derived using v2 instead of v},
the resulting expression would have included the odd parts of A which would correspond
to rotational motion, V x v. These terms would have no hydrodynamical analog as they

would have reflected a violation of angular momentum conservation.

Analogous expressions for Eq. (5.10) are now provided in the language of hydrody—
namics. The expression for entropy production [152, p. 55] in terms of velocity gradients

is

 

._:1_ Bra-€212.23»):
a S—ZT .. (8xj+axi_3s”v v) +T(V v)2 +fi(VT)2. (5.11)

The coefficients can be extracted by comparing Eq. (5.10) to the first two terms in 5.11.

 

 

mnI‘r2
2 5.12
T] 30 ( )
man2
2 . 5.13
c 18 < >

91

Applying similar reasoning, an expression for the thermal conductivity can be derived
as well. First, the temperature gradient is related to the energy flow. Again, particles are
considered to be separated by r. If collisions occur between two particles at locations with

different temperatures, the average energy exchanged is

l
SE1 2 5 ,r-VT. (5.14)

Here C, represents the change in radial kinetic energy per particle per change in tempera-
ture,

c, = 15,. (5.15)

37
In the nonrelativistic limit, C, 2 1/2.

Since the exchange corresponds to moving an energy a finite distance over an effective
time given by the collision rate, the average momentum density can be defined in terms of

the energy flow.

 

 

nFC, 8T
Ml_ —7—<rirj)a—x;. (5.16)
nrc,r2
2 — VT 5.17
12 ( )
An extra factor of l / 2 was added to correct for double counting the collisions.
The energy flow can be related to the entrOpy production,
(15 l as
—— 2 d3 —— .
dt / xT at (5 l8)
3 l
2 —/d'x—V-M (5.19)
T
nFC r2
= d3 ’ VT 2 5.20

where the continuity equation has been applied. This can be compared with the last term

in Eq. (5.11) to obtain the thermal conductivity x.

_nI‘Crr2
7‘“ 12 '

 

(5.21)

92

The forms for the three coefficients, '1], C, and x, fundamentally differ from the forms
that result from considering a nonzero mean free path If. The viscous coefficients that result
from finite Z are independent of density and scale inversely with the cross section [152]. The
coefficients arising from a nonzero interaction range r scale as the square of the density and
are proportional to the cross section. Thus, nonlocal interactions provide viscosities that are
important in wholly different conditions than those where the finite mean free path plays an
important role. For a rapidly expanding system, finite-range interactions play an important
role when the interaction range multiplied by the velocity gradient provides velocities of
similar magnitudes to local thermal velocities. Such conditions exist in the first one fm/c

of highly relativistic hadronic collisions.

The expressions here derive from a very specific picture, nonrelativistic particles mov-
ing on straight-line trajectories punctuated by sharp collisions when the separation equals r.
However, all interaction at a finite distance should result in viscous behavior. Relating the
distance r to the cross section might involve a detailed microscopic evaluation of the col-
lision kernel. This is especially true for relativistic motion. Despite the complications, the
distance r2 can indeed be expected to be of similar magnitude to collision cross sections.
It would be interesting to discern whether the ratio between the shear and bulk viscosities

varies for different scattering models.

5.2 Analyzing the Collision Kernel

In Subsection 5.1.2 a finite range of interaction, characterized by a length scale, r, was
shown to generate viscous behavior. This length scale is determined by the cross section,
r2 ~ 0 / 1:, but the constant of proportionality is not trivially determined and can depend on
seemingly arbitrary aspects of scattering algorithms. In the following subsection, we illus-

trate how viscous heating can also be directly related to the collision kernel by considering

93

a simple example of one-dimensional expansion with a Bjorken spacetime geometry. In the
subsequent subsection predictions based on the form of the scattering kernel are compared

with results from simulations based on the same kernel.

5.2.1 Viscous Heating in a Bjorken Expansion

As a simple example, a one-dimensionally boost-invariant system of infinite extent in all
dimensions is considered again. The number of particles is kept fixed such that the density
scales inversely with 1,

72(1) = $d2N/(dAdn). (5.22)

where the number of particles per area per rapidity interval, dZN/(dAdn), is fixed. To

simplify matters particles are considered massless.

Here, the rate at which a particle at n 2 O has collisions with other particles, ch /d1:,
and the rate at which it gains or loses energy from such collisions, d (E 1) / dt, are calculated.
The particles are assumed to be in local thermal equilibrium, and the collision kernel is

analyzed to determine the rates.

Referring to the colliding partners with the subscript ‘2’, the collision and heating rates

per particle are

 

 

 

 

ch d2N 3 o 2 d3N
d’t _ ZRdAdT] [(1112 rzdrzd p2 5(1'C —‘L') 9 (II— -b )m, (523)
d<E1> dzN / 3 O 2 d3N
22 d' d d 5 c— 9 ——b 5E —. 5.24
dt ndAdn T12 r2 r2 P2 (T T) (n ) 'd3p2 ( )

Here, SE] is the average change in energy of the particle ‘1’ due to the collision and r2 and
n2 describe the position of the second particle. The impact parameter for the two-particle

collision is b and IC is the collision time.

The strategy employed here is to calculate the collision time, IC. The first particle is

at a position x1 2 (t,0,0,0) with four momentum p1. The second particle is at a position

94

x2 2 (tcoshnz, r2,0,tsinhn2), with a four-momentum p2. Once 16 is understood in terms

of r2, 712, p1 and p2, the delta function can be replaced in the expression above,

gr:

at}, (5.25)

5(TC - I) —-> 5(rc - r2)

where rc is the position required to make the collision occur at I. By substituting the delta
function with re for the delta function with To the integrals in Eq. (5.23) and Eq. (5.24) can

be simplified and solved numerically.

The first step that must be performed is to find “cc in terms of r. The prescription
for finding Tc is somewhat arbitrary due to covariance issues as discussed in 3.2.5.2. A
prescription must be arbitrarily chosen to find “cc in terms of £1 and 2:12. For the purposes
of this study, the prescription according to Eq. (3.51) is employed to simplify the algebra.

Using this choice,

CH2 2Cr’ C
fizfi+° +2‘+2 5%)

2
C0 = YIPIJ. +Y2P2I

 

C1 = Y10£1+ SIYIPII + Y20‘2 + 92Y2P21

The coefficients are defined

(11 = El! ’ [91,32 5 (12 Z E21+P2,;Z (5-27)
51: 2ZP2..~/(P1 192), 52 = -22P1.z/(P1'P2)

Y1 = 2P2.x/(P1'P2) , Y2 = -2P1,x/(P1'P2)

Here t, r and z are the positions describing the first particle in a frame where the particles

are centered about 1] 2 0 and r 2 0.

71,: 1‘23, I 2 tcosth’, Z = 2TSinhTII, r' = r32 (5'28)

95

In terms of these new variables, the rates above can be rewritten,

 

 

 

 

dN d2N 3 o 2
d1 —1rl6 dAdn dZN’fdnr ”dr d p2—8(t— tc)@ (fi_b ) (5.29)
21 o
216 d ’ ’d ’d3 —5r ————9 ——b2 5.0
d<El> 11] I I 3 21 2
--——-—21 d rd ’d --————8 ___—8 —--— E 5.
dt 671212122111 11 ' ”(132p ('J’ rC)C0r’+C1 (:5; b)5 1 ( 31)

When solving for re, there are two solutions to Eq. (5.26). If solutions for both positive and

negative rc are used, the factor 161! should be reduced to 81:.

These results can be expressed as averages over p1 and p2

 

 

 

 

 

dNC _ d2N , r. o 2
d(El) d21v 2
= 32 d2 ' —8E — — b . 5.33
dt mdAdn n <C0rcr: +C1 '>@ (n ) ( )

If scattering angles are chosen with equal probability forwards and backwards, the average

change in the two energies is
E2 — E 1
2 .

 

SE] 2 (5.34)

The nonlocal aspect of the collision kernel should contribute to heating the particles in
their local frame. From physical arguments, the nonlocal contribution to heating is expected

to scale with temperature, density, time and cross section in a simple manner.

arr:i _ 3 _d(E1)B d2N 027“
dt 4n d1? dAdn T3 ’

 

(5.35)

where E, is the mean transverse energy per particle and B is a dimensionless constant.

The simple scaling derives from the hydrodynamically motivated Eqs. (5.11) and
(5. 12). The collision rate is expected to scale proportional to the density, which requires
the factor dzN/ (dAdn) and one inverse power of “c, and the cross-section. The squared-

velocity gradient suggests an extra factor of 1‘2, and the range of the interaction requires

96

an extra factor of o. The constant B is determined by the form of the differential cross sec-
tion, e. g., s-wave scattering would result in a higher B than a highly forward-peaked form.
Since one power of 0 comes from the range of interaction, B should scale inversely with the

sampling factor 3.. Finally, the effective mass should be proportional to the temperature.

 

 

 

 

 

UTIITrTTIIIIIIIIIIT
0.03 -------------
C)
”(3'— - 0 ~
'0 N O O
\2 o
/\ >\ ”o 'J
4‘0
“-1 <
V ‘0 -o _
‘0 \
n N2
t~ ‘0 ~ 0:10 mb-
Vv
0.02 - —
OllljllliilllJJLLlLl

 

O 10

5
T (fm/c)

Figure 5.1: The heating due to nonlocal interactions as calculated numerically from the
collision kernel (circles) is scaled in such a way that it would be constant if the simple
hydrodynamic scaling arguments were valid. The dashed line represents the asymptotic
value.

The heating rate due to nonlocal interactions is calculated from Eq. (5.33) by numer-
ically analyzing the collision kernel and is displayed in Figure 5.1 after being scaled by
the temperature, time and cross section. Had the simple scaling argument been correct the
ratio would have been a constant B. However, due to higher order corrections in 1 /‘t, the
ratio varies as a function of 13. The ratio approaches a constant for large 12, “£2 >> 6. The
scaling fails when 1: becomes less than the range of the interaction, which for this example

is ~ 1.0 fm due to the cross section of 10 mb. The departure of the ratio from a straight line

illustrates the limitation of simple hydrodynamic arguments to describe the viscous heating

97

from nonlocal interactions.

5.2.2 Comparison to Numerical Results

Finally, numerical results are presented from a model based on the Gromit framework. A
boost-invariant and in the transverse directions closed system are simulated with cyclic
boundary conditions in n and in the transverse coordinates, respectively. The viscous heat-
ing rate observed in the numerical calculation are compared with the viscous heating rate
expected from the scaling arguments expressed in Eq. (5.35), where the coefficient B was
determined from analyzing the collision kernel from the last section. The temperature was
set to 400 MeV at a time 1.’ = 0.1 fm/c, and the cross section was chosen to be 10 mb. A

simple s-wave form was used for the angular dependence of the cross section.

 

4 I IIjTWITI I Tllllll' j FF

\\N.S O GROMlT—T _

.3 — \ -
\

.. 00\ -l

N
l

\ o

"‘..‘E3jorken

         
  

+non|ocol

    

 

 

Ei('r)/Ei(25 fm/c)

 

1 11111111 1 llllllll

"‘ 10° 10‘
T (fm/c)

Figure 5.2: The mean transverse energy for particles participating in a boost-invariant
Bjorken expansion is displayed as a function of the proper time (circles). Also displayed
are the Bjorken hydrodynamic result (dotted line), the Navier—Stokes correction which ac-
counts for viscous shear arising from a finite mean free path (dashed line) and the correction
due to nonlocal interactions as expected from simple scaling arguments (solid line). Nonlo-
cal interactions are important at small times when the velocity gradients and collision rates
are high.

 

98

The resulting mean transverse energy is displayed in Figure 5.2 as a function of ’c. The
initial heating derives from the nonlocal nature of the interactions. Longitudinal cooling
ultimately dominates the behavior as the nonlocal contribution to the heating falls roughly

atr‘3.

The dotted line in Figure 5.2 describes the evolution of the transverse energy in the limit
of ideal (nonviscous) hydrodynamics. In that limit, the stress-energy tensor has a simple
form,

TC“B = tEuO‘uB + P(u°‘uB — gaB). (5.36)

In the Bjorken limit, av/Bz = 1/1: and the evolution of the energy density becomes

 

a P+ e
—e = — . 5.37
at 1' ( )
For the massless case, P : 6/3, which gives the result,
1 -—4/3
8(1) : 8(1'0) (—) (5.38)
To

The velocity gradient generates a shear which contributes an additional term to the

stress-energy tensor.

1 Ba all 2
TéB=n(—gaY—uazlY)(§é+5£—§ $38-14). (5.39)

In a simulation, Tle can be determined by evaluating the stress-energy tensor.

T” T...
_Ll _ TZZ : 2%. (5.40)

2 r
The components of the stress-energy tensor are extracted by analyzing the momenta of

particles as measured in the local rest frame.
1 Pin
T--=— _, 5.41
r, V 2, Ep ( )

The Navier-Stokes evolution of the energy density is governed by the equation,

a P-l-E 4nNS

—8— —.
a1.- 1 3 12

 

(5.42)

99

For massless particles interacting with a constant cross section, dimensional arguments
force the viscosity to rise linearly with I since the mean free path should grow with I due
to the density falling as l / I.

7le = CNSET, (543)

where CNS is a constant, determined only by the form of the cross section. The energy
density then follows the Navier Stokes form [85]. This form should be valid unless the vis-
cous contribution to the stress energy tensor approaches the same order as the equilibrated

pressure [52].

1: ~(4/3)(l—CN5)
8(I) = 8(I) (a) . (5.44)

The Navier Stokes result is represented by a dashed line in Figure 5.2. The value of CNS
was determined by evaluating the asymmetry of the stress-energy tensor in the simulation at
large times. Running simulations with a large sampling ratios generated results in excellent

agreement with the Navier Stokes result.

The inclusion of nonlocal effects is responsible for the discrepancy between the sim-
ulation results in Figure 5.2 and the Navier-Stokes results. Eq. (5.12) suggests that the

nonlocal correction to the viscosity should scale proportional to 1'2,

nlrl : CINE/1:7 (545)

where C”, is independent of I and scales with (32 as explained in the previous section.
Now, the evolution of the energy density can be determined by solving the hydrodynamical

equations of motion,

a 4 8 4 8
— =—— l— — —C,,— 5.46
BIE 3( CNS)T+3 113 ( )
This equation has a simple solution.
4/3( l_C/VS)
T 2Cn1 l l
: — _ . 5.47
8(7) 8(TO) (To) exp( 3 (T8 12)) ( )

100

This form is shown with the solid lines in Figure 5.2. The constant C", was determined
by the asymptotic limit of Figure 5.1. The effect of nonlocal interactions is somewhat
overestimated at small times by Eq. (5.47) as would be expected by considering Figure 5.1

which shows that that growth of the viscous heating at small I is somewhat slower than the

 

 

 

naive expectation that it scales as I’3.
A 4 I II‘IIITTI I lllll'll I II
0 ‘.
\ _ o GROMlT—‘r _
E ‘-
H— 3 ._ _.
L0
9', _\\N.S. .. _
—r o yBJorken
LrJ 2 — O .. —
\ 0cm©\
A _ 9 '. _
r~ 9
v '
_l 1 _— N_
LL] 1 r r rrrrrl r 1 1 erul 1 L

 

“ 10° 10‘
T (fm/c)

Figure 5.3: The mean transverse energy is plotted as a function of time for the case where
a forward-peaked cross section as described in Eq. (3.35) is implemented. The results
from Gromit-I (circles) are well described by the Navier-Stokes (dashed line) correction
to the Bjorken solution (dotted line). As compared to s-wave scattering, the Navier-Stokes
viscosity is increased while the nonlocal contribution to the viscosity becomes negligible.

10

The effect of nonlocal interactions is lessened for scatterings that are more forward
peaked. Figure 5.3 illustrates the behavior of E, as a function of I for a screened Rutherford
scattering as given in Eq. (3.35). Here, the screening mass It is chosen to provide a cross
section of 10 mb. As compared to the s-wave scattering result in Figure 5.2, the effect of
nonlocal interactions is diminished while the Navier-Stokes viscosity is increased. This is
expected since the mean free path is effectively increased while the energy transfer inherent

to collisions is decreased.

101

Finally, the effect of nonlocal interactions that might be expected in the earliest mo-
ments of the big bang should be mentioned. In that case cross sections should become
perturbative and particles should be approximately massless. In this case cross sections
should scale as 0.2/T2, where or is the unified coupling constant. Since T would scale as
l /I, the Navier-Stokes viscosity and the nonlocal viscosity should both scale identically
with I. Since the nonlocal terms scale as 02 while the Navier-Stokes terms scale as 0'“,
the terms would differ in importance by a constant proportional to a6. Thus, if the sys-
tem becomes perturbative, the nonlocal terms are expected to be negligible compared to
the Navier-Stokes terms. However, the physics of high-energy hadronic collisions is far

from perturbative, and the large coupling constants magnify the importance of the nonlocal

terms relative to Navier-Stokes terms.

102

Chapter 6

Correlation Function for Particles from
Resonance Decays

The Hanbury Brown—Twiss (HBT) effect was first proposed in 1954 as a means to measure
the angular size of distant stars by exploiting the quantum correlations of photons [88]. The
same authors then demonstrated their technique in 1956 by determining the angular size of
the star Sirius [89]. Independently in the early 1960’s, the I-IBT effect was discovered in

proton—anti-proton collisions, by Goldhaber et al. [79].

Since then, analyses of two-boson correlations have provided tangible information re-
garding the spacetime development of hadronic reactions [91, 155]. Pions, kaons, and pho-
tons have all been exploited for their bosonic nature which results in a positive correlation at
small relative momentum. Numerous other correlations, involving nucleons [90, 21, 70, 75]
or light nuclei [25, 26] that are correlated due to the strong or Coulomb interaction as well
as identical-particle statistics, have also been analyzed and have given further information
regarding collision dynamics. Source sizes and time scales have been extracted from colli-
sions covering a wide assortment of reactions, from heavy-ion collisions at a few MeV per
nucleon, where time scales of thousands of fm/c have been determined, to e+e‘ annihila-

tions, where lifetimes of a fraction of l fm/c have been observed [78, 19, 47, 3].

The comparison of theoretically predicted correlation functions with experimental re-

103

sults provides an important test of the dynamical properties of reaction models. Most
models provide semi-classical information about the source function S(p,x), the proba-
bility of emitting a particle of momentum p from the spacetime point x. By convoluting
the source functions for particles of momenta p1 and p2 with the squared relative wave
function |¢q(x1 — x2)|2, one is able to predict the correlation function C(p1,p2). Source
functions are usually obtained from semi-classical simulations, where the source points are
associated with the last point of interaction [126]. Particles from resonances are usually
assumed to be emitted according to an exponential decay law, with the characteristic time
usually chosen to be independent of the energy of the resonance. Quantum considerations
have been explored by Lednicky and Progulova [106] and by Bertsch, Danielewicz and

Herrmann [30].

In this chapter, two goals are pursued. First, the importance of quantum treatments is
quantified by comparing to semi-classical forms for a simple thermal model. Although the
presented formalism is not much different from that discussed previously in the literature
[106, 30], the differences with semi—classical treatments have not previously been studied
quantitatively. Quantum corrections are found to become important when kinematics con-
strain the resonances to be off mass shell. Second, an alteration is proposed to the methods
for extracting correlations from classical simulations to better account for quantum effects.
This modification is shown to be able to account for the quantum corrections by incorporat-
ing information regarding the off-shell energy of the decaying resonance. In this study, any
interaction between the particles aside from the constraints imposed by symmetrization is

neglected.

104

6.1 Derivation of Corrections to the Correlation Function

6.1.1 Correlations from Direct Sources

The two-particle correlation function for particles with momenta p1 and p2 is usually [91,

155] defined as
dN(2)/(d3prd3P2)
(Wm/£13m) (dN“’/d3p2)°

 

C(pr,p2) E (6.1)

The one-particle probabilities in the denominator can be determined by matrix elements
Tf (x) where f describes the remainder of the system, andx is the point at which the particles

had the final interaction with the system. Without loss of generality

 

 

dN“) . . 2
f
The definition of the source function is
S(p,x) E 2 [d45x T;(x+ sx/2)rf(x —— 5x/2)e_ip'&, (6.3)
f
which leads to the simple relation
dN“)
2Ep : [d4x S(p,x) (6.4)
d3p
POZEp

 

The source function can be interpreted as the probability per unit spacetime for creating a

particle of momentum p.

Since source functions can be extracted from semi-classical simulations or thermal
models, it has proven useful to also express two-particle probabilities from Eq. (6.1) in
terms of source functions. The two-particle probability requires a two-particle matrix el-
ement Tf(2)(xa,xb). With the assumption of independent (or uncorrelated) emission, the

two-particle matrix element factorizes [30].

T}2’(xa,x,,) ——> Tfa(xa)be(xb), (6.5)

105

where, a and b label independent sources. The two-particle probability then becomes

de

2E 2E ——

2

=2

aibrfflrfl)

iPr -Xa+iP2'Xb + eiP2°xa+iPr 2%}

 

 

/ d x0614 xb Tfa (xa) be( X),

”We

: 2{‘/d4x(1d4x’) SG(pl,xa)Sb(p2’xb) + [d4xad4Xb Sa(p-,XQ)Sb(p_,Xb)ei(p2’P1)‘(Xa—Xb)},

where [9’ 2 (p1 + p2) / 2. The expression above can be rewritten conveniently in terms of

the Fourier transform of the source function

 

101,691,122) = [d4x 50,1, (p‘ :m ,x) e‘1P2“P'>*, (6.6)
which leads to
(1M2)
(2Er)(252)m = 2{Ia(pr,pr)lb(p2,p2)+la(pr,p2)1b(p2,pr)}- (6.7)
a,b

The correlation function then takes on a simple form

02 1a (P1,P2)1b(P2,P1)
(6.8)

 

C - l
(P13P2)_ +0 2210 (p1 p1)Ib(P2aP2)

which is similar to that originally mentioned by Shuryak [139].

Simulations typically provide a sampling of the on-shell source function. The applica-
tion of Eq. (6.6) in simulations is made difficult because the source functions are evaluated
at p’o ;é E( Pr+ p2 ) / 2 meaning that they require off-shell information. The above formalism

can be related to simulations through the smoothness approximation [125, 7],

+ +
50(P1 p2,xa) 5b (pr mm?) —9 5a (P1,xa)5b(P2,xb)- (6-9)

 

 

2 2

For thermal sources, a form for the off-shell behavior of the source function can be justified

and the smoothness approximation can be averted. A quantum theory would provide the

106

T matrices that would also allow off-shell evaluation of the source function. In fact, some
formalisms have been developed where classical simulations are augmented by converting
the point particles into wave packets [119, 164]. This also allows to forego the smoothness
approximation, but at the price of inserting an ansatz for the quantum behavior that has
some peculiarities. The issue of the smoothness approximation will be side-stepped in
this study, so as to focus on quantum aspects associated with the propagation of off-shell

particles.

When calculating correlation functions from simulations, particles from resonances are
usually included in the source function by using the spacetime points from which the res-
onances decay. The lifetime of the resonance affects the correlation function through the
exponential decay that is simulated in the transport model. As will be demonstrated below,
the exponential decay law is modified if the dynamics emit resonances with a particular
mass or range of masses. In this case, the form of the source function becomes nonexpo-

nential as will be explained in the following.

6.1.2 Correlations from Resonant Sources

A formalism is presented for calculating two-particle correlation functions from resonance
decays given that one or both of the pions might originate from a resonance. The result
will depend on the source function that represents the creation of the resonance rather than
the source function that represents the points at which the final-state pions are created. The
evolution and decay of the resonance will be accounted for through the quantum propagator
of the resonance. The spacetime point at which the resonance decayed to produce the
final-state pions will be treated as an intermediate quantum step in the evolutionary path
between the initial creation of the resonance and the asymptotic momentum states of the
decay products. The expressions derived here thus incorporate an integration over all points

at which the resonance might have decayed relative to the points at which the resonance is

107

created.

Three examples will be considered, a hypothetical scalar particle A that decays into
two pions, the vector meson p, which also decays into two pions, and the vector meson
(t), which decays into three pions. The matrix element for creating the pion and a second

particle with momentum k through a scalar resonance A is

Tn,fA.k(r) = g/d4xA GA(x—.rA)TA,L,(xA)eik x, (6.10)

where GA is the Fourier transform of the propagator for the resonance and fA refers to
the state of the remainder of the source. Here T,t is effectively the T-matrix element for
emission of the pion, while TA is the T-matrix element that would describe emission of the

resonance if the resonance were stable.

Following the same method as in the previous subsection, Tn can be used to create the
pion source function using Eq. (6.3) to express TA" (x) TA (1” ) in terms of the source function
of the resonance. The resulting expression for the pion source function can be used to

generate IA, as defined in Eq. (6.6), which is sufficient to calculate correlation functions,

 

d3k . 1+ 2
1A(PI,P2) =82/QE/d4x/r €XP11(P2-Pr)'XAlSA(p 2 p +k,XA)

 

x cm +k)G:r(P2+k)a (6.11)
i
G = , 6.12
o M
IIA(m“) I MAF,\i-—A. (6.13)
cu; m

Here, m2 2 p3,, MA and FA are, respectively, the pole mass and width of the resonance. The

2

relative momentum of the outgoing pions in the frame of the resonance is q2 = "22/4 — m,[

and q}? is the same quantity for an on-shell resonance. It should be emphasized that the
interference term in Eq. (6.11) can be calculated without reference to the direct source of

the pions, as the source function of the resonance becomes the required input.

108

In order to understand the role of the propagators, the emission of a pion pair with mo-

menta p and k through a single, scalar resonance with momentum pA = p + k is considered

dN

 

 

 

   

2Ep-d—3p: [A (p, p), (6.14)
d d3k 2
dN— z—lE—pz—E 11% S (p+k,x) 2 g. 2 (6.15)
k |(p+k)2-MA+1HA((p+/<)2)l
d3“ S( )d 211m 1 (616)
= x ,x m - . .
2EA A pA 7t mZ—Mi—inMpE‘)

The spectral function that describes the density of states of resonances of invariant masses

m,
l 1 dn
——Im 7 . 2 = —-2—.
211 m2 "MZ —11_IA(pA) dm

 

(6.17)

Thus, it is evident that the source function does not provide any information regarding the
mass or width of the resonance. However, combined with the spectral function, which de-
rives from the product of propagators, it provides the probability of creating the resonance
at space time point x with momentum p and with invariant mass m. For the direct case
considered in the last subsection, the source function was always evaluated on shell, i.e.,

the spectral function was effectively a delta function.

The same calculation can be performed for vector resonances such as the p meson. In

that case, the coupling of the pions to the vector meson is
Lin! 2 i}t(n8u1t) p". (6.18)

The expression for 19 is similar to 6.11 with the source function and propagator accounting

for the vector nature of the p,
’1) (P1 P2) =4 /%k/d4xp 61113111172 -pr) 161)](171— k)aGgB(Pr +k)

x5"

P1 + P2
BY

2 + k,xp) GEY8(P2 + k)(Pz - k)s- (6.19)

109

The propagator for the vector resonance is

 

_ dB 01 B 2
G 32%) i 28 grip/1’2 . (6.20)
P _MA+lnp(P)

For the vector case the self-energy scales differently as a function of the resonance mass

than in the scalar case,

M
mm )= Mpr‘pgu—p (6.21)
R

where the same notation as in Eq. (6.13) was used.

As a third example, the propagation is considered of an (1) meson, which is also a vector

resonance, but decays into three pions through
. 0 _
Lint : zkewww"a"n+8§n awn . (6.22)
In this case the expression for 1w becomes even more complicated than the p example.

1( )—K dBkd—Bl/afix ex[( — )x]
to P1 P2 TEkZE, m P P2 P1 (1)

 

2
x6235(p2 + k +1)25W§, p12“ kV’zé’. (6.23)

+
xeawép‘l‘k 1§G$B(pr+k+1) 60,0“ ”Jr/(+13%)

The expression for the self-energy is also somewhat more complicated.

d3k d3l 2 2 2 2
11mm )—B/———8( (m—Ek—El) —Ek+|—2klcose)m |kxl|. (6.24)
215,. 2E,

where 0 is the angle between R and I and B is an uninteresting constant fixed by the con-
dition that II(,,(M(20) : Mwl‘w, with Mu, and I}D being the mass and the width of 0), respec-
tively. After applying the delta function, there remains a fairly complicated expression,

which was evaluated numerically in the subsequent considerations of the to case.

6.1.3 Correlation Functions from Monte Carlo Simulations

Simulations of heavy-ion collisions usually provide the creation points of pions along with

their outgoing momenta. Neglecting other interactions besides symmetrization, correlation

110

weights for pions of momenta p1 and p2 originating from spacetime points xa and xb are
usually determined by calculating the average symmetrization weight of all pairs satisfying

the imposed binning or acceptance [126].

2%)fd4xad4xb Sa(P1,xa)5b(P2,xb)Wa(P1,P2)Wb(P2,P1)

 

C , 2 1+ 6.25
(Pl p2) bed4xad4xb Sa(pl ,xa)Sb(p2,Xb) ( )
E 1+ (Wa(P1»P2)Wb(P23P1))- (6-26)
Inspection of Eqs. (6.6), (6.8), and (6.9) reveals the weight for direct sources,
$94812»
w , : ex i — ox ——"———, (6.27)
d(Pl p2) P1(P1 122) a] 509”")
wff” (P1 ,p2) = exp[i(pr — 1221de (6.28)

where the difference between the two forms is the assumption of the smoothness approxi-

mation [125,7] in the Eq. (6.28).

If the particle originates from the decay of a scalar resonance A, the weight takes on a
different form as can be seen from inspecting 6.11,

SA(fliz_fl+k,x/\) (pl +k)2_M/2\ “in/11(Pl +k)2] ei(pl-—p3)'1A (6 29)
5.0241,“) (p2+k>2—M3,—in11(p2+k)21 '

 

v121071.192):

Here the resonance was created at xA and decayed into pions of momenta p1 and k. The
spacetime coordinate of the decay does not enter the weight as all decay points have been
considered. If the decay is of a vector resonance such as a p meson, the weights are slightly

different,

50(P1,P2,k) (PI + ’02 — Mg-inp[(p1+ ’02]

 

 

__ i(Pl-P2111)
W 121,122 -— . o 6’ (630)
M ) 56(P13Prrk) (P2+k)2-M5-lnpl(P2+/<)“l
+
5p (pr.p2.k.xp) = (P1 — k)0L 5813 (p1 p2 +k,Xp) (P2 — k)B. (6.31)

In the derivation of SD the two pions involved in the resonance decay are assumed to have

equal mass.

111

Any sort of resonance can be included in this manner, including resonances that decay
into three or more bodies. One such example is (1), which decays into three pions, one of
each species. Labeling the momenta of the two pions, whose symmetrization is ignored, as

k and l, the weights turn out to be

W(1)(Plep2) = exP[i(pr —- p2) ~xm] (6.32)

X510(PlaP23kal) (171+k-H)2 —M120 - intqur +k+l)zl
Sw(p11p11kvl) (p2+k+l)2 —M&)-iI—Iw[(p2+k+l)21,

 

 

| + I I I
S(I)(plrp21krlrx(r)) : 81V: p’llkvlé (;B (p1 2 p2 +k+lrxm) EB#IVI§I Pg kv I: ' (633)

Thus, weights can be used to calculate correlation functions for the decay of any reso-
nance in a rather straightforward manner. The formalism coherently accounts for all points
at which the resonance may have decayed, but requires information regarding the points at
which the resonances were created as well as information about the accompanying particles
in the decay. The only difficulty comes in assigning the ratios of the source functions, i.e.,
the smoothness problem. For a thermal source, the source functions assume a simple form,

in the scalar case

SAM) ~ poe‘PO/T (6.34)

and in the vector case

SSMPaX) N Poe—pO/T(—801B + PaPB/P2)- (6-35)

The Boltzmann factor cancels out of the ratios. Thus, for the thermal example, the weights
become a product of four factors, a phase arising from the points at which the resonance
R (= A, p, (1)) is created, a ratio of energies, a spin factor, and a ratio of propagator denom-

i nators.

(mgr; +k’) mphpz) (P1 +102 -143, —irIR1(pr +1021

(P1+k')'" xR(pr.pr) (p2+k')2—M,2.—irIR[(pz+k’)zl
X ei(P1-P2)-IR

 

h
w); )(Prer) =

(6.36)

9

112

where It refers to the frame of the thermal source and k’ is equal to either k for a two-body

decay, or to k +1 for a three-body decay.

The spin factors xR depend on the sort of resonance being considered.

XA(P1»P2)= 1 (6.37)
Xp(P1,P2) = (P1 — 10%-8613 + PaPp/PZ) (P2 — ’06, (6.38)
x..,(pr.p2) = cam pl‘k"!é (ms + para/1'52) 85M; 125‘ché , (6.39)

where [3 =(P1+p2)/2+k’.

However, this technique cannot be easily applied to nontherrnal sources because the
off-shell behavior of the source function is not always known. This problem also confronts
calculations with direct sources, and forces either the invocation of the smoothness approx-

imation, or the assumption of some form for the off—shell behavior of the source function.

In analogy to the smoothness approximation, the thermal weight wgh) might be used

assuming a reference frame n" or simply neglect the ratio of energies in 6.36. In fact,
the ratio xR(p1 , p2) /X,R(Pl , p() can also be neglected to a reasonable approximation as the
structure of the correlation function derives largely from the last factor in Eq. (6.36), the

ratio of propagator denominators.

6.2 The Importance of Proposed Modification

Calculations of the quantum type described in Subsection 6.1.2 should be compared with
semi-classical calculations where the resonance is assumed to propagate classically and de-
cay according to an exponential form exp(—t / I). For the purposes of this comparison, two
Simplified systems are modeled, one of decaying p resonances and a second of decaying
(D’s. For each case resonances are produced and decayed with a Monte Carlo procedure

aCCording to a thermal distribution characterized by a temperature of 150 MeV. The mass

113

of p is 770 MeV and the width is chosen to be 150 MeV, while the mass of (1) is 783 MeV

and the width is 8.4 MeV.

6.2.1 Formalism for Semi-Classical Models

In the semi-classical descriptions, the correlation weights are determined by calculating the

expectation value
W'SC)(P1,P2,XR) = (dim—”2)” (640)
= €XP1i(P1- P2) 'XRl /d4(x ”XR) 531x — XR — "R(t — t1%)]
. l
X €XP11(P1 - P2) ' (X - ””71; (”(1314t — tR)/(YRT)1('3(t - IR)

mR/IR
mR/TR + iPR ° (P1 — P2),

 

=exp[i(pl —p2)-xR] (6.41)

where YR is the Lorentz factor due to the motion of the resonance. Here, wlsc) assumes an
exponential form for the pion emission, which is characterized by a lifetime I. The same

form for exponential decays was developed by Padula and Gyulassy [118].

Prescriptions for the energy dependence of the lifetime I(m) are arbitrary and can be
chosen freely. Three possibilities are investigated: (1) The lifetime is chosen such that
m/ I = [I(mz). This choice would be motivated by the form of the propagator. (2) The
lifetime is chosen to correspond to the average emission time as described in Subsection
6.2.3, except that the relativistic generalization of Eq. (6.54) is used I = 2m 8(m2 — M121 -

in)“. (3) A fixed lifetime l/F is used.

If the resonance distribution function in coordinate space is independent from that in
momentum space, the interference term in the correlation function factorizes into a term

stemming from the spacetime extent of the resonance source itself and one that arises from

114

the decay process.

C(p13p2)- 1 = 9i{<expli(pr - p2) - (XA —x3)l> (C'(pr,p2) - 1)} (6.42)

Here, xA and )1}; refer to points at which the resonances are created, and (exp[i(p[ —— p2) -
(xA — xg)]) represents the weighted average using the product of the source functions as the
weight:

(CXPMPI - P2) ° (XA — x3)]) = IJR(P1 — P2) l2, (6.43)

where

§Id4xASAKPr + P2)/2 + kaAleXPlKPr — P2) 'xAl
;fd4X/15AKP1 + P2)/2 + 191,901]

 

JR(PI — P2) = (6.44)
The reduced correlation function, C’ (p1, p2) — l, is similar to the average of the weights
<WaWb> from Subsection 6.1.3, only with the factors of exp[i(p1 — p2) -xR] removed from
the weights in Eq. (6.36). Since C’ contains all the relevant information about the decay,

we will focus on the reduced correlation function for our comparisons.

6.2.2 Numerical Results

To model the uncorrelated emission of pion pairs, thermal particles of momentum k (and l

for the (1)) are created by Monte Carlo, then a weight is added to the particles

7_ ___ Ek{+EI} +Er XR(PiaPi) (6 45)

EdEzlEr |(p.- +k’)2 “M121+inl(l’i+kl)2]12,

 

where i : 1,2 and the braces indicate terms that appear for (1) only. This weight accounts for
the spectral function of the resonances as described in Eq. (6.17) and for the spin factors XR
as in Eqs. (6.37)-(6.39). Note that the weights 21,2 are merely used to generate resonances
and their products, and are not related to the correlation weights. If these weights were
included through a keep/reject prescription, they would not need to appear in any of the

following expressions.

115

 

Once the pions are statistically generated, the average weights described previously are

calculated to generate the correlation functions,

a h ,b
ZZ‘fWgh)’ (P1,P2)§Zgwiet ) (P2,Pr)

" 6.46
y], ézg ( )

€021,192): H91

For this comparison wR is either the weight for p, wp, or the one for (1), wm. Note that a and

b label individual resonances.

 

 

 

 

  

 

 

 

 

 

1.0 ' I I I I I I I I I I I I I I L
-Thermol p Source .
- :, T=150 MeV‘
0.5
|
A
O“ 0.0
v 1.0
\ quantum ‘
O _ 0 0 0 0 semi—closs.(1) 4
~ 0 U D D semi—closs.(2) 4
- A A A A semi—closs.(3) ‘
0.5 - ‘
1.
L
0-0 IJLLIIIIJLIII
0 250 500 750

q (MeV/C)

Figure 6.1: Reduced correlation function for a source of p mesons at a temperature of
150 MeV. The average momenta of the two pions is fixed at 200 MeV/c. By factoring out
the spacetime dependence of the p, the manifestation of the p lifetime is singled out. The
exact quantum treatment is shown to differ from the three semi-classical treatments that are
described in the text. The failure of the semi-classical descriptions arises from the fact that
the p width is sufficiently large for the thermal source to effectively emphasize a specific
region of off-shellness.

116

Subsection 6.2.3.

C /(q)-1

1.0

0.5

rho
00

0.5

0.0

Figure 6.2: Same as Figure 6.1,

 

      
      
  

I I I I I I I I I I I I I I L
hermol p Source _

T2150 MeV:

 

 

 

 

 

I 5. quantum _
r 0 O 0 O semi-class.(1) ~
1 D D D U semi-closs.(2) ‘
h A A A A semi-closs.(3)_
1....11 1 r . 1.1
0 500 1000 1500

q (MeV/C)

except that p’s have aVerage momentum P = 800 MeV/c.

Figures 6.1 and 6.2 display the reduced correlation function for the cases where (Pl +
p2) / 2 = 200 MeV/c and 800 MeV/c, respectively. The upper panel of each figure displays
results for the case where p; is parallel to p2 while the lower panel displays the results for
the case where the two momenta are perpendicular. All three semi-classical results exhibit

significant deviations from the quantum calculations.

Figures 6.3 and 6.4 display the same information but for a thermal source of (1) mesons.
In this case, the semi-classical treatments are significantly more accurate. This could be
explained by the width of the (1) being much smaller than the temperature, which allows

the spectral function of the (1) to be sampled evenly. This notion will be explored further in

The overall width of the correlation functions in the upper panels of Figs. 6.3 and 6.4

117

 

 

 

  
   
 

 

 

 

1.0 I r T I I I I I I I I I I 1 L
~' hermal a) Source -
; ' T=150 MeV j
(15
1 0C) __

/\ ' ’
3 130 i 1 h
\ quantum _
O 0 0 O 0 semi-class.(1) .
» U U a U semi—class.(2) a
’ A A A A semi—class.(3) ‘

DIS —
0'0 TIL—fl 1* L—l_‘1_l_—l‘—l*I—*l—wl “I “I T

 

0 25 50 75
q (MeV/c)

Figure 6.3: Same as Figure 6.1, but using a thermal source of (0’s at a temperature of
150 MeV. In this case, since the width of the (1) is much less than the temperature, semi-
classical treatments work remarkably well.

can be understood by noting that the correlation’s width should be determined by the condi-
tion AEI = 1, where AE x mm. The width of the correlation function for the lower panels,
where pl .1. p2, is more complicated since it is more sensitive to the spatial movement of
the resonance while it decays. The correlation functions for the p in Figs. 6.1 and 6.2 are
even more complicated since they extend to large relative momenta where AE approaches

q. Given the complicated kinematics involved, it is not surprising that the result is sensitive

to the exact form for the semi-classical treatment.

In the case of the p, none of the three semi-classical prescriptions for I(m) provides a
consistently good approximation to the quantum result for all scenarios shown in Figs. 6.1

and 6.2. They all exhibit significant deviations from the quantum result for high relative

118

 

 

1.0 FT I I I I TTI ITI I I r I L
- hermal a) Source-

 

  
  
   

 

 

 

 

: T: 150 MeV :
0.5 :- q 11 (p1+p2)‘_
A
1— A -+
T— )- AA 4
I ,_ .1
A 0-0 1— :1— 1_ _1—1—1-1— _1—1—1 _1
31.0 TUITIIIIIIIIInfi—E
\ quantum -
C.) o a a O semi—01038.0) -
- a U D D semi—class.(2) -
. A A A A semi-c1ass.(3) ‘
0.5 b
' W
0-0 L1.1—1—1_-_1-_-1—'1—-1_--1_--1—n1-—1_-1—-1_1“I'd
O 50 100 150
q (MeV/C)

Figure 6.4: Same as Figure 6.3, except that (11’s have average momentum P = 800 MeV/c.

momentum of the pion pair, above 200 MeV/c. In the case of the (1), where only small
relative momenta, q < 200 MeV/c, are relevant, prescription (1) seems to best reproduce

the quantum result, as can be seen in Figs. 6.3 and 6.4.

It should be emphasized that the correlation from the nonzero extent of the resonance
source function has been factored out in this calculation. The deviations of the semi-
classical results occur for relative momenta of a hundred MeV/c or more. In a heavy-ion
collision the factor (exp[i(p1 — p2) - (xA — x3)]) in Eq. (6.43) would tend to zero for relative
momenta much greater than 50 MeV/c due to the large spatial sizes of the emitting regions.
Thus, the form of C’ becomes irrelevant for higher relative momenta unless the source sizes

are small, e.g., pp collisions.

These findings imply that the semi-classical treatments work quite well for the larger

119

sources considered with heavy-ion collisions. However, for smaller sources, especially
when the resonance widths are comparable to the temperature and resonances are produced
far off-shell, the semi-classical treatments significantly deviate from the quantum result.
This can be linked to the failure of the usual exponential decay law when off-shell reso-

nances are involved which will be investigated further in the following.

6.2.3 Limit of a Narrow Resonance

In this discussion, the source function is considered for a narrow Breit-Wigner resonance
A. Then, according to Eqs. (6.2) and (6.10), the probability of emitting a pion pair with
momenta k and p is

dN - . _ 2
(ZEk)(2Ep)m = ;lg/d‘ixd‘ixA GA(x—xA)TA,fgk(xA)e'(p+") x (6.47)

___. gZ /d4x,,d4x SA(p+k,xA)K(p+k,x—.IA). (6.48)

Here, K represents the probability of a resonance carrying momentum p + k propagating

from xA to x. It may be expressed in terms of the Wigner transform of the propagators,

(148 .
K (IN) = f (27% e'sq'xGMp + 54/2)GA (p — 811/2), (6.49)
which is in general a complicated function. In the limit of a narrow resonance, the propa-

gators can be expressed in a simplified form,

1 i
2MA p0 — 15,, + ir/2'

 

01(1)) = (6.50)

The narrow resonance limit is adopted here to illustrate the quantum nature of the propa-
gator. In addition, particles are assumed to be emitted with equal and opposite momenta.
Integrating Eq. (6.49) over spatial coordinates gives the probability of the resonance prop-
agating for a time t

_(EAEL)

2FM2 r sin
1115.: s———A/d3 K , =—e "I" , 6.51
M ,) n x (N) n (t)6 M ( )

120

where AE : Ek + Ep — MA is the off-shellness. When Eq. (6.51) is integrated over the

off-shellness, the expected exponential behavior is obtained.
[d(AE) 9mm) = re-F', (6.52)

whereas integrating over 1 describes the preference for emitting the particle with energy

close to on-shell.
1 172
11 A1?2 + (172)2'

 

[d1 mam) = (6.53)

The oscillating term sin(2AEt), which is responsible for preferentially emitting reso-
nances with small AE, also governs the distribution of emission times. Classical simu-
lations, which are typically based on Monte Carlo algorithms, cannot easily incorporate
regions with negative probabilities as suggested by this form. The mean propagation time
in a transport simulation could be altered to match the mean time of the quantum propaga-

tor,
_ fdt t 9((AE,1) _ 172

W— fdt R(AE,1) —AE2+(l“/2)2'

 

 

(6.54)

However, the second moment for the time would not correspond to that expected for an

exponential decay with the same average time.

31“2 — 4AE2

<22) — 1
T 801132 +r2/4)

2 74 2(1)? (655)

 

Not only is this form inconsistent with exponential decay, (12) can even become negative
for resonances far off shell. This illustrates, an a formal level, the need for performing the

quantum corrections described above.

121

Chapter 7

Balance Functions — A QGP Signal?

Balance functions were originally proposed to investigate hadronization in jets produced
in proton—proton collisions [59, 60]. By studying what Drijard et al. called the “associated
charge balance function” in a series of papers, they found that the leading hadrons carry
information about the quantum numbers of scattered partons [60]. They further concluded
that jet fragmentation is universal, i.e., that the “charge adjustment between the hadronic

fragments and the partons is independent of the latter’s quark content” [61].

Similar studies were performed with balance functions in the e+e" ——> 2 jets annihila-
tion processes. Leading particles in opposing jets were shown to be correlated across a
long range [38]. These long range correlations were then employed to demonstrate that the
primary quark of a jet ends up in the leading hadron [29, 8]. Short range correlations due to
local quantum number conservation were used to study jet fragmentation. KI? strangeness
correlations and heavy hadron correlations were found to be a cleaner signal than elec-
tric charge correlations because, respectively, much fewer strange quarks are produced [8]
and momentum distributions of heavier particles are less smeared [9]. These studies have

largely confirmed the fragmentation models.

122

 

7 .1 Balance Functions for Relativistic Heavy Ion Collisions

7.1.1 Theory

Bass, Danielewicz, and Pratt proposed balance functions as a signal that can distinguish
between early-stage and late-stage hadronization [24]. The balance functions are calculated
for a charge a and a balancing anti-charge b, which are created pairwise. They are defined

as

l

$(A21AI) = E{

 

 

<~<A31A1))— 00541)) (womb) — 005A?»
<N(A‘r’>> + 001)) 1 ’ ‘7'”

where A? signifies some condition involving particle type a, e.g. “11+ detected anywhere”,
N (A?) is the number of particles of type a that satisfy condition A; in one event, and
N (Ab , A?) is the number of particle pairs meeting both conditions simultaneously in a single

event. The other terms are defined accordingly.

The width of the balance function can be used as an indicator for the average phase-
space separation of charge—anti-charge pairs at freeze-out. This information can be used to
infer the time scale at which hadronization takes place in a relativistic heavy-ion collision
as much of the charge in a relativistic heavy-ion collision is created at hadronization. Early
hadronization, at I < l fm/c, [131, 100], would result in a broader balance function because
charge—anti-charge pairs that are created close to each other at an earlier time have a longer
time to diffuse apart by scattering off other particles. In addition, the velocity gradients in
the collision region are higher at early times, 1112/ dz = 1 / I, which contributes significantly
to the separation of the pair. A broad balance function indicating an early hadronization

would rule out the possibility that a QGP is formed in relativistic heavy-ion collisions.

Late-stage hadronization is consistent with the formation of a QGP that would exist
for some time during the collision before the partons are confined into hadrons. This late—

stage creation of charge—anti-charge pairs leaves them less time to separate from each other

123

with significantly smaller velocity gradients in the collision region. This scenario would,
thus, lead to a narrower balance function. However, a narrow balance function does not
necessarily constitute a unique signal of delayed hadronization as there might exist some

exotic mechanism that constrains charge—anti-charge pairs to stay close to one another.

Besides from the diffusive separation of charge—anti-charge pairs described above, the
balance function width receives another contribution from the charge—anti-charge creation
mechanism. An example for such a model is the widely-used Lund string fragmentation
model [15, 13]. As the string fragments in the model, quark—anti-quark and diquark-anti-
diquark pairs are produced and combined with other quarks or diquarks to form hadrons.
The balancing charges of the quarks, thus, end up in hadrons separated by a distance on the
order of 0.5 fm, when the hadrons are released from the string fragmentation. This initial
spatial separation of the charge-anti-charge pair can lead to large separations in phase space
in the final state because of the velocity gradients, which are particularly large during the

early stages of the collision.

When using a detector with perfect acceptance, the extracted balance function Eq. (7.1)
is normalized to unity, even if the multiplicity of the balancing charges a and b differ. Bass,
Danielewicz, and Pratt ensured the latter property in the balance function by introducing

the second term, which was missing from previous definitions:

1

/dArd/\1’B(A2,Ar)=§{

 

 

MbMa — Ma(Ma — 1) + MaMb — Mb(Mb — 1)

Ma Mb } z 1. (7.2)

The normalization is different, when extracting balance functions from experimental data
obtained with limited detector acceptance. This and other effects of limited detector accep-

tance on the balance function can be corrected for as discussed in [95].

124

7 .1.2 Experiment

The STAR experiment at RHIC was the first to determine a balance function for a rela-
tivistic heavy-ion collision [27]. The electric charge balance function as well as the n+n‘
balance function were studied. Condition A1 of Eq. (7.1) was taken to indicate “anywhere
in the detector”, while A2 signified “separated from the first particle by a pseudo-rapidity

gap of An”, where the pseudo-rapidity is defined as

1 1p 1 ‘1' pz
11 = —log———. (7.3)
2 ”’1 _' [72
In this case, the balance function can be written simply as a function of relative pseudo-

rapidity only,

93(An)=1{

 

 

(7.4)

(Nbu(AT1)>—' (Naa(An)) + (Nab(An)) _ (Nbb(Ay)>}
2 .

(Na) (Nb)

The STAR collaboration found that the the balance function narrows from peripheral
to central collisions [28, 136]. This finding is consistent with the notion that the dynam-
ics in a relativistic heavy-ion collision differs significantly between central and peripheral
collisions. The narrow balance function in central collisions and wider balance function
in peripheral collisions are consistent with the notion that the production of the QGP is
more likely in the former than in the latter. To further substantiate this hypothesis, the ex-
perimental balance function width is compared to that extracted from the Hij ing model
[151, 86]. Hijing is essentially a Glauber-type model but also incorporates soft, non-
perturbative QCD phenomenology. As a model, Hij ing does not contain the production
of a QGP. The balance functions are calculated from the Hij ing model simulations and
compared to the experimental data without further modification. The extracted balance
function widths are essentially flat as a function of impact parameter as the dynamics of the
superposition of pp collisions does not change from central to peripheral collisions. For

peripheral collisions, the experimental data were matched by the model calculations within

125

 

 

I All Charged Fairs |

 

 

 

 

 

 

 

 

 

 

 

 

 

008 I
_
. I :5) - .. 2. 53"“ it'- 31
0.7%. ..... r "3,”; .7. ............ r “‘52:? ....................... ‘ ................................................ - :3.
.c "' . . : :
g : : : : : : :
g 0.7 [__u_*..- ...... . .................... . .................... . ................... ...... .......
g _ +
‘g’ 055 __ ............... . .......... 1 ......... i ..... i .............. i .................... i ........................................ i .................... i...
n- F i ' + ;
a . r
5 0.5 _. ............... .................... .................... ........ “F“ ........................ #- .................
g + Amway];NN = 200 GeV
_g. HIJING, Au+Au. W: 200 GeV - g
0.55 ' p+p' G = 200 GeV ...................................... i-
: = * 2 2t
STAR Pfeliminary . . i
o 5 1 LLI I I I l I I I I ILL l L L I I L J l I 1 L I l I l I I l l
' 0 50 100 150 200 250 300 350
"part

Figure 7.1: Width of the electric charge balance function for gold on gold collisions at
W = 200 GeV as a function of the number of participants [136]. A peripheral collision
involves a small number of participants, whereas a more central collision involves more
participants. The narrowing of the balance function with increasing centrality is consistent
with an increasingly delayed hadronization.

errors. For central collisions, however, the experimentally determined balance function was

significantly narrower than the balance function from H ij ing.

The result described above indicates that peripheral relativistic heavy-ion collisions
follow the dynamics of well-known pp reactions, whereas central collisions of heavy ions
create some new phenomena. Yet, there needs to be more certainty in the model calcula-
tions to rule out any source, other than late-stage hadronization, for the narrowing of the
balance function in central collisions. Such studies are essential to firmly establish signals
of the QGP. One issue, to be investigated in the following, is the question how the lack of

hadronic rescattering in the Hij ing model affects the extracted balance function.

126

7 .2 Model Calculations with Hadronic Rescattering

The Hi j ing model, which was used to extract a balance function as theoretical comparison
to the experimental balance function, does not include hadronic rescattering and hence
does not accurately describe the hadronic evolution of a relativistic heavy-ion collision.
The difference noted between the model calculations and the experimentally determined
balance function may conceivably be attributed to this deficiency in the model. To move
further toward establishing the existence of a QGP, balance functions should be calculated
from a simulation, in which Hij ing is complemented with a simple model for hadronic
rescattering. In practice, an event is created first with Hij ing and then fed into a Gromit
model with a set of particle types corresponding to Hi j ing’s and with elastic and resonance
interactions between those particles. 24,400 simulated events are used to calculate the

n+n' balance function in relative rapidity and in
2
qinv = (P1 - P2)2 (7.5)
as well as the electric charge balance function in relative pseudo—rapidity.

With perfect detector acceptance, the n+n" balance function in relative rapidity turns
out to be broader with hadronic rescattering than without, as can be seen in Figure 7.2. The

statistical errors in Figure 7.2, which are smaller than the symbols, are calculated according

 

 

 

 

10
l (MA/A)» +<Naa(A1’)> (Nab(A)> ‘1' (th(A)>
623 Ay : ’ + . (7.6)
1 ) 2m\/ (W (W
The width of the balance function, shown in Table 7.1, can be quantified With the aid of
XiAYIBi

A)‘ = --——-. (7.7)
< > 2131'

where Ayi is the mid-value in the rapidity difference bin. The statistical error of this balance

function width is obtained with

_ 2113x2531; 21 13.115312 215312
MA» — (M) \/(ZiA)’iBi)2 — 2(2131)(214yi31) (213112.

 

 

 

(7.8)

127

 

0.7 T I i i j

 

 

 

 

Hijingr T a
06 10 0° Hljing+Gromit in q
0 O
O
0.5 L 00
1590mm 00
I; 0.4 ' mum
s “1380
m 0.3 ~ 5111111
0013
0.2 h 03% -1
O In III
0.1 ~ 0° “’08” '3

 

 

Figure 7.2: Pion balance function in relative rapidity from Hij ing and Hij ing + Gromit
without acceptance cuts.

 

 

 

 

 

 

[ Variable | Model 1 Width 1
. . Hijing 0.893 10.016
K (it
“I” 1y Hijing+Gromit 1.160i0.016
2 Hij ing 0.718 :1: 0.013 Gev2
‘1'") Hij ing+Gromit 0.646 :1: 0.012 GevT

 

 

 

 

 

Table 7.1: Summary of balance function widths without acceptance cuts.

Results for the balance function in qiznv are shown in Figure 7.3 and Table 7 .1. They
indicate that the width of the balance function decreases by a small amount due the inclu-
sion of hadronic rescattering. This narrowing can be attributed to the cooling induced by
the hadronic rescattering. It should be noted that with infinite hadronic rescattering, i.e.,
vanishing mean free path for hadrons, the width of the balance function in qizm, would be
even narrower. In that case, the diffusive contribution to the width of the balance function

would vanish, leaving only the contribution due to the particle creation mechanism.

The pronounced peak in the balance function around ‘1an z 0.52 GeV in Figure 7.3 is a

manifestations of a stron n+1t‘ correlation due to a resonance deca . The sition of the
g Y P0

128

 

_L
O
_(
.1

 

 

 

 

 

 

Hijiné L a
9 r in Hijiry+Gromit 1:1 ‘
8 m "
,1: 1
o in
... 5 ' Go ‘
(NO-E 5 _ a .
5 '° 9m
4 P am 000 q
. a
3 029° 1
2 ~ Dag-n00 .
in a
I“ I
1 ~ I: BIB ‘
II
o . . . L {uh'Pfl'Ur-glsl

 

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
qlznv [GeVz]

Figure 7.3: Pion balance function in qiznv from Hij ing and Hij ing + Gromit without ac-
ceptance cuts.

peak corresponds to an invariant mass of around 773 MeV, close to the pole mass of the p
mesons. Hence, the correlations in the peak arise from the decay of the p0 meson into a

n+1t‘ pair.

In summary, it is found that the balance function as a function of relative rapidity is
widened, when hadronic rescattering is included, but slightly narrowed as a function of
(112w These seemingly contradictory results can be explained as a result of the cooling of

the system. For small relative rapidity Ay
(1.2... w m?(Ay)2, (7.9)

where the transverse mass is defined by

111,: MEL—pg. (7.10)

Hadronic interactions cool the particles as the system performs work on the longitudinal
expansion of the collision region, thus leading to a lower transverse mass. If the average

separation of the particle—anti-particle pairs in q?nv decreases less than the transverse mass,

129

then, by Eq. (7.9), the separation in rapidity has to increase. This explains the different
findings for the balance function widths and implies that the balance function in qizm, is the

cleaner observable.

7.2.1 Accounting for Experimental Acceptance

The discussion in the previous subsection assumed a perfect detector, neglecting the un-
avoidable limited acceptance of real experiments. The effects of a limited acceptance are
studied here by passing particles through a simulated acceptance filter of the STAR de-
tector. One of the main features of this acceptance filter is to reject particles with large
pseudo-rapidity It” _>_ 1.3 and/or high transverse momenta p, = W Z 700 MeV.
Only particles that pass this filter will then be used in the calculation of the balance func-
tions. These calculations start form the same 24,400 simulated events as above in order to

arrive at comparable results.

 

0.45 r . f

r

 

 

 

 

 

 

 

Hiiing I o
0.4 L? g HllinyGromlt m ‘
0.35 - ‘l‘ e,
0.. . l 3*
0.25 - ,.
9 ell?
g 0.2 - ,. 53 - «
m 0.15 ~ I? " 9%, l
0.1 L {>93
0.05 " Ea
0» aaaufllaaaaa
O 0.5 1 1.5 2 2.5
AV

Figure 7.4: Pion balance function in relative rapidity from Hij ing and Hij ing + Gromit
with STAR acceptance cuts.

When the STAR acceptance cuts are applied, the n+1t‘ balance function in relative

130

 

rapidity, shown in Figure 7.4 and Table 7.2), displays nearly the same width, with or without
hadronic interactions. This indicates that the same pions that are rejected by the acceptance
filter, are also linked to the widening of the balance function in relative rapidity, when a
perfect acceptance is assumed. The rejected pions are characterized by high p, , and/or large
pseudo-rapidity and are most affected by the cooling due to hadronic rescattering. This
underscores the hypothesis stated in the previous subsection that cooling is responsible for

the widening of the balance, when hadronic interactions are included.

 

 

 

 

 

 

 

[ Variable LModel 1 Width I
. . Hij ing 0.541 :l: 0.009
R dt
3‘” 1’ Hijing+Gromit 0.550 t 0.009
. . Hij ing 0.628 d: 0.010
P d - dt

56" ° rap‘ 1 y Hi jing+Gromit 0.675 i 0.010

2 Hij ing 0.262 :l: 0.005 GeV2

qmv Hi j ing+Gromit 0.163 :1: 0.005 GeV

 

 

 

 

 

 

Table 7.2: Summary of balance function widths with STAR acceptance cuts.

 

 

Hijing I o
Hijing+Gromit In

 

 

 

L

)

2
lnv

B(q
00

 

 

 

Figure 7.5: Pion balance function in qfiw from Hi j ing and Hij ing + Gromit with STAR
acceptance cuts.

The balance function in qiznv for pions filtered by the STAR acceptance confirms the

131

previous finding that the width of the balance function is not affected by hadronic interac-

tions, as can be seen in Figure 7.5 and Table 7.2. The same is found again when calculating

 

0.4 . . .

 

 

 

 

 

 

 

Hijing L a?
0.35 g illngGromlt El «
aa-éiQré
0.25 ~ I}? g, . M
‘ a:
E ”'2' ‘96“
E 0.15 . ‘ €1$fi .
01- 03$?
0.05 r 4’ Q E i 3 -
0» '3uuaam
0 0.5 T 11.5 2 21.5

An

Figure 7.6: Pion balance function in relative pseudo-rapidity from Hij ing and Hij ing +-
Gromit with STAR acceptance cuts.

the electric charge balance function as shown in Figure 7.6.

7 .3 Corrections for Net Charge

The definition of the balance function Eq. (7.1) assumes that all charges are produced
in pairs, i.e. there is no net charge in the system. Nevertheless, due to the presence of
unbalanced protons and neutrons in the colliding nuclei in a relativistic heavy-ion colli-
sion, there are net charges like electric charge and baryon number. Much of the quan-
tum numbers of the colliding nuclei in a relativistic heavy-ion collision at RHIC energies
are carried by leading baryons. For instance, the ratio of the proton yield to the anti-
proton yield in the mid-rapidity region was measured by the STAR collaboration to be
0.65 i 0.01 (stat) i 0.07(Syst) at RHIC energies [6], indicating that few initial protons are left

in the center of the collision region. In a typical experiment with relativistic heavy-ion col-

132

lision, the detector acceptance does not cover the leading baryons so that those net charges
do not affect the balance functions. This can be seen in Figure 7.6, which displays the

electric charge balance function that was calculated using the STAR acceptance filter.

Any correlation among the excess charges would dominate the balance function if
Eq. (7.1) were used to calculate the balance function with a perfect detector acceptance.

This effect is illustrated in Figure 7.7, where the electric charge balance function is shown

 

0.5 I r F r r t

 

 

 

Hfiifng L fl 0
0.4 ~ lflnQGromlt In ~

     

 

0.3 F

0.2 P,,,_ 7

0.1

8(A n)

0 ..

-O.1 ~
-0.2 ~
-0.3

 

 

>-
t.
...
L.

l l l

0 2 4 6 8 10 12 14 16 18 20

-0.4

 

An

Figure 7.7: Uncorrected balance function in relative pseudo-rapidity from Hi jing and
Hi j ing + Gromit for all electric charges.

to be heavily dominated by the pp" contribution. These correlation can be seen in isolation
in the p13 balance function shown in Figure 7.8. The p13 balance function exhibit strong
correlations at AT] :- 0 and at A1] : 10 as the protons from the incident nuclei are mainly
concentrated around 1] 2 21:5 at RHIC energies. A balance function such as shown in Figure
7.7 does not reveal the pair correlation of interest and is little use as a signal of late-stage

hadronization.

The complication of net charges as pointed out above can amount to a serious imped-
iment for the use of balance functions, even if the acceptance cuts of real detectors effec-

tively removed the leading baryons which carry most of the net charge. After all, there

133

 

 

.-°~"‘--... fiijingJ I -—
' " HijingGromit -------- «

 

 

 

 

 

3(A n)

 

 

.5 ' 1 L L r r 1 r r r

0 2 4 6 8 10 12 14 16 18 20
An

 

Figure 7.8: Uncorrected balance function in relative pseudo-rapidity from Hijlng and
Hij ing + Gromit for protons and anti-protons.

still are contributions from net charges in the mid-rapidity region because it is not entirely
baryon free, i.e. there are more baryons than anti-baryons, as evident from the ratio of
the proton yield to the anti-proton yield. Furthermore, net charges do not only exist for
Coulomb charge or baryon number, there is also a small excess of negative pions over pos-
itive pions. This discrepancy arises as the initial nuclei contain more neutrons than protons
which leads to a negative net isospin. In order to conserve isospin more negative pions
are produced than positive pions. Therefore, removing the contributions of excess charges

would significantly improve the balance functions as a signal.

To derive such a corrected balance function, it is assumed that there are two independent
sources, one source that creates a certain particle type b “unpaired”, i.e., individually and
independently from other particles, and another source that creates particles in charge—
anti-charge pairs ab. It is also supposed that the particles in the final state can be associated

uniquely with either source. In this case, the particle numbers from Eq. (7.1) can be written

134

in term of separate contributions from each source.

Na = N5, (7.11)

Nb : N}: +N“, (7.12)

where N}; is the number of particles from the unpaired source and N5 = N5 the number
of particles that were created in pairs. The pair counting terms in Eq. (7.1) count pairs
irrespective of their origin. However, for the purposed of this discussion, the counted pairs
need to be separated into classes based on the original source of the particles. In short

notation

N... 2 N55, (7.13)
Nab = N5: +Nfg, (7.14)
A7,, = N5: + N55, (7.15)
Nbb = N5: + N55 + N5; + Ngg, (7.16)

where NZ: refers to the number of particle pairs ba, where b originated from the unpaired
source and a from the paired source and the other terms are defined accordingly. The goal
is to define a balance function that only contains contributions from paired terms, which is,
therefore, independent of net charges. In addition, the corrected balance function should
yield the same result as in Eq. (7.1) in the case where the unpaired source does not produce

any particles.

First, the two different denominators in Eq. (7.1) are replaced by a single denomina-
tor Nf to remove the particle number from the unpaired source. The combination of the
four pair counting terms in Eq. (7.1) should be preserved as closely as possible, since
they remove other event correlations like anisotropic flow and jets. Instead of altering the
pair counting term, the unwanted terms are subtracted from the numerator (N (Ab ,A‘l1 )) —

(N(A",A‘l’)) + (N(A‘2’,A’,’)) — (N(Ab,A’l’)). The contribution of N5: equals N5: because

135

paired particles ab are created with the same probability distribution, regardless whether
they are of type a or b. The two terms cancel in the balance function. The same holds for
NZ: and NEE, which leaves N3}; as the sole term that needs to be subtracted. Symbolically,

NE}; can be written as

NEE = (N3) ®(NZ) (7.17)
i 0%.
r j z

to signify that unpaired particles are checked against other unpaired particles when building
the histogram for the number of pairs. The expression j 75 i is a reminder that a particle
must not be correlated with itself. Since the emission of unpaired particles is assumed to be
independent, the second set of particles may be taken from a different event with the same

multiplicity of particle types a and b

N3}; = (N8) ® (N3) (7.18)
V
from one event, 1' from another event, j¢r

Practically, it is impossible to discern which source, paired or unpaired, a particle originated
from. The distribution of unpaired particles b is equal to the distribution of all particles b

minus that of paired particles b or a. Therefore

N5}; = whom.) a (NbeNa) . (7.19)
\—-\/—" \-\/—’
from one event, 1' from another event, jyér'

The event mixing is necessary to avoid reintroducing the pair correlations between particles
from the same event. The expression from Eq. (7.19) can be rewritten in term of particle

pairs from mixed events
N33 ——> <N’(A’2’1A’.’)>— <N”(A”.A‘.’)>+<N”( 5A?» — <N’(A",Al’)> (7.20)

where N’ (Ab, Alf) counts the number of pairs from mixed events that satisfy the two condi-

tions simultaneously, N” ( g, A?) means the same, except that particles with the same index

136

 

are not counted. Finally, the corrected balance function is given by

1 { <N<A3,Ar>> — <~<A5Ar>> + <~< 5A?» - <~<A5A€>>

 

B<A2IA1)= 5 WW,» (7.21)

 

_ <N"(Ab,Ar)> — <~"<A5A1>>+<N'<A5A’r>> — <~'(A3,At>> }
<N<Arr>

To test this new expression, the corrected balance function is calculated for the worst

case of pp" balance functions, where the net charge plays the most important role. The re-

 

o.5 I I l I f

 

Hiring' ‘ _
HflinyGromlt .......

 

 

 

0.4 ‘

0.3

3(A n)

0.2

0.1 -

 

 

 

 

l l l l l l l l l

02468101214161820
An

 

Figure 7.9: Corrected pp" balance function in relative pseudo-rapidity from Hij ing and
Hij ing + Gromit.

sults in Figure 7.9 show a dramatic difference to the balance functions calculated with the
original definition in Eq. (7.1), which were shown in 7.8. The effects of proton concentra-

tions at n = 35 have been removed completely.

With this result, the corrected form for the balance function can now be used to calculate
the electric charge balance function with a perfect detector acceptance, as illustrated in
Figure 7.10. The balance functions are evidently wider under the influence of hadronic
rescattering. This is consistent with the similar finding from section 7.2 about the n+n‘

balance function in relative rapidity.

137

0.55 T . . f

00 Hiiing o
0.5 "a o o Hing-imam“ n
0.45 . ° ‘
0.4 ~ ‘

0.35 can u .

 

 

 

 

 

 

 

 

E 0.3- o oo ‘
< D
E 0.25“ ODD "
Dr:

02» 0000 -

ODD
0.15- o n0 .
0.1~ 0 09 an.

O
0.05“ ca"
0 4 r r r 1
0 0.5 1 1.5 2 25

An

Figure 7.10: Corrected electric charge balance function in relative pseudo-rapidity from
Hij ing and Hij ing + Gromit.

Since in relativistic heavy-ion collisions, the multiplicity of 1t“ is greater than 10*, the
influence of the net pion charge on n+1t‘ balance function needs to be investigated. The
calculations from Section 7.2 are repeated with the expression for the corrected balance
function and shown in Figures 7.11 and 7.12. Little change is found in the shape of the
n+7t‘ balance function. The calculated widths of the corrected balance functions, (Ay) =
1.121 for Hij ing and (Ay) : 1.379 for Hijing + Gromit, are larger than the widths of
the balance functions in the original definition, Table 7.1. However, the absolute width is
of little interest in this study. The main focus is the comparison of the relative balance
function width with hadronic rescattering to the case without. Even with the corrected

balance function, the width is still larger with hadronic scattering than without.

138

 

0.7 T I r T

Hijing, corr. —

0., _ [ileng .

 

 

 

 

8(A y)

 

 

 

 

Figure 7.11: Comparison between corrected and original 71+TI_ balance function in relative
rapidity from Hij ing.

0.5 . . .
0.45 ix
0.4 -
0.35 L
0.3 .
0.25 ~
0.2 ~
0.15 L
0.1 ~
0.05 -

 

 

Hijing+Gromit, corr. —
ijing+Gromlt ----- ‘

 

 

 

B(A v)

 

 

 

 

Figure 7.12: Comparison between corrected and original n+7t‘ balance function in relative
rapidity from Hij ing + Gromit.

 

Chapter 8

Conclusions

Experiments with relativistic heavy-ion collisions offer the possibility of studying some
predicted, but difficult to observe, QCD phenomena, two of which were discussed in this
thesis: the phase transition from hadronic matter to the quark-gluon plasma (QGP) and the
restoration of chiral symmetry. Neither effect is directly accessible by experiments due to
the small scale of the collision region in both space and time. Hence, the existence of the
phase transitions must be established from the information contained in the collision debris
recorded by the detectors. The complexity of both the underlying theory of QCD and the
dynamics in relativistic heavy-ion collisions make theoretical treatments from first princi-
ples intractable. This difficulty is further compounded by the distortion of any effect, which
is predicted to result from the phase transitions, in the collision dynamics. Consequently,
much theoretical insight on relativistic heavy-ion collision is deduced from analytical as
well as Monte Carlo models. In this thesis, new models were contributed for the study of

relativistic heavy-ion collisions and phenomena associated with them.

In Chapter 2, alternative sources for anomalous isospin fluctuations other than disori-
ented chiral condensates (DCC) were studied by considering the influence of total isospin
conservation, quantum symmetrization, and resonance decays. For this purpose, novel ex-

pressions were obtained for the multiplicity distributions and the isospin fluctuations for

140

 

a canonical ensemble, in which total isospin as well as additive quantum numbers are ex-
actly conserved. The formalism was extended to include both pions and resonances and
can account for Bose-Einstein symmetrization of the pion wavefunction. Numerical cal-
culations were then performed to study the effects of total isospin conservation, quantum
symmetrization, and resonance decays on the width of the multiplicity distribution, which
can be squared to obtain the isospin fluctuations. The widths of the multiplicity distribu-
tions are found to be largely dominated by the behavior of the tails, thus making it im-
perative to perform exact calculations. Such calculations were made possible by recursion
relations that circumvent the summation over the immense number of partitions in the par-

tition functions.

It was found that conservation of total isospin and its projection has little effect on the
width of the multiplicity distributions, when the systems are larger than a dozen particles.
At high phase-space densities, inclusion of Bose-Einstein symmetrization leads to a multi-
plicity distribution that is much broader than a random distribution. However, addition of
resonances more than compensates for this broadening and narrows the multiplicity distri-
bution below the width of the random distribution. Both effects are small when the phase

space density is below 0.1 fm'3.

Another indicator for the restoration of chiral symmetry was studied in Section 4.3.
Resonance mass distributions are believed to broaden and shift toward lower invariant mass
in an environment with restored chiral symmetry. The p meson with its relatively short
lifetime would be a good probe for such modifications to the resonance properties, if the
decaying p mesons from the initial state are replenished through formations from two pions.
A model based on the Gromit framework was used to gain insight into such a regeneration
of p mesons in the collision dynamics. It was found that the p is sufficiently regenerated
to supply probes at all stages of the collision. Furthermore, the regeneration of p mesons

is sensitive to the existence of hadrons at the early, hot and dense phase of a relativistic

141

heavy-ion collision. A larger number of p mesons reconstructed from lepton pairs from
a relativistic heavy-ion collision than from a cocktail of nucleon-nucleon collisions would
indicate the existence of hadrons at early stages, which is inconsistent with the formation

of a QGP.

The Gromit framework for microscopic transport models based on the Boltzmann trans-
port equation was introduced in Chapter 3. It provides a general framework that can be used
for any degrees of freedom and any interaction; it is expandable, modular, transparent, and
easy to debug and maintain. Gromit, therefore, provides a common framework that can be
used to study the different stages of a relativistic heavy-ion collision at RHIC energies with

their vastly different dynamics.

Like some other Boltzmann transport models the Gromit framework introduces an arbi-
trary model parameter, the sampling factor A, which has a profound effect on the extracted
results. These effects were investigated in Chapter 4 for simplified models and shown for
particle spectra and elliptic flow, to vanish in the Boltzmann limit, A. —> 00. Large sampling
factors were also found to eliminate sensitivities to a variety of arbitrary choices inherent
to Boltzmann-type simulations and to solve a variety of problems related to acausal prop-
agation. This is linked to the reduction of the interaction range between particles in the

simulation that restores the locality of interactions.

In nature, nevertheless, particles indeed interact over a finite range, either by exchange
of off-shell particles or through a mutual interaction through classical fields. Therefore,
Chapter 5 sought to understand the degree to which these effects are physical as opposed
to representing numeric artifacts. Nonlocal effects were interpreted in terms of viscous
parameters, allowing the incorporation of nonlocal effects into hydrodynamic models in a
straight-forward manner. Knowing the viscous parameters also provides criteria for tuning

Boltzmann algorithms so that they are consistent with quantum transport considerations.

142

The tuning can be accomplished by either changing the sampling factor or by adjusting
the scattering algorithms. The simple manner in which the viscous parameters scale with

density and cross section should simplify such a procedure.

As transport theories address the first one fm/c of a relativistic heavy-ion collision,
the role of nonlocal interactions becomes increasingly important. For times above 2 fm/c,
it is unlikely that the nonlocalities play any significant role as the effects scale as 1’3.
Since the nuclei pass one another at RHIC in less than 0.2 fm/c, whereas cross sections
approach a square fm, nonlocal effects might provide a nonnegligible source of stopping as
the viscous drag converts longitudinal collective velocity to heat. The role for such effects

in the stopping phase at LHC collisions should be even greater.

The influence of resonance decays on correlation function were studied in Chapter 6. A
new method was proposed to correctly calculate correlation functions from semi-classical
models, exploiting the source function of the resonance and its creation point in spacetime
as opposed to the creation points of the final-state pions. This method is easily applicable to
generate correlation functions from the event histories of simulations. When using direct
pions, the creation points of the final-state pions provide all the necessary information
for creating correlation functions. By considering the creation points of the resonances
that decayed into the final-state pions, all spacetime points at which the resonance might
have decayed were coherently accounted for by modifying the prescription for generating

correlation weights.

It was found that the proposed modification is important only for sources that are
quantum-mechanical in nature. If the source is large and the product of the momentum
and spatial uncertainties are large, Ap Ax >> Ii, the behavior of the correlation function is
dominated by the exponential term that is determined by the points at which the resonance

is created. Quantum considerations in resonant decays could play an important role when

143

considering the decay of small sources that push the limits of the uncertainty principle.
However, such sources are also accompanied by questions regarding the quantum nature of
the source functions responsible for initial creation of the resonances, i.e., the smoothness
approximation might not be justified. For such problems, unless the off-shell behavior of
the source functions is known, as in the case of a thermal model, the treatments presented

here address only half the problem.

In Chapter 7, the Gromit framework was also used to explore the influence of hadronic
interactions on balance functions, a new observable that has been proposed as a possible
signal of late-stage hadronization. It was found that hadronic interactions widen the n+n’
balance function in relative rapidity, but not in qiznv. These seemingly contradictory results
can be explained as a kinematic consequence of a cooling of the collision region due to

hadronic rescattering.

Furthermore, it was demonstrated that previous definitions of balance functions are
applicable only for a vanishing net charge. A new formalism for constructing balance
functions in the presence of net charges was introduced and studied numerically. Though
not entirely efficient, the new balance function was shown to remove even the dramatic
effects of the net electric charge of the protons in a relativistic heavy-ion collision. Finally,
for n+7t‘ the width, not so much the shape, of the balance functions were found to be
somewhat affected by the new formalism. This can have implications for other studies but

does not alter the conclusions about the effects of hadronic rescattering mentioned above.

144

 

Appendix

145

 

Appendix A

Parameter Definitions for Numerical
Simulations

A.1 Initialization File for Pionwind Study

# specifications for pionwind calculation
particles {
add Species N { mass
add Species pi { mass

0.938; maxXS= 10; pngdz 2212; baryon: 1;}
0.1396; maxXS= 10; pngd= 211; }

}

collisions {
add PionWindModel {
xs = 10;

# collisions 2 pi pi, pi N; # comment out to disable N N
}

}

run {

t0 = l; tf = 64;

rand = R250 { }

events = 4;

testParticles = 32;

initCond = BjorkenInitCond {
T = { pi20.165, N=O.l80 }
N = { pi=2400, N=24O }

posDist 2 Radial { etamin = —2; etamax = 2; rmax = 5; }
}
generator: TauEtaGenerator {

boundary_cond= simple;

mesh_min= {—10,—lO,—4};

mesh_max= {lO,lO,4}

collision_file= col out;
}

146

’ ,‘o 12.-D S- .Q‘ ‘1‘)

 

 

 

analysis {

add Histogram {
variable = mt - m; weight = 1/(2*PI*mt);
species = pi;
xrange = [0, 2]; xstep = 0.05; format 2 "%.3f %10.4f %10.4f";
acceptance = Interval { range = [-0.5, 0.5]; variable = y; }
file = "pi.mt";

}

add Histogram {
variable 2 mt - m; weight = 1/(2*PI*mt);
species = N;
xrange = [0, 2]; xstep = 0.05; format = "%.3f %10.4f %10.4f";

acceptance = Interval { range = [-0.5, 0.5]; variable = y; }
file = "N.mt";

 

147

A.2 Initialization File for Balance Function Study

loadlib ""/1ib/libhijing.so";
particles {
add PDG {
names 2 gamma, D, D_s, pi, K, eta, rho, omega, K*, eta’,
f_0, a_O, phi, h_l, b_l, a_1, f_2,
N, Lambda, Sigma, Xi, Delta, Sigma*, Xi*, Omega
}
}
decays {
load "tables/resonances tbl";
l
collisions {
add ResonanceDecayInverter { } # 76 resonance formation channels

add AddQuarkModel { } # 2775 elastic channels
1
run {

t0=1;

tf = 1000;

rand 2 R250 { }
events 2 100;
initCond = HijingInitCond {
bmin = 0; bmax = 0; # min and max impact parameter
targ A; targA = 197; targZ = 79; # target (Au)
proj A; projA = 197; pron 2 79; # projectile (Au)
frame = CMS; # frame of calculation, do not change
eframe = 200; # cms energy
particlesToPropagate = hadrons; # propagate hadrons or partons
hadronDecay = 1; decay hadrons via jetset
jetQuenching : 0; jet quenching, turn off
nuclearShadowing = 1; nuclear shadowing is on
softRad = 1; soft radiation is on
baryonJunction = 0; baryon junctions are temporarily off

###=fl=#

}

generator: TauEtaGenerator {

boundary_cond= simple;
min_cell_width= {.625, .625, .625};
mesh_min= {-10,-10,—6};
mesh_maxz {10,10,6}

}

}

analysis {
add CustomOutput { when: initial; file: initial.out; }
add CustomOutput { when: final; file: final.out; }

148

 

 

Bibliography

149

Bibliography

[1] S. Abachi et al. Observation of the top quark. Phys. Rev. Lett., 74:2632-2637, 1995.

[2] F. Abe et al. Observation of top quark production in p’p collisions with the collider
detector at ferrnilab. Phys. Rev. Lett., 74:2626—2631, 1995.

[3] P. Abreu et al. Two-dimensional analysis of the Bose-Einstein correlations in e+e‘
annihilation at the 20 peak. Phys. Lett., B4712460—470, 2000.

[4] K. H. Ackermann et al. Elliptic flow in Au + Au collisions at 31010 = 130 GeV.
Phys. Rev. Lett., 86:402—407, 2001.

[5] C. Adler et al. Elliptic flow from two- and four-particle correlations in Au + Au
collisions at ,/sNN : 130 GeV. [nucl—ex/0206001], 2002.

[6] C. Adler et al. Mid—rapidity anti-proton to proton ratio from Au + Au collisions at
,/sNN = 130 GeV. Phys. Rev. Lett., 86:4778—4782, 2001.

[7] J. Aichelin. Can one extract source radii from transport theories? Nucl. Phys,
A617z510, 1997.

[8] H. Aihara et al. Observation of strangeness correlations in e+e‘ annihilation at
\/E = 29 GeV. Phys. Rev. Lett., 5322199, 1984.

[9] H. Aihara et al. Study of baryon correlations in e+e“ annihilation at 29 GeV. Phys.
Rev. Lett., 57:3140, 1986.

[10] N. S. Amelin et al. Collectivity, energy density and baryon density in pb on pb
collisions. Phys. Lett., B26lz352—356, 1991.

[11] A. A Amsden, G. F. Bertsch, F. H. Harlow, and J. R. Nix. Relativistic hydrodynamic
theory of heavy ion collisions. Phys. Rev. Lett., 35:905—908, 1975.

[12] Anthony A. Amsden, Francis H. Harlow, and J. Rayford Nix. Relativistic nuclear
fluid dynamics. Phys. Rev., C15z2059—2071, 1977.

[13] B0 Andersson, G. Gustafson, G. Ingelman, and T. Sjostrand. Parton fragmentation
and string dynamics. Phys. Rept., 97:31, 1983.

[14] B0 Andersson, G. Gustafson, and B. Nilsson-Almqvist. A model for low p(t)

hadronic reactions, with generalizations to hadron - nucleus and nucleus-nucleus
collisions. Nucl. Phys, B281z289, 1987.

150

[15] Bo Andersson, G. Gustafson, and T. Sjostrand. A model for baryon production in
quark and gluon jets. Nucl. Phys, B197z45, 1982.

[16] A. A. Anselm. Classical states of the chiral field and nuclear collisions at very high-
energy. Phys. Lett., B2171169—172, 1989.

[17] A. A. Anselm and M. G. Ryskin. Production of classical pion field in heavy ion
high-energy collisions. Phys. Lett., 8266:482-484, 1991.

[18] N. Armesto, C. Pajares, and D. Sousa. Analysis of the first RHIC results in the string
fusion model. Phys. Lett., B527z92—98, 2002.

[19] P. Avery et al. Bose-Einstein correlations in e+e‘ annihilations in the Upsilon re-
gion. Phys. Rev., D32:2294——2302, 1985.

[20] A. M. Baldin, L. A. Didenko, V. G. Grishin, A. A. Kuznetsov, and Z. V. Metreveli.
Universality of hadron jets in soft and hard particle interactions at hi gh-energies. Z.
Phys, C33:363, 1987.

[21] J. Barrette et al. Two-proton correlations from 14.6-A-GeV/c Si + Pb and 11.5-A-
GeV/c Au + Au central collisions. Phys. Rev., C60z054905, 1999.

[22] S. A. Bass and A. Dumitru. Dynamics of hot bulk QCD matter: From the quark-
gluon plasma to hadronic freeze-out. Phys. Rev., C612064909, 2000.

[23] S. A. Bass et al. Microscopic models for ultrarelativistic heavy ion collisions. Prog.
Part. Nucl. Phys, 41:225-370, 1998.

[24] Steffen A. Bass, Pawel Danielewicz, and Scott Pratt. Clocking hadronization in
relativistic heavy ion collisions with balance functions. Phys. Rev. Lett., 85:2689—
2692, 2000.

[25] E. Bauge et al. Observation of a saturation in the time scale for multifragment
emission in symmetric heavy-ion collisions. Phys. Rev. Lett., 70:3705-3708, 1993.

[26] L. Beaulieu et al. Signals for a transition from surface to bulk emission in thermal
multifragmentation. Phys. Rev. Lett., 84:5971—5974, 2000.

[27] M. Belt Tonjes, G.D. Westfall, and AM. Vandemiolen. Clocking hadronization
at RHIC. In G.D. Westfall and W. Bauer, editors, Proceedings of the 17th Winter
Workshop on Nuclear Dynamics. EP Systema Bt., Debrecen, Hungary, 2001.

[28] Marguerite Belt Tonjes. Using the Balance Function to Search for Late Hadroniza-
tion in Au+Au Collisions at a Center of Mass Energy of 130 GeV per Nucleon Pair.
Ph.D. thesis, Michigan State University, 2002.

[29] Christoph Berger et al. A measurement of charge properties of quark jets at PETRA.
Nucl. Phys, B214:189, 1983.

[30] G. F. Bertsch, P. Danielewicz, and M. Herrrnann. Hanbury-Brown-Twiss analysis in
a solvable model. Phys. Rev., C492442, 1994.

[31] G. F. Bertsch and S. Das Gupta. A guide to microscopic models for intermediate
energy heavy ion collisions. Phys. Rept., 1602189—233, 1988.

151

.1..- no”

 

 

[32] G. F. Bertsch, H. Kruse, and S. D. Gupta. Boltzmann equation for heavy ion colli-
sions. Phys. Rev., C29z673—675, 1984.

[33] Jean-Paul Blaizot and Andre Kr'zywicki. Soft pion emission in high-energy heavy
ion collisions. Phys. Rev., D462246—251, 1992.

[34] N. Borghini, P. M. Dinh, Jean-Yves Ollitrault, A. M. Poskanzer, and S. A. Voloshin.
Effects of momentum conservation on the analysis of anisotropic flow. Phys. Rev.,
C662014901, 2002.

[35] Nicolas Borghini, Phuong Mai Dinh, and Jean-Yves Ollitrault. Flow analysis from
multiparticle azimuthal correlations. Phys. Rev., C642054901, 2001.

[36] Nicolas Borghini, Phuong Mai Dinh, and Jean-Yves Ollitrault. A new method

for measuring azimuthal distributions in nucleus nucleus collisions. Phys. Rev.,
C632054906, 2001.

[37] A. S. Borisov et al. Observation of a high-energy cosmic-ray family caused by a
CENTAURO-type nuclear interaction in the joint emulsion chamber experiment at
the pamirs. Phys. Lett., B190z226, 1987.

[38] R. Brandelik et al. Evidence for charged primary partons in e+e"‘ —> two jets. Phys.
Lett., B1002357, 1981.

[39] David A. Brown and Pawel Danielewicz. Partons in phase space. Phys. Rev.,
D58z094003, 1998. '

[40] Frank R. Brown et al. On the existence of a phase transition for ch with three light
quarks. Phys. Rev. Lett., 65:2491—2494, 1990.

[41] A. Capella, U. Sukhatme, 0] Tan, and J. Tran Thanh Van. Dual parton model. Phys.
Rept., 236:225—329, 1994.

[42] K. C. Chase and A. Z. Mekjian. Exact methods for expectation values in canonical
fragmentation models. Phys. Rev., C52:R2339—R234l, 1995.

[43] Sen Cheng and Scott Pratt. Quantum corrections for pion correlations involving
resonance decays. Phys. Rev., C63z054904, 2001.

[44] Sen Chen g and Scott Pratt. Isospin fluctuations from a thermally equilibrated hadron
gas, 2002. nucl-th/0207051, submitted to Phys. Rev. C.

[45] Sen Cheng et al. The effect of finite-range interactions in classical transport theory.
Phys. Rev., C65:024901, 2002.

[46] S. A. Chin. Transition to hot quark matter in relativistic heavy ion collision. Phys.
Lett., B78z552—555, 1978.

[47] S. K. Choi et al. A measurement of Bose-Einstein correlations in e+e‘ annihilation
at TRISTAN. Phys. Lett., B355z406—414, 1995.

[48] John C. Collins and M. J. Perry. Superdense matter: Neutrons or asymptotically free
quarks? Phys. Rev. Lett., 3421353, 1975.

152

[49] Fred Cooper and Graham Frye. Single particle distribution in the hydrodynamic and
statistical thermodynamic models of multiparticle production. Phys. Rev., D10: 186,
1974.

[50] Fred Cooper, Graham Frye, and Edmond Schonberg. Landau’s hydrodynamic model
of particle production and electron-positron annihilation into hadrons. Phys. Rev.,
D11:192,1975.

[51] Michael Creutz. Quarks, gluons and lattices. Cambridge University Press, Cam-
bridge, New York, Melbourne, 1983.

[52] P. Danielewicz and M. Gyulassy. Dissipative phenomena in quark gluon plasmas.
Phys. Rev., D31253—62, 1985.

[53] Pawel Danielewicz. Determination of the mean-field momentum-dependence using
elliptic flow. Nucl. Phys, A673z375—410, 2000.

[54] Pawel Danielewicz and Scott Pratt. Delays associated with elementary processes in
nuclear reaction simulations. Phys. Rev., C532249—266, 1996.

[55] S. das Gupta and A. Z. Mekjian. Phase transition in a statistical model for nuclear
multifragmentation. Phys. Rev., C5721361—l365, 1998.

[56] C. Derreth, W. Greiner, H. T. Elze, and J. Rafelski. Strangeness abundances in anti-p
nucleus annihilations. Phys. Rev., C31: 1360—1364, 1985.

[57] J. Dias de Deus and S. J adach. On the universality of quark jet fragmentation. Phys.
lett, B70:73, 1977.

[58] H. J. Drescher, M. Hladik, S. Ostapchenko, T. Pierog, and K. Werner. Parton-based
Gribov-Regge theory. Phys. Rept., 350:93—289, 2001.

[59] D. Drijard et al. Density, charge and transverse momentum correlations of particles
in nondiffractive proton-proton collisions at \/3 = 52.5 GeV. Nucl. Phys, B155z269,
1979.

[60] D. Drijard et al. Quantum number effects in events with a charged particle of large
transverse momentum. Part 1. Leading particles in jets. Nucl. Phys, B 1562309, 1979.

[61] D. Drijard et al. Quantum number effects in events with a charged particle of large
transverse momentum. Part 2. Charge correlations in jets. Nucl. Phys, B166z233,
1980.

[62] W. Ehehalt and W. Cassing. Relativistic transport approach for nucleus nucleus
collisions from SIS to SP8 energies. Nucl. Phys, A602:449—486, 1996.

[63] H. T. Elze and W. Greiner. Finite size effects for quark-gluon plasma droplets. Phys
Lett., B179z385—392, 1986.

[64] H. T. Elze, W. Greiner, and J. Rafelski. On the color singlet quark-glue-plasma.
Phys. Lett., B124z515, 1983.

[65] Hans-Thomas Elze, Walter Greiner, and Johann Rafelski. Color degrees of freedom
in a quark-glue plasma at finite baryon density. Z. Phys, C24136l—365, 1984.

153

 

[66] Enrico Fermi. High-energy nuclear events. Prog. Theor. Phys, 52570—583, 1950.

[67] Z. Fodor and S. D. Katz. Lattice determination of the critical point of QCD at finite
t and u. JHEP, 032014, 2002.

[68] Z. Fodor and S. D. Katz. A new method to study lattice QCD at finite temperature
and chemical potential. Phys. Lett., B534287—92, 2002.

[69] Barry A. Freedman and Larry D. McLerran. Ferrrrions and gauge vector mesons at
finite temperature and density. 3. The ground state energy of a relativistic quark gas.
Phys. Rev., D1621169, 1977.

[70] S. J. Gaff, A. Galonsky, C. K. Gelbke, T. Glasmacher, M. Huang, J. J. Kruse, G. J.
Kunde, R. Lemmon, W. G. Lynch, M. B. Tsang, J. Wang, and P. D. Zecher. Time

scales from two-neutron intensity interferometry for the reaction Ar—40 + Ho-165 at
E/A = 25 MeV. Phys. Rev., C58z2161—2166, 1998.

[71] Sean Gavin and Joseph I. Kapusta. Kaon and pion fluctuations from small disori-
ented chiral condensates. Phys. Rev., C652054910, 2002.

[72] Klaus Geiger and Bemdt Miiller. Dynamics of parton cascades in highly relativistic
nuclear collisions. Nucl. Phys, B3692600—654, 1992.

[73] Murray Gell-Mann. A schematic model of baryons and mesons. Phys. Lett., 8:214—
215, 1964.

[74] Murray Gell-Mann and M Levy. The axial vector current in beta decay. Nuovo Cim.,
162705, 1960.

[75] R. Ghetti et al. Neutron-neutron intensity interferometry in E/A = 45 MeV Ni-58 +
Al-27, Ni-nat and Au-197 reactions. Phys. Rev., C622037603, 2000.

[76] S. L. Glashow. Partial symmetries of weak interactions. Nucl. Phys, 222579—588,
1961.

[77] Roy J. Glauber. Coherent and incoherent states of the radiation field. Phys. Rev.,
131:2766—2788, 1963.

[78] Gerson Goldhaber. Multi-pion correlations in e+e" annihilation at SPEAR. In
J. Dias de Deus and J. Soffer, editors, Proceedings of the Int. Conf on High En-
ergy Physics, Lisbon, 1981. European Physical Society, Erice, 1982.

[79] Gerson Goldhaber, Sulamith Goldhaber, Won-Yong Lee, and Abraham Pais. Influ-
ence of Bose-Einstein statistics on the antiproton proton annihilation process. Phys.
Rev., 120:300—312, 1960.

[80] K. Goulianos. Diffractive interactions of hadrons at high-energies. Phys. Rept.,
1011169, 1983.

[81] D. E. Groom et al. Review of particle physics. Eur: Phys. J., C1521, 2000.

[82] D. J. Gross and Frank Wilczek. Asymptotically free gauge theories. 1. Phys Rev.,
D823633—3652, 1973.

154

 

[83] D. J. Gross and Frank Wilczek. Ultraviolet behavior of non-Abelian gauge theories.
Phys. Rev. Lett., 30:1343—1346, 1973.

[84] D. J. Gross and Frank Wilczek. Asymptotically free gauge theories. 2. Phys. Rev.,
D92980—993, 1974.

[85] Miklos Gyulassy, Yang Pang, and Bin Zhang. Transverse energy evolution as a test
of parton cascade models. Nucl. Phys, A6262999—1018, 1997.

[86] Miklos Gyulassy and Xin-Nian Wang. HIJING 1.0: A Monte Carlo program for par-
ton and particle production in high-energy hadronic and nuclear collisions. Comput.
Phys. Commun., 832307, 1994.

[87] M. Y. Han and Yoichiro Nambu. Three-triplet model with double SU(3) symmetry.
Phys. Rev., 1392B1006—B1010, 1965.

[88] R. Hanbury Brown and R. Q. Twiss. A new type of interferometer for use in radio
astronomy. Phil. Mag, 45:663—682, 1954.

[89] R. Hanbury Brown and R. Q. Twiss. A test of a new type of stellar interferometer on
Sirius. Nature, 178:1046—1048, 1956.

[90] D. O. Handzy et al. Understanding proton emission in central heavy-ion collisions.
Phys. Rev. Lett., 75:2916—2919, 1995.

[91] Ulrich Heinz and Barbara V. Jacak. Two-particle correlations in relativistic heavy-
ion collisions. Ann. Rev. Nucl. Part. Sci, 49:529, 1999.

[92] D. Horn and R. Silver. Ann. Phys, 662509, 1971.

[93] P. Huovinen, P. F. Kolb, Ulrich W. Heinz, P. V. Ruuskanen, and S. A. Voloshin.
Radial and elliptic flow at RHIC: Further predictions. Phys. Lett., B503z58—64,
2001.

[94] Sangyong Jeon and Joseph Kapusta. Linear extrapolation of ultrarelativistic nucleon
nucleon scattering to nucleus nucleus collisions. Phys. Rev., C562468—480, 1997.

[95] Sangyong Jeon and Scott Pratt. Balance functions, correlations, charge fluctuations
and interferometry. Phys. Rev., C652044902, 2002.

[96] Ferencz Jiittner. Das Maxwellsche Gesetz der Geschwindigkeitsverteilung in der
Relativtheorie. Ann. Phys. und Chemie, 34:856—882, 1911.

[97] David E. Kahana, Declan Keane, Yang Pang, Tom Schlagel, and Shan Wang. Col-
lective flow from the intranuclear cascade model. Phys. Rev. Lett., 74:4404-4407,
1995.

[98] Frithjof Karsch. Thermodynamics of 2 and 3 flavour ch. AIP Conf. Proc., 6022323—
332, 2001.

[99] J. H. Koch, N. Ohtsuka, and E. J. Moniz. Nuclear photoabsorption and compton
scattering at intermediate-energy. Annals Phys, 154299—160, 1984.

155

[100] P. Koch, B. Muller, and J. Rafelski. Strangeness in relativistic heavy ion collisions.
Phys. Rept., 1422167—262, 1986.

[101] P. F. Kolb, Ulrich W. Heinz, P. Huovinen, K. J. Eskola, and K. Tuorrrinen. Centrality
dependence of multiplicity, transverse energy, and elliptic flow from hydrodynamics.
Nucl. Phys, A696: 197-2 15, 2001 .

[102] P. F. Kolb, P. Huovinen, Ulrich W. Heinz, and H. Heiselberg. Elliptic flow at SP8
and RHIC: From kinetic transport to hydrodynamics. Phys. Lett., B5002232—240,
2001.

[103] G. Kortemeyer, W. Bauer, K. Haglin, J. Murray, and S. Pratt. Causality violations in
cascade models of nuclear collisions. Phys. Rev., C5222714—2724, 1995.

[104] L. D. Landau. On the multiparticle production in high-energy collisions. Izv. Akad.
Nauk SSSR Ser Fiz., 17:51—64, 1953.

[105] C. M. G. Lattes, Y. Fujimoto, and S. Hasegawa. Hadronic interactions of hi gh-energy
cosmic ray observed by emulsion chambers. Phys. Rept., 652151, 1980.

[106] R. Lednicky and T. B. Progulova. Influence of resonances on Bose-Einstein correla-
tions of identical pions. Z. Phys, C55:295, 1992.

[107] Bao-An Li and Che Ming Ko. Formation of superdense hadronic matter in high-
energy heavy ion collisions. Phys. Rev., C5222037-2063, 1995.

[108] Zi-wei Lin and C. M. Ko. Partonic effects on the elliptic flow at RHIC. Phys. Rev.,
C652034904, 2002.

[109] Denes Molnar and Miklos Gyulassy. New solutions to covariant nonequilibrium
dynamics. Phys. Rev., C62:054907, 2000.

[110] Denes Molnar and Miklos Gyulassy. Elliptic flow and freeze-out from the parton
cascade mpc. Nucl. Phys, A6982379—382, 2002.

[111] Dénes Molnar and Miklos Gyulassy. Saturation of elliptic flow at RHIC: Results
from the covariant elastic parton cascade model MPC. Nucl. Phys, A6972495—520,
2002.

[112] Azwinndini Muronga. Second order dissipative fluid dynamics for ultra- relativistic
nuclear collisions. Phys. Rev. Lett., 882062302, 2002.

[113] J. L. Nagle. First results from RHIC: What are they telling us? In Eric Norman,
Lee Schroeder, and Gordon Wozniak, editors, Nuclear Physics in the 21st Century,

volume 610 of AIP Conference Proceedings, pages 70—81. American Institute of
Physics, College Park, 2002.

[114] Y. Nara, N. Otuka, A. Ohnishi, K. Niita, and S. Chiba. Study of relativistic nuclear
collisions at AGS energies from p + Be —> Au + Au with hadronic cascade model.
Phys. Rev. , C612024901, 2000.

[115] Y. Nara, S. E. Vance, and P. Csizmadia. A study of parton energy loss in Au+Au
collisions at RHIC using transport theory. Phys. Lett., B5312209—215, 2002.

156

 

 

[116] K. Niyogi, P. Ghose, and A. K. Ray. Renormalizability and lepton - hadron uni-
versality in pion beta decay in a gauge field theory. Phys. Rev., D1322403—2412,
1976.

[117] Jean-Yves Ollitrault. Flow systematics from SIS to SP8 energies. Nucl. Phys,
A6382195c—206c, 1998.

[118] S. S. Padula and M. Gyulassy. Pion interferometry of O + AU at 200-A/GeV. Nucl.
Phys, A4982555c-560c, 1989.

[119] Sandra S. Padula, Miklos Gyulassy, and Sean Gavin. Pion interferometric tests of
transport models. Nucl. Phys, B3292357, 1990.

[120] Y. Pang, T. J. Schlagel, and S. H. Kahana. Cascade for relativistic nucleus collisions.
Phys. Rev. Lett., 68:2743—2746, 1992.

[121] C. Pinkenburg et al. Elliptic flow: Transition from out-of-plane to in-plane emission
in Au + Au collisions. Phys. Rev. Lett., 83:1295—1298, 1999.

[122] I. U. Pomeranchuk. Dokl. Akad. Novk. USSR, 78:889, 1951.

[123] A. M. Poskanzer et a1. Centrality dependence of directed and elliptic flow at the
SPS. Nucl. Phys, A661:341—344, 1999.

[124] S. Pratt. Pion lasers from high-energy collisions. Phys. Lett., B3012159—164, 1993.

[125] S. Pratt. Validity of the smoothness assumption for calculating two- boson correla-
tions in high-energy collisions. Phys. Rev., C56: 1095, 1997.

[126] S. Pratt. What we are learning from correlation measurements. Nucl. Phys,
A6382125, 1998.

[127] S. Pratt and K. Haglin. Hadronic phase space density and chiral symmetry restora-
tion in relativistic heavy ion collisions. Phys. Rev., C5923304—3308, 1999.

[128] S. Pratt and V. Zelevinsky. Explaining Centauro events by formation of pions in the
isospin singlet channel. Phys. Rev. Lett., 722816—819, 1994.

[129] Scott Pratt. Canonical and microcanonical distributions for Fermi systems. Phys.
Rev. Lett., 8424255—4259, 2000.

[130] Scott Pratt and Subal Das Gupta. Statistical models of nuclear fragmentation. Phys
Rev., C62:044603, 2000.

[131] Johann Rafelski and Bemdt Muller. Strangeness production in the quark-gluon
plasma. Phys. Rev. Lett., 48:1066, 1982. Erratum, ibid 5622334, 1986.

[132] Krishna Rajagopal and Frank Wilczek. Emergence of coherent long wavelength
oscillations after a quench: Application to QCD. Nucl. Phys, B4042577—589, 1993.

[133] Krishna Rajagopal and Frank Wilczek. Static and dynamic critical phenomena at a
second order QCD phase transition. Nucl. Phys, B3992395-425, 1993.

157

 

[134] Ralf Rapp and Edward V. Shuryak. Resolving the antibaryon production puzzle in
high-energy heavy-ion collisions. Phys. Rev. Lett., 8622980—2983, 2001.

[135] Ralf Rapp and Edward V. Shuryak. Regeneration of anti-protons in ultrarelativistic
heavy-ion collisions. Nucl. Phys, A6982587—590, 2002.

[136] Lanny Ray. Correlations, fluctuations and flow at star. In 16th International Confer-

ence on Ultrarelativistic Nucleus-Nucleus Collisions, Nantes, France. 2002. To be
published.

[137] RHIC Transport Theory Collaboration. http://wwwphyduke.edu/rttc/.

[138] Asher Shor and Ronald S. Longacre. Effects of secondary interactions in proton -
nucleus and nucleus-nucleus collisions using the HIJET event generator. Phys. Lett.,
B2182100, 1989.

[139] Edward V. Shuryak. The correlation of identical pions in multibody production.
Phys. Lett., B442387—389, 1973.

[140] Edward V. Shuryak. Instantons in quark plasma, multi — baryon hadrons and neutron
stars. Phys. Lett., B81265, 1979.

[141] Josef Sollfrank et al. Hydrodynamical description of 200-A-GeV/c S + Au colli-
sions: Hadron and electromagnetic spectra. Phys. Rev., C552392—410, 1997.

[142] H. Sorge. Highly sensitive centrality dependence of elliptic flow: A novel signature
of the phase transition in QCD. Phys. Rev. Lett., 82:2048—2051, 1999.

[143] H. Sorge, H. Stocker, and W. Greiner. Relativistic quantum molecular dynamics
approach to nuclear collisions at ultrarelativistic energies. Nucl. Phys, A4982567c—
576C, 1989.

[144] An Tai and Ben-Hao Sa. LUCIAE 3.02 A new version of a computer program
for firecracker model and rescattering in relativistic heavy- ion collisions. Comput.
Phys. Commun., 1162353—365, 1999.

[145] D. Teaney, J. Lauret, and Edward V. Shuryak. Flow at the SP8 and RHIC as a quark
gluon plasma signature. Phys. Rev. Lett., 86:4783—4786, 2001.

[146] D. Teaney, J. Lauret, and Edward V. Shuryak. Hydro + cascade, flow, the equation
of state, predictions and data. Nucl. Phys, A6982479—482, 2002.

[147] Derek Teaney and Raju Venugopalan. Classical computation of elliptic flow at large
transverse momentum. Phys. Lett., B539253-58, 2002.

[148] E. A. Uehling. Transport phenomena in Einstein-Bose and Fermi-Dirac gases. II.
Phys. Rev., 462917—929, 1934.

[149] E. A. Uehling and G. E. Uhlenbeck. Transport phenomena in Einstein-Bose and
Fermi-Dirac gases. 1. Phys. Rev., 432552—561, 1933.

[150] Alexander Volya, Scott Pratt, and Vladimir Zelevinsky. Multiple pion production
from an oriented chiral condensate. Nucl. Phys, A6712617—643, 2000.

158

 

[151] Xin-Nian Wang and Miklos Gyulassy. HIJING: A Monte Carlo model for multiple
jet production in p p, p A and A A collisions. Phys. Rev., D4423501-3516, 1991.

[152] S. Weinberg. Gravitation and Cosmology. John Wiley & Sons, New York, 1972.

[153] S. Weinberg. Conceptual foundations of the unified theory of weak and electromag-
netic interactions. Rev. Mod. Phys, 522515—523, 1980.

[154] K. Werner. Strings, pomerons, and the VENUS model of hadronic interactions at
ultrarelativistic energies. Phys. Rept., 232287—299, 1993.

[155] Urs Achim Wiedemann and Ulrich Heinz. Particle interferometry for relativistic
heavy-ion collisions. Phys. Rept., 319:145—230, 1999.

[156] G. Wolf, G. Batko, W. Cassing, U. Mosel, K. Niita, and M. Schafer. Dilepton pro-
duction in heavy ion collisions. Nucl. Phys, A517:615—638, 1990.

[157] Cheuk-Yin Wong. Dynamics of nuclear fluid. VIII. time-dependent Hartree-Fock
approximation from a classical point of view. Phys. Rev., C2521460—1475, 1982.

[158] E. E. Zabrodin, C. Fuchs, L. V. Bravina, and Amand Faessler. Elliptic flow at collider
energies and cascade string models: The role of hard processes and multi-pomeron
exchanges. Phys. Lett., B5082184—190, 2001.

[159] Bin Zhang. ZPC 1.0.1: A parton cascade for ultrarelativistic heavy ion collisions.
Comput. Phys. Commun., 109:193—206, 1998.

[160] Bin Zhang, Miklos Gyulassy, and Che Ming Ko. Elliptic flow from a parton cascade.
Phys. Lett., B455245—48, 1999.

[161] Bin Zhang, Miklos Gyulassy, and Yang Pang. Equation of state and collision rate
tests of parton cascade models. Phys. Rev., C5821175~1 182, 1998.

[162] Bin Zhang, C. M. Ko, Bao-An Li, and Zi-wei Lin. A multi-phase transport model
for nuclear collisions at RHIC. Phys. Rev., C612067901, 2000.

[163] Bin Zhang and Yang Pang. Frame dependence of parton cascade results. Phys. Rev.,
C5622185—2190, 1997.

[164] J. Zimanyi and T. Csbrgé'. Multi-boson correlations using wave-packets. 12 Exact
results. Heavy Ion Phys, 92241, 1999.

159