MICROSCOPIC CALCULATION
OF CHARGE BALANCE FUNCTIONS
By
Vinzent Steinberg

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Physics—Master of Science
2015

ABSTRACT
MICROSCOPIC CALCULATION OF CHARGE BALANCE FUNCTIONS
By
Vinzent Steinberg
Charge balance functions have been proposed to study the quark-gluon plasma hypothesized to form in relativistic heavy-ion collisions [1]. Motivated by previous work [2],
a microscopic model is developed that considers the evolution of thermal quark-antiquark
pairs in the quark-gluon plasma to calculate hadronic balance functions based on the
chemical properties of the plasma. This involves modeling the generation of the quarks
given the hydrodynamic properties of the quark-gluon plasma combined with what we
know about its chemistry from lattice quantum chromodynamics, as well as modeling the
separation of these quarks and their hadronization. After feeding the resulting hadrons
into a cascade where they scatter and decay, the balance function for any combination of
hadrons species can be calculated.

ACKNOWLEDGEMENTS

I would like to thank all the people who contributed in some way to the work in this
thesis. First and foremost, I thank my advisor, Professor Scott Pratt, who taught me how
to think like a theorist. His help and support have been outstanding, especially during the
final stages of this thesis. His contributions to the work presented here are indispensable.
I would also like to thank my committee members Professor Carl Schmidt and Professor
Gary Westfall for their interest in my work. A special thanks goes out to Debbie Barratt
who was of great help whenever I had administrative problems. I am grateful for the
discussions, insight and friendship my fellow grad students shared with me. They shaped
my experience at MSU.
I would like to thank my family for their unyielding support in everything I do.
I am grateful for the financial support provided by Michigan State University and the
German National Academic Foundation.

iii

TABLE OF CONTENTS
LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Chapter 1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vi
1

Chapter 2 Theory . . . . . . . . . . . . . . . . . . . . . . . . .
2.1 Required Concepts . . . . . . . . . . . . . . . . . . . . .
2.1.1 Bjorken Coordinates . . . . . . . . . . . . . . . .
2.1.2 Modeling Heavy-Ion Collisions . . . . . . . . . .
2.1.2.1 Canonical Picture . . . . . . . . . . . . .
2.1.2.2 Hypersurface Element . . . . . . . . . .
2.2 Generating Quarks . . . . . . . . . . . . . . . . . . . . .
2.2.1 Charge Correlation . . . . . . . . . . . . . . . . .
2.2.2 Quark Pair Production Rate . . . . . . . . . . . .
2.3 Separating Quarks . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Random Walk . . . . . . . . . . . . . . . . . . . .
2.3.2 Causal Diffusion Equation . . . . . . . . . . . . .
2.3.3 Modeling Relativistic Diffusion . . . . . . . . . .
2.3.3.1 Derivation of a Relativistic Propagator
2.3.3.2 Comparison of Propagators . . . . . . .
2.3.3.3 Random Walk vs. Propagator . . . . . .
2.4 Generating Hadrons . . . . . . . . . . . . . . . . . . . .
2.4.1 Hadron Production Rate . . . . . . . . . . . . . .
2.4.2 Maxwell-Jüttner Distribution . . . . . . . . . . .
2.4.3 Hadron Production Rate – Continued . . . . . .
2.4.3.1 Non-Relativistic Limit . . . . . . . . . .
2.4.3.2 Ultra-Relativistic Limit . . . . . . . . .

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

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

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

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

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

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

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

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

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

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

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

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

5
6
6
8
8
8
9
9
12
13
14
15
16
16
18
18
18
18
23
25
26
27

Chapter 3 Implementation . . . . . . . . . . . . . . . . .
3.1 Generation of Quark Pairs . . . . . . . . . . . . . .
3.2 Diffusion of Quarks . . . . . . . . . . . . . . . . . .
3.3 Hadronization . . . . . . . . . . . . . . . . . . . . .
3.4 Hadron Cascade . . . . . . . . . . . . . . . . . . . .
3.5 Calculating Balance Functions . . . . . . . . . . . .
3.5.1 Contribution above Freezeout Temperature
3.5.2 Contribution below Freezeout Temperature

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

28
28
30
33
37
37
38
38

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

Chapter 4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
iv

Appendix A: Notation and Conventions . . . . . . . . . . . . . . . . . . . . . . . . 44
Appendix B: Normalization of Maxwell-Jüttner Distribution . . . . . . . . . . . . 46
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

v

LIST OF FIGURES
Figure 1.1 Diagram summarizing the proposed model to calculate the charge
balance function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

Figure 2.1 Plot of the diffusion propagators as functions of position x and time t
for different space and time step sizes ∆x and ∆t. The continuous lines show
the classical Wiener propagator, the dashed lines the generalized relativistic
one. The dots show the corresponding binomial probabilities to diffuse
from position 0 to x given a random walk. To guide the eye, the dots are
connected by dotted lines. Note that the dots have to be compared to the
area under the corresponding segment of the continous propagators. As
expected, for smaller ∆x the propagators are narrower. Also, the difference
between the relativistic and the non-relativistic propagator becomes smaller
with decreasing ∆x. This is to be anticipated, because the diffusion speed
∆x/∆t = 0.5 is smaller than before (∆x/∆t = 1). . . . . . . . . . . . . . . . . 19
Figure 2.2 Plot of probability p to end up at position x after a one-dimensional
random walk. The blue line is a histogram for x after the specified amount
of steps. The green line is the classical propagator integrated in time, the red
line the relativistic one integrated in time. While the relativistic propagator
is narrower and actually vanishes at the maximal possible x, the results
are very similar. For a high amount of random steps, the random walk
histogram is very close to the propagators. . . . . . . . . . . . . . . . . . . . 20
Figure 3.1 Charge fluctuation χ over entropy s as a function of temperature T
for various quark pairs. The dotted vertical line represents the freezeout
temperature at 165 MeV. The solid lines represent the results of lattice QCD
simulations of a quark-gluon plasma [32]. The dashed lines show the results of assuming a hadron gas in thermal equilibrium (2.13). Note that
for high temperatures χ/s is roughly constant according to lattice QCD,
and the offdiagonal element (χus ) vanishes as expected for the parton gas
model (2.14). The hadron gas model roughly agrees with the lattice data
for small temperatures. For high temperatures it does not, in that limit we
expect the plasma to behave like a parton gas where all light species are the
same, because their masses are siginifcantly smaller than the temperature
(mu ≈ md ≪ ms ≈ 100 MeV). This expectation is confirmed by the lattice
data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
Figure 3.2 Temperature T as a function of space ( x, y) in the transverse plane.
The solid line represents the freezeout temperature at 165 MeV. This was
taken at the beginning of an Israel-Stewart hydrodynamics simulation (see [33]
for the implementation). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
vi

Figure 3.3 Example of a random walk as implemented here. This shows the
diffusion trajectories of a random sample of quarks from the first population
until they freeze out. The trajectories are projected to the transverse plane. . 34
Figure 3.4 An illustration of the modified Bresenham algorithm used. Each
square cut by the diffusion trajectory has to be considered, because it is a
possibility for a freezeout. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

vii

Chapter 1
Introduction
The quark-gluon plasma is a phase consisting of partons that is believed to form at high
temperatures and densities, comparable to the first few microseconds after the Big Bang. In
relativistic heavy-ion collisions at RHIC (Relativistic Heavy Ion Collider) and at the LHC
(Large Hadron Collider) matter is created with energy densities near 10 GeV/fm3 . This
is above the threshold for defining individual hadrons. The chemistry of an equilibrated
quark-gluon plasma is well understood: There are 3 light quark species and 8 gluons with
a total of 52 degrees of freedom. Such a plasma can be modeled using lattice quantum
chromodynamics. However, there has been little experimental evidence that the matter
created in heavy-ion collisions has chemical properties close to those of equilibrated matter.
It is not possible to measure the quarks in the plasma directly, which makes it challenging
to infer the chemistry of the partonic matter.
A promising approach at looking through this haze of hadronization are charge balance
functions. They measure the pairwise correlations between charges as a function of relative rapidity. The correlations are driven by the local charge conservation in heavy-ion
collisions. Hadrons are produced in pairs of opposite charge that are separated in rapidity
due to thermal motion and diffusion into regions of different collective velocity. If their
creation is delayed they will have less time to separate, yielding a contribution to the

1

balance function for small relative rapidities, whereas hadrons created early contribute for
larger relative rapidites.
Balance functions have first been suggested for investigating hadronization in jets produced by proton-proton [3, 4] and electron-positron [5] collisions. Later it was proposed to
use them for testing the quark-gluon plasma hypothesis [1], in particular for understanding when the quarks are created and how they are diffused. The existence of a quark-gluon
plasma implies delayed hadronization, which would make the balance functions narrower
for the reasons discussed above.
The STAR collaboration measured balance functions for gold-gold, deuteron-gold and
proton-proton collisions at RHIC [6–10]. They are defined for each event as

Bαα¯ (∆y) =

1
2

Nαα¯ (∆y) − Nαα (∆y) Nα¯ α (∆y) − Nα¯ α¯ (∆y)
+
Nα
Nα¯

(1.1)

where α is a hadron species and α¯ the corresponding anti-species (for example π + and
π − ). Nαβ is the number of αβ pairs separated in rapidity by ∆y. Nα is the total number
of α particles. For central collisions the balance functions were narrower, as expected if
there is a quark-gluon plasma. Only for high impact parameters was the width in agreement with models based on independent nuclear scattering [11, 12]. This confirmed the
predictions made in [1]. Additionally, the experimental results have been compared with a
microscopic hadronic model and a blast-wave model, where only the latter could correctly
describe the spectra [13].
Similarly, balance functions have been measured for proton-proton, carbon-carbon and
lead-lead collisions at CERN SPS [14]. It was shown that the width of the balance function
decreased with increasing centrality and nucleus size, again confirming the predictions
stemming from the assumption of delayed hadronization due to the existence of a quarkgluon plasma.
The small width of the balance function can also be explained with a coalescence

2

rehadronization model [15]. Furthermore, pion balance functions have been calculated in
good agreement with the STAR data using a thermal model that considers decays [16].
Balance functions can be binned by the azimuthal separation ∆φ. This has been proposed to constrain the amount of transverse flow at freezeout [17, 18] and can be calculated
from two-particle angular correlations measured at RHIC [19].
The models discussed so far did not consider balance functions for different species.
The first one that did assumed the existence of two charge creation waves: one when the
quark-gluon plasma forms and one when it hadronizes [2]. For a chemically equilibrated
quark-gluon plasma one would expect that most of the electric charge is created in the
second wave, while most of the strangeness is created in the first one. Because pions carry
most of the charge and kaons most of the strangeness, this implies that the K + K − balance
function is broader than the π + π − balance function. For protons the balance function
should be dominated by baryon conservation. Because of the drop of baryon susceptibility at hadronization, the second wave contributes negatively to the p p¯ balance function,
resulting in a dip at small relative rapidity. Fitting the model to preliminary STAR results [10] showed that the data exposes a chemical makeup that agrees within a 20% error
margin with results from lattice simulations of a quark-gluon plasma. By considering all
three balance functions, it was made clear that assuming the creation of charge at one time
is inconsistent with the data.
The goal of this thesis is to develop a model in the spirit of [2] that uses a microscopic
simulation of the quark pairs instead of the simple parametric model assuming two charge
creation waves. This should be more accurate, since quark pairs are continously produced
during the evolution of the reaction. Such a model could improve our understanding of
the microscopic properties of the quark-gluon plasma, for example it might be possible to
determine its diffusion constant.
The canonical picture of a relativistic heavy-ion collision is the following: During the
collision a quark-gluon plasma is formed that can be described as a hydrodynamic region.

3

lattice QCD
quark generation

diffusion

hadronization

cascade

hydrodynamics

Figure 1.1

Diagram summarizing the proposed model to calculate the charge balance function.

The plasma evolves and cools down, until it is cold enough for hadronization to take place
(at ca. 165 MeV). The resulting hadrons scatter and decay in a cascade before they can be
observed in a detector. In this picture the proposed model can be roughly divided into the
following steps:
1. Generate up, down and strange charges in the hydrodynamic phase.
2. Model the separation of the charges by diffusion.
3. Translate the charges into hadrons at the interface of the hydrodynamic region and
the cascade.
4. Scatter and decay hadrons.
After performing those steps, the balance functions can be directly calculated from the
hadrons. This thesis tries to address step 1, 2 and 3, for which the theory is discussed
in chapter 2. The implementation details of the model and how everything fits together is
explained in chapter 3.

4

Chapter 2
Theory
There are two contributions to the balance function that have to be considered:
1. The contribution due to the quark pairs generated above the freezeout temperature that
freeze out to correlated hadrons.
2. The contribution due to uncorrelated hadrons below the freezeout temperature that
decay and annihilate each other while being spread in space by collisions.
Modeling the first contribution is the topic of this chapter. For the second contribution an
existing cascade code can be used [20], which will be discussed in more detail in chapter 3.
This chapter provides the background and theoretical description of a method for
generating and propagating correlated quark-antiquark pairs during the evolution of the
heavy-ion collision. We then show how those correlated quarks can be projected onto
correlated hadrons, and ultimately into measurable correlations of the final state.
In the first section a very convenient set of coordinates is introduced, which are used
for defining the balance functions. They also give some insight into the hydrodynamic
evolution of the plasma.
The quark pairs contribute to the balance functions in three steps:
1. The pairs are generated according to thermal equilbrium (section 2.2).
5

2. The pairs separate in a diffusive process (section 2.3).
3. The individual quarks hit the freezeout surface and hadronize (section 2.4).
This yields hadrons that can be correlated by tracing back the original quark pairs, enabling
the calculation of balance functions.

2.1

Required Concepts

2.1.1

Bjorken Coordinates

In highly relativistic heavy-ion collisions it can be observed that in the center-of-momentum
frame small boosts should not influence the observables, because the colliding nucleons
are at much larger rapidities than the thermal velocities of the created particles [21]. In
such a boost-invariant system the flow is not accelerated in the longitudinal direction, because there are no gradients able to cause an acceleration. This means that the longitudinal
component of the velocity of a hydrodynamic fluid element is approximately constant.
Motivated by these considerations, we define the proper time τ and the spatial rapidity η:

τ:=

⇔

η : = tanh−1

t2 − z2

z = τ sinh η

t = τ cosh η

z
t

(2.1)
(2.2)

Under boosts in z-direction the proper time τ is invariant and the rapidity η is shifted
by a constant. As we will see later in this section, local quantities depend only on τ and
are independent of η.
For small η and small z, Equation 2.2 can be approximated by

⇔

z
t

τ=t

η=

z = τη

t=τ.
6

(2.3)
(2.4)

This corresponds to the ultrarelativistic limit, where z will be small due to the Lorentz
contraction.
The hydrodynamic equation for the entropy density s is given by the continuity equation:
0 = ∂µ (suµ ) = s∂µ uµ + uµ ∂µ s ,

(2.5)

where u is the 4-velocity of the fluid. This can be rewritten in Bjorken coordinates using
that τ is the time in the rest frame and that uz = τz (no acceleration):
uµ ∂µ = ∂τ

∂µ uµ =

1
,
τ

(2.6)

resulting in a much simpler equation for the entropy:

0=

s
+ ∂τ s
τ

(2.7)

This yields a solution that does not depend on η:
τ
s(τ ) = s(τ0 ) 0
τ

(2.8)

If we assume a non-interacting gas of massless partons, then the density of each parton
species scales with the entropy. In an isentropic expansion of this gas the number of each
species must thus be conserved until hadronization.
By assuming an equation of state relating the energy density ǫ and the pressure p (for
instance p = ǫ/3 for an ideal gas of massless partons), it can be shown that the temperature
and energy density similarly only depend on τ.

7

2.1.2

Modeling Heavy-Ion Collisions

2.1.2.1

Canonical Picture

When simulating relativistic heavy ion collisions, it is appropriate to model the quarkgluon plasma using hydrodynamics as long as the collisions are sufficiently rapid so that
all the species move with a single collective velocity and a single local kinetic temperature.
At some point of the reaction (or far away enough from the center) the plasma is cold
enough that those conditions are no longer met. The hadrons will then move independently, and a microscopic model based on binary interactions should take over. Those
models are usually called "hadronic cascades". For many particles (or high oversampling)
they approach a Boltzmann description [22].
2.1.2.2

Hypersurface Element

The decoupling temperature defines a surface separating the hydrodynamic description
and the cascade interface. At this surface hadrons are emitted from the hydrodynamic
region into the cascade.
We define a hypersurface A using a scalar function C of the 4-vector1 x:
A : = { x ∈ R 4 | C ( x ) = 0}

(2.9)

The hypersurface is described by local quantities, usually the temperature or density, e.g.
C (x) = T (x) − T0 . If C is positive inside the hydrodynamic region and negative outside
of it, then the corresponding hypersurface element per 4-volume is
Ωµ (x) = −∂µ θ C (x) = −δ C (x) ∂µ C (x) .

(2.10)

Ω is orthogonal to the surface. Multiplying C by a positive constant does not change Ωµ .
1 The following notation will be used in this document: 4-vectors are typeset in bold (x), 3-vectors have
an arrow on top (x) and vector components use the normal font style (t, x, y, z).

8

For C (x) = T (x) − T0 this definition implies that Ωµ is parallel to −∂µ T.
Given a 4-momentum distribution f (p) at the boundary, the number of particles emitted through a hypersurface element Ω is given by the Cooper-Frye formula [23]:

dN =

d4 x

Ω(x) · p
d3 p
f (p)
Ep
(2π )3

(2.11)

We will need this later to calculate how many hadrons are to be generated during freezeout.

2.2

Generating Quarks

In this section a formula for the source function of thermal quark pairs is derived. These
correlated pairs will ultimately be used to generate the charge balance correlations seen in
experiment. For this, it is helpful to introduce charge fluctuation for quarks which can be
related to the charge correlations of hadrons. The source function enables us to generate
quark pairs in the quark-gluon plasma given its hydrodynamic evolution.

2.2.1

Charge Correlation

Consider the correlation of two quark charges a, b separated by the spatial rapidity η. Here
we only care about up, down and strange quarks, since the other quarks are too heavy to
exists in significant quantities. Then their charge fluctuation is given for a box of volume V
by
χ ab :=

Q a Qb − Q a
V

Qb

(2.12)

For zero net charge, the last term vanishes. This is a reasonable assumption for many
measurements of heavy-ion collisions at LHC or at the highest RHIC energies.
For a non-interacting hadron gas we can assume that only quarks confined to the same

9

hadron are correlated, yielding

χ ab =

∑
α

nα qα,a qα,b ,

(2.13)

where the sum goes over all hadron species α and nα is the hadron number density.
For a non-interacting parton gas, the quarks are only correlated with themselves. Assuming unit charge, we get

χ ab = δab

n a + n a¯

.

(2.14)

The charge correlation for a given difference in spatial rapidity η is

gab (η ) := ρ a (0)ρb (η ) − ρ a (0)

ρb (η )

(2.15)

where ρ a is the density of a charges. That is, the difference of the number density of
a quarks and the number density of a¯ quarks.
If the net charge is zero, the last term vanishes as before. χ ab corresponds to the charge
correlation at η = 0:
ǫ

χ ab =

0

gab dη

for small ǫ .

(2.16)

For an equilibrated system χ ab can be taken from calculations based on the grand canonical
ensemble, such as lattice gauge theory.
If the local correlation is equilibrated and if the spatial extent of the correlation is
smaller, the charge flux can be expressed in terms of the correlation. Since the net charge
is zero
∞
0

gab dη = 0 ,

(2.17)

gab dη = −χ ab .

(2.18)

this implies
∞
ǫ

10

Since ρα refers to the density of a conserved charge, it may be predicted by a model, i.e.
quantumchromodynamical simulations on a lattice. This will be discussed in more depth
in section 3.1.
Once we know gab , we need to understand how this translates into a correlation between hadrons. We quantify this correlation between hadron species α, β separated by
spatial rapidity η:

Gαβ (η ) :=

nα (0) − nα¯ (0) n β (η ) − n β¯ (η )

(2.19)

This is directly related to the quark charge correlation:

gab (η ) =

∑ Gαβ (η )qα,a q β,b

(2.20)

αβ

where qα,a is the a-charge of hadron species α. As shown by Pratt [24], for vanishing net
charge we can relate Gαβ to the properties of the quarks:
Gαβ (η ) = 4

∑
abcd

1 (0) g ( η ) n q χ −1 ( η )
nα qα,a χ−
β β,d cd
bc
ab

(2.21)

This equation links the measurable charge correlations of the hadrons with the properties
of the quark-gluon plasma, which cannot be directly observed.
The charge balance function can be calculated from the correlation function by dividing
by the multiplicity:
Bαβ (η ) =

Gαβ (η )
n β + n β¯

(2.22)

where n β is the number of β hadrons per unit rapidity. It represents the additional
¯ Note that the normalization is
density of seeing an α − α¯ given one observed an extra β − β.
slightly different from the one used for the experimentally measured balance function (1.1).
However, for vanishing net charge we often have nα

11

≈ nα¯ , which makes the two

definitions (2.22) and (1.1) identical.
Balance functions have been measured for various species [6–10, 14, 19].

2.2.2

Quark Pair Production Rate

In the previous section it was shown how gab could be used to caculate correlations between hadrons. To calculate gab , we need to model the creation of balancing charges and
their separation. In this section we derive a formula for the pair production rate. We will
assume that we are given the hydrodynamic evolution of the the quark-gluon plasma.
That is, we know the solution to the hydrodynamic equations, in particular: the entropy
density s, the particle number density n and the fluid 4-velocity u as functions of the
4-position x.
Those solutions depend on the initial conditions of the hydrodynamic problem, which
are essentially the parameters of the heavy-ion collision, such as beam energy, centrality
or the type of the colliding nuclei.
For sources σ and ν, continuity demands for the entropy current s
∂µ (suµ ) = σ

⇔

uµ

∂µ s
σ
= − ∂µ uµ
s
s

(2.23)

⇔

uµ

∂µ n
ν
= − ∂µ uµ
n
n

(2.24)

and the particle current n:
∂µ (nuµ ) = ν

Combining (2.23) and (2.24) gives
uµ ∂µ

n
n
=
s
s

uµ ∂µ n uµ ∂µ s
−
n
s

12

=

n
s

ν σ
−
n
s

.

(2.25)

Assuming there is a particle source (ν = 0), but entropy is conserved (σ = 0), this gives
suµ ∂µ

n
=ν.
s

(2.26)

Identifying n with −χ ab , we obtain the number of correlated quark pairs
dNab = ν d4 x = −suµ ∂µ

χ ab 4
d x.
s

(2.27)

In the fluid rest frame this simplifies to

dNab = −s∂t

χ ab 4
d x.
s

(2.28)

This means that there are no pairs created if the charge fluctuation per entropy is constant
in time.
χ
The ratio sab as a function of temperature T can be obtained from the results of lattice

quantum chromodynamics calculations or by assuming a thermalized hadron gas. The
temperature as a function of spacetime is given by the hydrodynamical model. Combining
both, we can calculate the rate of pair generation for each hydro cell. All this will be
discussed in more detail in chapter 3.

2.3

Separating Quarks

In the previous section we derived which quark pairs should be generated where. Each
pair is created at one point in spacetime and the two partner quarks are correlated. Since
their separation is zero, this implies a contribution to the two-particle correlation function
corresponding to the type of the pair at η = 0. However, due to their considerable thermal
energy, the two quarks will move apart, yielding contributions at η = 0.
We will model this separation as a diffusive process. Similarly to the hydrodynamic

13

model of the quark-gluon plasma, this approach is only valid if the particles are dense
enough to collide with each other multiple times. Their mean free path should be smaller
than the size of the fireball. The aim of this section is to find a way to diffuse the generated
thermal quarks that is consistent with relativity.
The first subsection shows how a random walk and the non-relativistic diffusion equation are related. This gives a way to express the diffusion constant in terms of the random
walk parameters. In the second subsection, a simple causal generalization of the diffusion
equation is derived. However, this approach has some issues which are discussed. The
last subsection compares an approach remedying those issues (as suggested in [25]) to a
very simple relativistic random walk. This serves as a verification that the random walk is
a valid approach for modeling the relativistic diffusion of the quarks.

2.3.1

Random Walk

Consider a one-dimensional random walk with spatial step size l and time step size τ. Let
p be the probability to diffuse to the right. (1 − p is then the probability to diffuse to the
left.) Then the probability to be at position x at time t + τ is given by:

⇒

Pt+τ ( x ) = pPt ( x − l ) + (1 − p) Pt ( x + l )

(2.29)

Pt+τ ( x ) − Pt ( x ) = pPt ( x − l ) + (1 − p) Pt ( x + l ) − Pt ( x )

(2.30)

Assuming isotropy, we have p = 0.5:

Pt+τ ( x ) − Pt ( x ) =

1
P ( x − l ) + Pt ( x + l ) − 2Pt ( x )
2 t

(2.31)

Those differences can be rewritten as derivatives

τ∂t P =

1 2 2
l ∂ P,
2 x

14

(2.32)

or
∂t P = D∂2x P

where

D :=

l2
1
= vl,
2τ
2

v :=

l
.
τ

(2.33)

D is called diffusion constant. This connects the discrete random walk to the continuous
diffusion equation.
Let x0 denote the initial 4-position of the random walk. Solving the classical diffusion
equation for P(x|x0 ) = δ(x − x0 ) gives rise to the Wiener propagator
P ( x | x0 ) =

1
4πD (t − t0 )

exp −

( x − x0 )2
4D (t − t0 )

for

t > t0 .

(2.34)

In one dimension P(t, x |t0 , x0 )dx gives the probability that a particle that was at position x0
at time t0 can be found in the element [ x, dx ] at time t > t0 . Equation 2.34 shows that this
probability is never zero, even for | x − x0 | > c(t − t0 ). Thus this propagator violates
causality. The following two subsections discuss possible solutions to this problem.
Note that the integral of the propagator over x is independent of time. If we interpret
P as a particle density instead of a probability, then this integral describes the net number
of particles.

2.3.2

Causal Diffusion Equation

One possible way to make the diffusion equation causal is to introduce a relaxation time tr .
Assume the phase space distribution f (x, p) cools towards an equilibrium f e :
dt f = ∂ t f + v ∇ f = −

15

f − fe
tr

(2.35)

When considering small tr and oscillating perturbations of f with small wave numbers,
this implies the following generalization of Fick’s law2 :

tr ∂ t j + j = − D ∇ n ,

(2.36)

where n is the particle density. Combining this with the continuity equation,

∂t n + ∇ j = 0 ,

(2.37)

tr ∂2t n + ∂t n = D ∇2 n

(2.38)

gives the causal diffusion equation:

Information propagates at the speed v =

√

D/tr . For tr = 0 we get the ordinary diffusion

equation. This equation is however not causal, because information travels at infinite
speed.
There are some shortcomings when applying this approach to relativistic diffusion:
The equation exhibits singular (i.e. δ-peaked) diffusion fronts moving with speed v. For
a massive particle this implies that a small fraction of the particle carries a large amount
of the kinetic energy (much larger than mc2 ), which is unlikely to occur in nature. See
Dunkel, Talkner, and Hänggi [25] for reference.

2.3.3

Modeling Relativistic Diffusion

2.3.3.1

Derivation of a Relativistic Propagator

The classical diffusion propagator (2.34) violates causality and is not Lorentz-invariant.
Dunkel, Talkner, and Hänggi [25] propose a relativistic generalization. A short reproduction of their approach is given here:
2 This is sometimes called Maxwell-Cattaneo relation or telegraph equation.

16

First, consider the classical case. The action per mass of a particle traveling from 4position x0 to x is given by
a(x, x0 ) =

1 t
v(t′ )2 dt ,
2 t0

(2.39)

where the velocity v is approximately constant between the scatterings during diffusion.
This becomes minimal if there are no collisions and the velocity is constant ∀ t′ ∈ [t0 , t]:
v=

x − x0
t − t0

⇒

a− (x, x0 ) =

1 ( x − x0 )2
.
2 t − t0

(2.40)

In the non-relativistic case the velocity is unbound, so the maximal action a+ will be
infinity. This lets us rewrite the propagator (2.34) as an integral over all possible actions:

a(x, x0 ) ∼

a+ (x,x0 )
a− (x,x0 )

exp

a
2D

da

(2.41)

Generalizing (2.39), we get the relativistic action

a(x, x0 ) = −

t
t0

1 − v(t′ )2 dt′ ,

(2.42)

yielding new maximal and minimal a± in analogy to before:
a− (x, x0 ) = −

(t − t0 )2 − ( x − x0 )2 = −|x − x0 | ,

a+ (x, x0 ) = 0 .

(2.43)
(2.44)

Plugging this into (2.41) gives a relativistic, manifestly Lorentz-invariant generalization of
the Wiener propagator:

P ( x | x0 ) ∼




exp

| x − x0 |
2D

− 1 for ( x − x0 )2 ≤ (t − t0 )2



0

otherwise

17

(2.45)

The normalization is time-dependent, see [25] for details.

2.3.3.2

Comparison of Propagators

Figure 2.1a compares the different propagators. It is visible that for large times the difference between the relativistic and the non-relativistic variant are small. For smaller
diffusion speeds, the difference becomes even smaller (see Figure 2.1b).
The discrete binomial propagator of a random walk agrees very roughly with the
propagators. A random walk does not violate causality, because it can respect the finite
diffusion speed by choosing only step sizes in space (∆x) and time (∆t) such that ∆x
∆t ≤ 1.
2.3.3.3

Random Walk vs. Propagator

In Figure 2.2 the two models (relativistic and non-relativistic propagator) are compared
with the histogram of a random walk. It can be seen that the models and the simulation
roughly agree, especially for a large amount of steps.
This justifies using the random walk to model relativistic diffusion.

2.4

Generating Hadrons

When a diffused quark hits the freezout surface, hadronization takes place. In this section
the production rate of such hadrons is derived, which will be used to generate the hadrons
in a Monte Carlo simulation. We know for each hadron from which quark pair it originated, and thus with which other hadrons it is correlated. This enables us to calculate the
corresponding balance function.

2.4.1

Hadron Production Rate

In order to calculate the number of hadrons ∆Nα emitted due to a quark diffusing through
the freezout surface, we need to consider the statistical thermal equilibrium first:
18

1.0

comparison of classical, relativistic and binomial propagator
for ∆x =1 and ∆t =1
t =1
t =4
t =10

0.8

t =1

propagator

t =4
t =10

0.6

t =1
t =4
t =10

0.4
0.2
0.0

10

5

0

5

10

x

(a) ∆x = 1 and ∆t = 1.

comparison of classical, relativistic and binomial propagator
for ∆ =0 5 and ∆ =1
1.0
x

.

t

=1
=4
t =10
t =1
t =4
t =10
t =1
t =4
t =10
t
t

propagator

0.8
0.6
0.4
0.2
0.0

10

5

0

5

10

x

(b) ∆x = 0.5 and ∆t = 1.
Figure 2.1 Plot of the diffusion propagators as functions of position x and time t for different space and time
step sizes ∆x and ∆t. The continuous lines show the classical Wiener propagator, the dashed
lines the generalized relativistic one. The dots show the corresponding binomial probabilities
to diffuse from position 0 to x given a random walk. To guide the eye, the dots are connected
by dotted lines. Note that the dots have to be compared to the area under the corresponding
segment of the continous propagators. As expected, for smaller ∆x the propagators are narrower.
Also, the difference between the relativistic and the non-relativistic propagator becomes smaller
with decreasing ∆x. This is to be anticipated, because the diffusion speed ∆x/∆t = 0.5 is smaller
than before (∆x/∆t = 1).

19

random walk with ∆x =1 vs. models
random walk
Wiener process
relativistic diffusion

0.40
0.35
0.30

p

0.25
0.20
0.15
0.10
0.05
0.00

10

5

0

x

5

10

(a) 10 random steps.

0.08
0.07
0.06

random walk with ∆x =1 vs. models
random walk
Wiener process
relativistic diffusion

0.05
p

0.04
0.03
0.02
0.01
0.00
0.01
100

50

0

x

50

100

(b) 100 random steps.
Figure 2.2 Plot of probability p to end up at position x after a one-dimensional random walk. The blue line
is a histogram for x after the specified amount of steps. The green line is the classical propagator
integrated in time, the red line the relativistic one integrated in time. While the relativistic
propagator is narrower and actually vanishes at the maximal possible x, the results are very
similar. For a high amount of random steps, the random walk histogram is very close to the
propagators.

20

Thermal equilibrium can be formulated as a constrained entropy maximization. Mathematically this requires Lagrange multipliers. For probabilities f i the Gibbs entropy is given
by ∑i f i log f i , which has to be maximized while conserving the following quantities:
1. total probability ∑i f i ,
2. average energy ǫ and
3. charge flux through the freezeout surface
with corresponding Langrange multipliers α, β, γ.

max − ∑
fi

⇒
⇔

i

f
f i log f i + α f i + β f i ǫi + γq i pi · Ω
ǫ

(2.46)

i

p ·Ω
0 = − ∑ log f i + 1 + α + βǫi + γq i
ǫi
i

(2.47)

pi · Ω
ǫi
i

(2.48)

∏ f i = exp
i

−1 − α − β ∑ ǫi − γq ∑
i

(Note that p · Ω = E p Ω0 − pΩ.)
This means that the requirement of the conservation of charge flux through the hypersurface introduces a factor of exp(γq p · Ω/E p ). We need to consider different quark
species a, so the γ term should depend on a.
Consequently we get the change in the number of hadron species α in terms of the
a-quark charge qα,a of the hadron and the Lagrange multiplier γ:

d(∆Nα ) = dNα

exp

∑ γa
a

≈ dNα ∑ γa
a

p·Ω
qα,a − 1
Ep

p·Ω
qα,a
Ep

21

(2.49)
(2.50)

(If we consider p · Ω/E p = 1 in the exponential correction, then we obtain Nα instead of
ξ α . The multiplier γa is then usually defined terms of the change in chemical potential ∆µ a
and temperature T:
γa

p·Ω/E p =1

∆µ a
T

(2.51)

∆µ a qα,a
T
a

(2.52)

=

This yields
d(∆Nα ) ≈ dNα ∑
instead.)
In our case, dNα can be calculated using the momentum distribution f α and the CooperFrye formula [23]:
1 d3 p
dNα =
p · Ω f α (p)
(2π )3 E p

where

p = ( E p , p)

(2.53)

Substituting this gives:
1
ξα : =
(2π )3
χˆ ab : =

d3 p

f α ( p)

p·Ω 2
Ep

where

∑ ξ α qα,a qα,b = χˆ ba

Ep
u·p
− T
− T
= ge
f α ( p) := ge

(2.54)
(2.55)

α

∆Nα = ξ α qα,a γa

(2.56)

Note that the normal phase space density contributes one factor of p · Ω/E p , the additional
Lagrange multiplier γ (corresponding to the conservation of charge through the freezeout
surface) contributes another one. This additional factor is the difference in the definition
of χ ab compared to subsection 2.2.1, where Nα as used instead of ξ α . The moment distribution f α chosen here is proportional to the degeneracy g and will be discussed in the
next section.

22

Conservation of b-quark charge qb demands that

∑ qα,b ∆Nα

(2.57)

≈ ∑ ξ α γa qα,b qα,a

(2.58)

= ∑ χˆ ba γa .

(2.59)

∆Qb =

α

α,a

a

Inverting the matrix product yields the Lagrange multiplier

γa =

∑ χˆ −ab1 ∆Qb ,

(2.60)

b

which can be substituted into the initial expression:
1 q ∆Q
∆Nα = ξ α ∑ χˆ −
b
ab α,a

(2.61)

ab

Here ∆Nα is the additional number of α hadrons due to an extra charge ∆Qb flowing
through the hypersurface. Note that for antiparticles α¯ we have ∆Nα¯ = −∆Nα , and that
∆N is independent of the magnitude of the hypersurface element, because ξ ∼ |Ω| and

| χ −1 | ∼ ξ −1 ∼ | Ω | −1 .

2.4.2

Maxwell-Jüttner Distribution

In the previous subsection an expression for the hadron production rate at the freezeout
surface was derived. To be able to calculate it, we need to choose an appropriate momentum distribution f α .
After the freezeout the densities should be low enough that we can assume the hadrons
are non-interacting particles in thermal equilibrium. In this case, the momentum follows
the Maxwell-Jüttner distribution [26], a relativistic generalization of the Maxwell-Boltzmann

23

distribution:
f (p) =

g

(2.62)

e(u·p−µ)/T + λ

where g is the degeneracy, µ the chemical potential in equilibrium and




1




λ = −1






 0

for Fermi-Dirac
(2.63)

for Bose-Einstein
for Maxwell-Boltzmann

statistics. Note that the degeneracy for a spin s particle is given by the number of possible
spin states:
g = 2s + 1

(2.64)

Assuming the ultrarelativistic case eu·p/T ≫ 1 and u · p ≫ µ , Equation 2.62 simplifies
to the Maxwell-Boltzmann distribution:
f (p) = g e−u·p/T

(2.65)

For pions, the Bose-Einstein distribution would differ at the 5% level. Integrating gives
the average particle number:

n=

f (p)

d3 p
gd
=
(2π )3
(2π )3

where

d := 4πTm2 K2

m
T

(2.66)

where K2 a modified Bessel function of the second kind, defined more generally [27, 28]
by the integral
1
∀a > − :
2

Ka (ζ ) :=

ζ a Γ 12
2
Γ a + 12

∞
1

e−ζy (y2 − 1)

a− 12

dy .

√
(Note that Γ 21 = π.) See section 4.2 for the derivation of the normalization d.
24

(2.67)

2.4.3

Hadron Production Rate – Continued

Now that we have picked a momentum distribution f α , we can calculate the integral from
subsection 2.4.1:

ξα : =

g
(2π )3

p·Ω 2
u·p
exp −
Ep
T

d3 p

Ω20 gd

g
+
=
3
(2π )
(2π )3

d3 p

(2.68)

pz Ωz 2 − E p
e T
Ep

(2.69)

Without loss of generality, we chose coordinates where

u = (1, 0, 0, 0)

Ω = (Ω0 , 0, 0, Ωz ) .

and

(2.70)

See Equation 2.66 for the definition of d. Since the integral is symmetric to the choice of z,
we can replace pz with p:
Ω20 gd

Ω2z g
+
ξα =
(2π )3 3(2π )3

d3 p

p2 − E p
e T
E2p

(2.71)

This integral is similar to the normalization of the Jüttner distribution. The same coordinate
transformation can be used here:
d3 p

∞
p2 − E p
T = 4π
e
0
E2p

∞ (γ2 − 1)3/2 − m γ
p4 − E p
3
T
e T dγ
e
dp = 4πm
γ
1
E2p

(2.72)

If it weren’t for the factor 1/γ, this could be expressed as a modified Bessel function
using (2.67).
Mathematica [29] could not compute this integral. SymPy [30] returns a result in terms
of Meijer G-functions [31]:
∞ (γ2 − 1)3/2 − m γ
e T dγ
γ
1

=



1
3 3,0 
G1,3 
8
− 3 , 0, 1
2

25

2

m2





4T 2

(2.73)

where


 a1 , . . . , a p
Gm,n
p,q 
b1 , . . . , bq



n
∏m
1
j =1 Γ ( b j − s ) ∏ j =1 Γ (1 − a j + s )

zs ds
z :=
p
2πi L ∏q
Γ
(
1
−
b
+
s
)
Γ
(
a
−
s
)
∏ j = n +1
j
j
j = m +1

(2.74)

for
• m, n, p, q ∈ Z :
• ak − b j ∈ N

m ∈ [0, q], n ∈ [0, p]

∀ k ∈ {1, . . . , n}, j ∈ {1, . . . , m}

• z=0 .
L represents the line in the complex plane along which the integration is performed. There
are three different options, see [31] for more details.
In practice it is difficult to find implementations that numerically calculate the Meijer
G-function. Since implementing them is non-trivial, the easiest approach to evaluate ξ α is
numerical integation. We are only interested in ξ α at the freezeout surface, so we can get
away with calculating the integral only once for every hadron species.

2.4.3.1

Non-Relativistic Limit

In the non-relativistic limit the temperature is much smaller than the mass of the hadron:

T≪m

⇒

p
≈v
Ep

(2.75)

We can also assume the hadrons behave like an ideal gas. The equipartition theorem then
gives
1
3
m v2 = T
2
2

⇔

26

v2 =

3T
m

(2.76)

Applying those results to (2.69):
Ω20 d
ξα
Ω2z
=
+
g
(2π )3 3(2π )3

≈

Ω20 d

+

p2
d3 p 2 f α ( p )
Ep

Ω2z d
v2
3
3(2π )

(2π )3
Ω20 d
Ω2z d T
+
=
(2π )3 (2π )3 m

2.4.3.2

(2.77)
(2.78)
(2.79)

Ultra-Relativistic Limit

In the ultra-realtivistic limit the temperatures is much higher than the hadron mass:

T≫m

⇒

p
≈1
Ep

(2.80)

This makes the integral very simple:

ξα =

≈

Ω20 gd

(2π )3
Ω20 gd

(2π )3

+

Ω2z
3(2π )3

+

Ω2z gd
3(2π )3

p2
d3 p 2 f α ( p )
Ep

(2.81)
(2.82)

Now we know how to calculate ξ α and thus ∆Nα from the diffused thermal quarks.
This concludes the theoretical chapter. In the next chapter we will apply these results to
implement a Monte Carlo simulation that outputs the generated hadrons, which can be
used to calculate balance functions.

27

Chapter 3
Implementation
This chapter discusses the implementation of the microscopic model derived in the previous chapter. Section 3.1 shows how the lattice quantum chromodynamics data can be
combined with the hydrodynamics data to calculate the quark pair source function and
how quarks are generated from this result. In section 3.2, technical implementation details
of the random walk algorithm and the freezeout detection are discussed. Section 3.3 describes how the quark freezeout events can be translated into hadrons. Section 3.4 gives
a brief description how the hadron cascade was implemented, and finally section 3.5 explains how the balance functions can be calculated.

3.1

Generation of Quark Pairs

To calculate the number of generated quark pairs (2.26), we need to know the charge
fluctuations χ ab . Quantum chromodynamics simulations using lattice gauge theory give
the charge fluctuation over entropy χ/s as a function of temperature T (see Figure 3.1
and [32]). The limits are what we would expect from subsection 2.2.1: For high temperatures we expect a parton gas (2.14); χ should be diagonal and constant in temperature.
Similarly for low temperatures there should be a hadron gas (2.13) dominated by mesons,
because they are lighter than baryons. This explains the negative off-diagonal term. Pions
28

are lighter and thus more common than kaons, explaining that χuu is higher χss .
Unfortunately the lattice QCD data is not available for all the off-diagonal elements
of χ. For low temperatures the densities are very low, so it makes sense to assume a hadron
gas that follows the thermal distribution discussed in subsection 2.4.2. The particle density
is then given by the Jüttner distribution (2.66). χ is then given by Equation 2.13. For high
temperatures we can assume a parton gas and the off-diagonal terms should vanish as
discussed above. Similarly the entropy can be calculated.
Those observations lead to the following approach to approximate χ:
• For diagonal terms, use the lattice data for temperatures above 200 MeV. Linearly
interpolate χ/s as given by the hadron gas at 165 MeV and the lattice data at 200 MeV
otherwise.
• For off-diagonal terms, set χ to zero above 200 MeV. Linearly interpolate χ/s as
given by the hadron gas at 165 MeV and 0 at 200 MeV otherwise.
The lattice data is only available for a finite sample of temperatures (that is, every
5 MeV between 130 and 400 MeV). Values for temperatures in between those samples can
be interpolated using cubic splines.
An Israel-Stewart hydrodynamics code [33] gives the temperature, the entropy and the
velocity as a function of space and time (see Figure 3.2 for an example). The initial conditions chosen when running the code represent the parameters of the collision (centrality,
collision energy, etc.). While the code is capable of running in three dimensions, for our
case we are interested in the relativistic limit, where the hydrodynamics is approximately
two-dimensional due to boost invariance. This means that we only need to consider the
hydrodynamic data in the transverse plane. Furthermore, the code can exploit elliptical
symmetry by only considering one octant of the coordinate system.
Given the temperature, we can now calculate χ using the approach described above.
Plugging this in, the fluid velocity and the entropy into Equation 2.26 gives the source
29

function for thermal pairs of a quark and an anti-quark.1 Assuming a Poisson distribution,
this serves as the number of particles that need to be generated.
A generated quark provides an additional charge to the net neutral quark gluon plasma.
The source function can be negative, in which cases negative weight particles can be
generated, corresponding to quark holes in the plasma. The partner quarks are generated
at the same point in spacetime and will be diffused apart as described in the next section.
For each timestep, each hydro cell and each quark species pair, we are given a positive
real number N representing the number of particles to generate (the source function). The
following algorithm2 is an efficient way to decide how many particles to generate:
1. Set ζ := 0.
2. Uniformly draw a random number r ∈ [0, 1]. Increase ζ by − log r.
3. If N > ζ, then create a quark pair and go to 2. Otherwise halt.
This way the quarks are generated according to a poissonian distribution. Their exact
position can be randomly picked from the corresponding hydro cell (using a uniform
distribution). Now that we have generated the quarks, we can go on to the next step and
diffuse them apart.

3.2

Diffusion of Quarks

As motivated in section 2.3, the thermal quarks are then diffused by performing a random
walk in three dimensions. (Because our hydrodynamics data is two-dimensional, we need
to project to the transverse plane if we need hydrodynamic quantities.) Relativity is respected by making sure that none of the steps are faster than the speed of light. The quarks
1 It is equivalently possible to generate pairs of quarks or pairs of anti-quarks, the weights would have
to be adapted. Consider for instance generating an up and an anti-down quark. Then χ ¯ will be positive
ud
and g + + should be as well. So the weight has to be 1. Analogously, when generating an up and a down
π π
quark the weight has to be -1.
2 See [22] for similar algorithm that scales well for many different species.

30

0.20

0.15
model
lattice QCD

χ s

hadron gas
0.10

quark pair
ss
us
uu

0.05

0.00

200

300

400

T [MeV]
Figure 3.1 Charge fluctuation χ over entropy s as a function of temperature T for various quark pairs. The
dotted vertical line represents the freezeout temperature at 165 MeV. The solid lines represent
the results of lattice QCD simulations of a quark-gluon plasma [32]. The dashed lines show the
results of assuming a hadron gas in thermal equilibrium (2.13). Note that for high temperatures
χ/s is roughly constant according to lattice QCD, and the offdiagonal element (χus ) vanishes
as expected for the parton gas model (2.14). The hadron gas model roughly agrees with the
lattice data for small temperatures. For high temperatures it does not, in that limit we expect the
plasma to behave like a parton gas where all light species are the same, because their masses
are siginifcantly smaller than the temperature (mu ≈ md ≪ ms ≈ 100 MeV). This expectation
is confirmed by the lattice data.

31

Figure 3.2 Temperature T as a function of space ( x, y) in the transverse plane. The solid line represents
the freezeout temperature at 165 MeV. This was taken at the beginning of an Israel-Stewart
hydrodynamics simulation (see [33] for the implementation).

32

are boosted into the restframe of the fluid (whose velocity is given by the hydrodynamics
data), displaced by a random step in a random direction, and boosted back into the lab
frame. See Figure 3.3 for an illustrating plot of such a random walk.
Anytime during this diffusion step, the temperature might fall below freezeout. To
detect this, a line is traced between the starting point and the end point of the diffusion.
The crossed hydro cells can be efficiently computed using a variant of the Bresenham
algorithm [34], which was designed to efficiently render lines on a raster. See Figure 3.4
for an illustration. This yields all the relevant temperatures to consider for freezeout. By
interpolating time linearily, we can check the temperature for all fluid elements of the
trajectory (within the limited resolution of the hydrodynamics data). In the end, this gives
a list of all quark freezeout events.

3.3

Hadronization

Whenever a quark (or quark hole) hits the freezeout surface, we get an additional quark
charge (positive or negative). This can be correlated to the generation of a hadron using
this equation (see subsection 2.4.1):
1 ∆Q
∆Nα = ξ α qα,a χˆ −
b
ab

where
ξα =

g
(2π )3

d3 p

p · Ω 2 − E p /T
e
Ep

and

(3.1)

χˆ ab =

∑ ξ α qα,a qα,b .

(3.2)

α

α represents the hadron species with all the associated properties (mass, charge, degeneracy). In theory, α can be any known hadrons species. In practice, we are considering
319 different species. While most of those are generated in negligible quantities, they are
still considered in the simulation. The change in b-charge ∆Qb is given by the quark freezeout event. For a negative weight quark (a quark hole), ∆Qb has the opposite sign compared

33

trajectories of diffused particles

15
10

d
s
u

y / fm

5
0
5
10
15

15

10

5

0

x / fm

5

10

15

Figure 3.3 Example of a random walk as implemented here. This shows the diffusion trajectories of a
random sample of quarks from the first population until they freeze out. The trajectories are
projected to the transverse plane.

34

Figure 3.4 An illustration of the modified Bresenham algorithm used. Each square cut by the diffusion
trajectory has to be considered, because it is a possibility for a freezeout.

to a positive weight quark. The a-charge qα,a is a property of the hadron species α. The surface element Ωµ must be proportional to the temperature 4-gradient ∂µ T. ∆Nα is invariant
to the normalization of Ωµ (see subsection 2.4.1). When assuming thermal equilibrium, f α
is given by the Jüttner distribution (see subsection 2.4.2), which depends on the mass of
the species and the (freezeout) temperature.
It is not clear at which temperature the freezeout occurs. Plausible values range from
130 MeV to 175 MeV [35], the exact value is essentially a parameter of the model. For
the calculations done here, a freezeout temperature of 165 MeV was chosen. It is worth
mentioning that the cascade can compensate for an earlier freezeout.3
Using all this, ξ α , χˆ ab , and thus ∆Nα can be calculated. Because we can pull the Ω out of
the integral and T is constant, the integral in ξ α needs to be calculated only once for every
1 has to be recalculated for every quark freezeout event, because it depends
species α. χˆ −
ab
3 An example: Choosing a freezeout temperature of 165 MeV yields about as many ρ mesons as protons.
Rho mesons are unstable and decay rather quickly when performing the hadron cascade. If we choose a
freezeout temperature of 130 MeV instead, then there are only few rho mesons compared to protons after
the freezeout.

35

on the hypersurface element Ω. Since χˆ ab is only a 3-by-3 matrix, inverting it is not very
expensive. The resulting ∆Nα can be used to generate hadrons in the same fashion as the
quarks where generated.
Additionally, a momentum from f α has to be drawn, which can be done using rejection
sampling: It is easy to generate a momentum that is exponentially distributed (by taking
the log of a uniformly random number). To honor the remaining coefficient in the integral
we need to find a weight to do rejection sampling. This can be done by factoring out an
arbitrary constant Ωmax , such that we obtain the weight
w :=

2
p·Ω
.
E p Ωmax

(3.3)

To minimize rejection (rejecting often hurts the performance of the program), we should
pick an Ωmax that maximizes w ∈ [0, 1]. For this we need to find the upper bound of the
following expression:
p
p·Ω
Ω = Ω0 −
= Ω0 −
Ep
Ep

m2 + p2

≤ |Ω0 | + |Ω| =: Ωmax
= |Ω · u| +

| p|

( Ω · u )2 − Ω2

|Ω| cos θ

(3.4)
(3.5)
(3.6)

This result can then be used to generate the correct momentum using the following algorithm:
Ep

1. Draw momentum from f ( p) ∼ exp(− T ).
2. Calculate w. Uniformly draw a number r ∈ [0, 1]. If r > w, reject the momentum and
go back to 1. Keep the momentum otherwise.
Now we have determined all the properties of the hadron we need: The momentum
was just calculated, and the position in spacetime is given by the freezeout event. Before
36

being able to calculate the balance functions that correspond to the experimental measurements, the hadrons need to scatter and decay in the cascade.
Efficient Calculation of ξ and χˆ

This paragraph discusses some technical details that

help to avoid calculating ξ and χˆ more often than necessary. While ξ α needs to calculated
for every hypersurface element Ω, it can be separated in parts depending only on α and
parts depending only on Ω:
ξ α =: Ω20 I1,α + Ω2z I2,α

(3.7)

ˆ
Similarly, we can separate χ:

χˆ ab =

∑ ξ α qα,a qα,b = ∑
α

α

Ω20 I1,α + Ω2z I2,α qα,a qα,b =: Ω20 A1,ab + Ω2z A2,ab

(3.8)

The scalars Ii,α and the matrices Ai can be calculated for every species once, as long as
the temperature is constant. Since we are only interested in calculating χˆ at the freezeout
surface, this condition is met.

3.4

Hadron Cascade

The hadrons that froze out after the quarks hit the freezeout surface do not yet correspond
to what one would observe in experiment. To get rid of unstable hadrons (e.g. rho mesons),
it is still necessary to simulate their decays and to scatter them with each other. For this the
b3d code described in reference [20] has been used. After this, the hadrons have reached
their final state and the balance functions can be calculated.

3.5

Calculating Balance Functions

As mentioned in the beginning of chapter 2, there are two contributions to be considered
in order to calculate the balance functions: The contribution due to the generated quark
37

pairs above the freezeout temperature (which has been discussed in this thesis), and the
contribution due to uncorrelated hadrons in the cascade that annihilate each other (which
has been implemented in the b3d code).
For comparison with experimental data, it is imporant to consider the selectivity and
efficiency of the corresponding detector. Because the model yields individual hadrons, it
is straightforward to apply these constraints and include the efficiencies in the calculation.

3.5.1

Contribution above Freezeout Temperature

We are mainly interested in the fate of the quark pair: how it separated and hadronized
and how the resulting hadrons were propagated through the cascade. This gives us the
hadronic balance functions, because we know which hadrons are correlated and how far
they separated by tracing them back to their quark origin.
For each quark a number of hadrons from 0 to ∞ can be generated. Since we expect
∆Nα to be small, this means that we are likely to get a lot of hadrons without a partner,
reducing the amount of useful data. In practice this means that about 80% of the generated
hadrons don’t have a partner.
The propagation of the correlated hadrons through the cascade and the calculation of
the resulting balance function contribution was handled by collaborators.

3.5.2

Contribution below Freezeout Temperature

It is possible to generate uncorrelated hadrons using the momentum distribution of a
relativistic gas (similar to what has been discussed in subsection 2.4.2). These hadrons can
be decayed and scattered as for the hadrons originating in the generated quark pairs from
the previous contribution. Annihilations need to be included, because they contribute to
the charge balance function. Again, the b3d code is used.
After that, the balance functions can be calculated by iterating over all possible hadron

38

pairs, calculating their separation in rapidity and filling histograms corresponding to the
pairs. If a negative weight particle is encountered, the count in the corresponding bin
has to be decreased instead of increased. Recalling the definition of the balance function
allows us to calculate the numerator directly from these histograms:

Bαβ (∆y) ∼ ( Nα − Nα¯ )( Nβ − Nβ¯ ) (∆y)

= Nα Nβ (∆y) − Nα Nβ¯ (∆y) − Nα¯ Nβ (∆y) + Nα¯ Nβ¯ (∆y)

(3.9)
(3.10)

Dividing this by the total number of β and β¯ hadrons then yields the the balance function.

39

Chapter 4
Conclusion
4.1

Summary

We have developed a model that will provide quantitative predictions of the number of
quark-antiquark pairs generated during the evolution of a relativistic heavy ion collision.
The production rate depends on results of lattice gauge theory for charge susceptibilities
and on the hydrodynamic evolution which provides the entropy density and temperature
profiles. The separation of the correlated quarks was modeled using a diffusive process
by employing a random walk for each individual quark. Care was taken to ensure that
causality was not violated in the description. These correlations were propagated through
a hadronization formula across an interface between the hydrodynamic description and
a hadronic cascade modeling the breakup and dissolution of the fireball. The correlations
from the initial charge conservation were tracked through the cascade and summed alongside correlations from decays and annihilations during the cascade. The final-state correlations were quantified in the form of charge balance functions, which can be compared
to experiment.
Once the approach is better tested and the responses of the experimental detectors are
accurately incorporated, the results from this approach should give profound insight into

40

the chemical evolution of the quark-gluon plasma. Charge balance functions have been
shown to be sensitive to both the chemical make-up of the quark-gluon plasma and the
diffusion. We hope that this approach will provide quantitative answers to the question
of whether the matter formed in relativistic heavy ion collisions had the same numbers
of up, down and strange quarks as the equilibrated matter studied in lattice gauge theory.
Further, one may even be able to extract the diffusion constant of the quark-gluon plasma,
a fundamental transport coefficient that has thus far not been addressable either through
experiment or by lattice calculations.

4.2

Future Work

The model shown here represents the field’s first attempt to microscopically model the
evolution of charge correlations. Most of the pieces of the approach are now in place, but
significant work is required to complete the model and to validate that it indeed faithfully
represents the equations described here. Once finished, the model will address a variety of
measurements in a way that provides insights into the chemical evolution of a heavy-ion
collision. A list of short-term goals is presented here:

Testing the Model

The model has not yet been extensively tested for self-consistency.

For instance, it should be ensured that the balance functions for electric charge and baryon
number integrate to unity.

Comparing with Experiment A detailed study how well the predictions agree with experimental measurements is necessary. For this, the experimental acceptances need to be
modeled in detail. The results of the model should be compared to the existing measurements of π + π − , K + K − , p p¯ and pK balance functions. Furthermore, balance functions
binned in relative azimuthal angle and by transverse momentum can be compared to
experimental results.
41

Diffusion Constant By simultaneously studying balance functions in relative rapidity
and azimuthal separation, it might be possible to distinguish the separation caused by
diffusions from the separation caused by string breaking or color-flux tube tunneling at
early times.

Freezeout Temperature The exact value of the freezeout temperature was chosen arbitrarily (within a plausible range). A parameter study could explore the relation between
the choice of this temperature and the resulting balance functions. The results are expected
to be sensitive to any change that alters the final-state yields.

Non-Central Collisions

By changing the initial condtitions of the hydrodynamics code,

centrality as a parameter could be studied. As discussed in the introduction, the balance
functions are expected to become broader for increasing impact parameter. The extent of
this effect could be tested.

More Lattice QCD Data

Currently, lattice quantum chromodynamics data is missing

for the ud off-diagonal terms in the charge susceptibility. This data would allow avoiding
to interpolate the values of the susceptibility in the region near the critical temperature.

42

APPENDICES

43

Appendix A: Notation and Conventions
Units The following units are used:

c = h¯ = k b = 1

(4.1)

Special Functions The Dirac δ-function is defined by the integral equation

f ( x )δ( x ) dx := f (0) .

(4.2)

The Heaviside step function can be given as

θ ( x ) :=

x

−∞

δ(y) dy =

(The value at 0 is an arbitrary convention.)





0





1

2






1

x<0
x=0 .

(4.3)

x>0

The gamma function is defined as
∞

Γ(t) :=

0

x t−1 e− x dx .

(4.4)

For t ∈ N this simplifies to
Γ ( t ) = ( n − 1) ! .

44

(4.5)

Minkowski Space 4-vectors are typeset bold, 3-vectors have an arrow:


x0





 
 
 x1 
 
x= 
 x2 
 
 
x3

x1



 
 
2
x=
x 
 
x3

(4.6)

Latin indices are spatial, greek indices can be temporal:

i ∈ {1, 2, 3}

µ ∈ {0, 1, 2, 3}

(4.7)

If a term contains the same index two times, it implies an implicit summation over this
index (Einstein sum convention).
The Minkowski metric is assumed to be




1

 −1

η=

−1



45

−1





 .




(4.8)

Appendix B: Normalization of Maxwell-Jüttner Distribution
We want to calculate the integral

d :=

d3 p exp −

u·p
T

.

(4.9)

d is Lorentz-invariant, so we can choose a reference frame where u = (1, 0, 0, 0):

d=

p2 dp d(cos θ ) dφ exp −

Ep
T

∞

= 4π

0

p2 e

Ep
− T
dp

(4.10)

Choosing the Lorentz factor as a new variable of integration
Ep
p 2
= 1+
m
m
p
dγ =
dp
mE p

⇔

γ:=

γ2 − 1

p=m

⇔

dp =

mγ
γ2 − 1

dγ

(4.11)
(4.12)

gives
d = 4πm2

∞
1

mγ

m

− γ
γ2 − 1e T dγ .
3

Now we can integrate by parts using ∂γ (γ2 − 1) 2 = 3γ
d = 4πm2

= 4πm2

(4.13)

γ2 − 1:

3 −mγ
m2 2
(γ − 1) 2 e T dγ
3T
m
T2
m
= 4πm2 TK2
· 3 2 K2
T
T
m

∞
1
m2

3T

giving the normalization.

46

(4.14)
,

(4.15)

BIBLIOGRAPHY

47

BIBILIOGRAPHY
[1] S. A. Bass, P. Danielewicz, and S. Pratt. “Clocking Hadronization in Relativistic
Heavy-Ion Collisions with Balance Functions”. In: Phys. Rev. Lett. 85 (13 Sept. 2000),
pp. 2689–2692. D O I : 10.1103/PhysRevLett.85.2689.
[2] S. Pratt, C. Ratti, and W. P. McCormack. “Chemical properties of super-hadronic
matter created in relativistic heavy ion collisions”. In: ArXiv e-prints (Sept. 2014).
arXiv: 1409.2164 [nucl-th].
[3] D. Drijard et al. “Density, charge and transverse √
momentum correlations of particles
in non-diffractive proton-proton collissions at s = 52.5 GeV”. In: Nuclear Physics
B 155.2 (1979), pp. 269–294. I S S N : 0550-3213. D O I : http://dx.doi.org/10.1016/
0550-3213(79)90269-4.
[4] D. Drijard et al. “Quantum number effects in events with a charged particle of large
transverse momentum: (II). Charge correlations in jets”. In: Nuclear Physics B 166.2
(1980), pp. 233–242. I S S N : 0550-3213. D O I : http://dx.doi.org/10.1016/05503213(80)90226-6.
+ −
[5] H.
√ Aihara et al. “Observation of Strangeness Correlations in e e Annihilation at
s = 29 GeV”. In: Phys. Rev. Lett. 53 (23 Dec. 1984), pp. 2199–2202. D O I : 10.1103/
PhysRevLett.53.2199.

[6] J. Adams et al. “Narrowing of the Balance Function with Centrality in Au + Au
√
Collisions at s NN = 130 GeV”. In: Phys. Rev. Lett. 90 (17 May 2003), p. 172301.
D O I : 10.1103/PhysRevLett.90.172301.
[7] STAR Collaboration et al. “Longitudinal scaling property of the charge balance
function in Au+Au collisions at s=200 GeV”. In: Physics Letters B 690 (June 2010),
pp. 239–244. D O I : 10.1016/j.physletb.2010.05.028. arXiv: 1002.1641 [nucl-ex].
[8] M. M. Aggarwal et al. “Balance functions from Au+Au, d+Au, and p+p collisions at
s NN =200 GeV”. In: Phys. Rev. C 82.2, 024905 (Aug. 2010), p. 024905. D O I : 10.1103/
PhysRevC.82.024905. arXiv: 1005.2307 [nucl-ex].
for charged
[9] G. Westfall and the STAR Collaboration. “Widths of the balance function
√
√
kaon pairs from Au + Au at s NN = 200 GeV and p + p at s = 200 GeV at
RHIC”. In: Journal of Physics G: Nuclear and Particle Physics 30.1 (2004), S345. D O I :
10.1088/0954-3899/30/1/040.
[10] H. Wang. “Study of Particle Ratio Fluctuations and Charge Balance Functions at
RHIC”. In: ArXiv e-prints (Apr. 2013). arXiv: 1304.2073 [nucl-ex].
[11] X.-N. Wang and M. Gyulassy. “HIJING: A Monte Carlo model for multiple jet production in pp, pA, and AA collisions”. In: Phys. Rev. D 44 (11 Dec. 1991), pp. 3501–
3516. D O I : 10. 1103 /PhysRevD.44 .3501. U R L : http:/ /link.aps .org/doi /10.
1103/PhysRevD.44.3501.

48

[12] S. A. Bass et al. “Microscopic models for ultrarelativistic heavy ion collisions”. In:
Progress in Particle and Nuclear Physics 41 (1998), pp. 255–369. D O I : 10.1016/S01466410(98)00058-1. eprint: nucl-th/9803035.
[13] S. Cheng et al. “Statistical and dynamic models of charge balance functions”. In:
Phys. Rev. C 69.5, 054906 (May 2004), p. 054906. D O I : 10.1103/PhysRevC.69.054906.
eprint: nucl-th/0401008.
[14] C. Alt et al. “System size and centrality dependence of the balance function in A + A
√
collisions at sNN = 17.2 GeV”. In: Phys. Rev. C 71 (3 Mar. 2005), p. 034903. D O I :
10.1103/PhysRevC.71.034903.
[15] A. Bialas. “Balance functions in coalescence model”. In: Physics Letters B 579 (Jan.
2004), pp. 31–38. D O I : 10.1016/j.physletb.2003.10.106. eprint: hep-ph/0308245.
˙
[16] P. Bozek,
W. Broniowski, and W. Florkowski. “Balance functions in a thermal model
with resonances”. In: Acta Phys.Hung. A22 (2005), pp. 149–157. D O I : 10.1556/APH.
22.2005.1-2.15. arXiv: nucl-th/0310062 [nucl-th].
˙
“The balance function in azimuthal angle is a measure of the transverse
[17] P. Bozek.
flow”. In: Physics Letters B 609 (Mar. 2005), pp. 247–251. D O I : 10.1016/j.physletb.
2005.01.072. eprint: nucl-th/0412076.
[18] B. Ling, T. Springer, and M. Stephanov. “Hydrodynamics of charge fluctuations
and balance functions”. In: Phys. Rev. C 89.6, 064901 (June 2014), p. 064901. D O I :
10.1103/PhysRevC.89.064901. arXiv: 1310.6036 [nucl-th].
[19] G. Agakishiev et al. “Anomalous centrality evolution of two-particle angular correlations from Au-Au collisions at s NN =62 and 200 GeV”. In: Phys. Rev. C 86.6, 064902
(Dec. 2012), p. 064902. D O I : 10 . 1103 / PhysRevC . 86 . 064902. arXiv: 1109 . 4380
[nucl-ex].
[20] S. Pratt and J. Vredevoogd. “Femtoscopy in relativistic heavy ion collisions and its
relation to bulk properties of QCD matter”. In: Phys. Rev. C 78.5, 054906 (Nov. 2008),
p. 054906. D O I : 10.1103/PhysRevC.78.054906. arXiv: 0809.0516 [nucl-th].
[21] J. D. Bjorken. “Highly relativistic nucleus-nucleus collisions: The central rapidity
region”. In: Phys. Rev. D 27 (1 Jan. 1983), pp. 140–151. D O I : 10.1103/PhysRevD.27.
140.
[22] S. Pratt. “Accounting for backflow in hydrodynamic-Boltzmann interfaces”. In:
Phys. Rev. C 89.2, 024910 (Feb. 2014), p. 024910. D O I : 10.1103/PhysRevC.89.024910.
arXiv: 1401.0316 [nucl-th].
[23] F. Cooper and G. Frye. “Single-particle distribution in the hydrodynamic and statistical thermodynamic models of multiparticle production”. In: Phys. Rev. D 10 (July
1974), pp. 186–189. D O I : 10.1103/PhysRevD.10.186.
[24] S. Pratt. “General charge balance functions: A tool for studying the chemical evolution of the quark-gluon plasma”. In: Phys. Rev. C 85.1, 014904 (Jan. 2012), p. 014904.
D O I : 10.1103/PhysRevC.85.014904. arXiv: 1109.3647 [nucl-th].

49

[25] J. Dunkel, P. Talkner, and P. Hänggi. “Relativistic diffusion processes and random
walk models”. In: Phys. Rev. D 75.4, 043001 (Feb. 2007). D O I : 10.1103/PhysRevD.75.
043001. eprint: cond-mat/0608023.
[26] F. Jüttner. “Das Maxwellsche Gesetz der Geschwindigkeitsverteilung in der Relativtheorie”. In: Annalen der Physik 339.5 (1911), pp. 856–882. I S S N : 1521-3889. D O I :
10.1002/andp.19113390503. U R L : http://dx.doi.org/10.1002/andp.19113390503.
[27] G. M. Kremer. Theory and applications of the relativistic Boltzmann equation. Feb. 2013.
U R L : http : / / www . ift . uni . wroc . pl / ~karp49 / LadekLectures2013 / Monday /
gmKremer1_49th.pdf.
[28] E. Weisstein. Modified Bessel Function of the Second Kind. May 2014. U R L : http://
mathworld.wolfram.com/ModifiedBesselFunctionoftheSecondKind.html.
[29] Wolfram Research, Inc. Mathematica 9.0. 2012.
mathematica/.

U R L:

http : / / www . wolfram . com /

[30] SymPy Development Team. SymPy: Python library for symbolic mathematics. 2014.
U R L : http://www.sympy.org.
[31] Wikipedia. Meijer G-function — Wikipedia, The Free Encyclopedia. [Online; accessed
5-June-2014]. 2014. U R L : https://en.wikipedia.org/wiki/Meijer_G-function.
[32] S. Borsányi et al. “Fluctuations of conserved charges at finite temperature from
lattice QCD”. In: Journal of High Energy Physics 1, 138 (Jan. 2012), p. 138. D O I : 10.
1007/JHEP01(2012)138. arXiv: 1112.4416 [hep-lat].
[33] J. Vredevoogd. “Relativistic Viscous Hydrodynamics for High Energy Heavy Ion
Collisions”. In: ArXiv e-prints (July 2013). arXiv: 1307.7677 [nucl-th].
[34] Wikipedia. Bresenham’s line algorithm — Wikipedia, The Free Encyclopedia. [Online;
accessed 23-June-2014]. 2014. U R L : https://en.wikipedia.org/wiki/Bresenham’
s_line_algorithm.
[35] P. Braun-Munzinger and J. Stachel. “The quest for the quark-gluon plasma”. In:
Nature 448 (July 2007), pp. 302–309. D O I : 10.1038/nature06080.

50