_. .: $2 a" r ... t4 . 4.7;: I ”V, 2r {WNW ». ‘u it. A .3. Ann “ z 3 a 1.. Y. ~ 4?? “WWW.“ . , v “In-ll :3 a. a. n. 4.3... 1A.... .3 ., La. .. S alien :fi Ln. 4...! .227: .3... .l .Rfi $214.93...) awn. rt ruin... rinsmum T‘s: "Am: 1,. b . .39: ‘1' My. 15‘. LlBRr‘f-‘tRY Michigan State University This is to certify that the thesis entitled Statistical Correlations in Relativistic Heavy Ion Collisions presented by Silvio Petriconi has been accepted towards fulfillment of the requirements for the MS. degree in Physics / r , Q Majoriyfofess‘or’s Signature 08/20/03 Date MSUls anAflhnatIve AdorVEqualOppoMlylnsflluflon 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 c:/CIFIC/DateDuo.p65~p. 15 STATISTICAL CORRELATIONS IN RELATIVISTIC HEAVY ION COLLISIONS By Silvio Petriconi A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Physics and Astronomy 2003 Thesis Adviser: Dr. Scott Pratt ABSTRACT STATISTICAL CORRELATIONS IN RELATIVISTIC HEAVY ION COLLISIONS By Silvio Petriconi In order to understand the state of matter produced in relativistic heavy ion collisions, it is crucial to understand the space-time characteristics of the emitting source. Whereas the momenta of particles can be directly measured, space-time information can only be inferred from two-particle correlation observables. In this thesis, two different kind of correlation observables will be discussed: Part I will investigate the two-particle correlation function for non-identical particles. It will be shown that the correlator is affected by Coulomb and strong final-state interactions in a way which allows extracting the size and lifetime of the source. ‘ Part II will discuss the modeling of particle multiplicities for charge balance functions. The thermal model of a canonical ensemble under charge conservation constraints will be used to generate particle multiplicities at various temperatures. A highly optimized algo- rithm will be developed for calculating the partition functions iteratively. Balance functions obtained using this thermal model and the blast-wave description will be compared to ex- perimental data, and good agreement will be seen. To my parents iii ACKNOWLEDGMENTS This thesis would not have been possible without the strong support of many people and organizations. I am deeply indebted to my adviser, Dr. Scott Pratt, for his strong academic support, his patience and his financial assistance. I wish to thank Dr. Wolfgang Bauer for giving me the opportunity of spending an enriching and interesting year of studies at MSU. I gratefully acknowledge the support of the Studienstiftung des deutschen Volkes which gave me orientation and motivation during the entire period of my studies. There have been many people who have been giving me understanding and moral sup- port whenever I needed it. Special thanks are due to my family, and to my friends Andreas Glaser and Kristina Wagner for giving me encouragement and motivation in some critical moments . iv Table of Contents List of Tables List of Figures 1 Introduction I Final-state Correlations of non-identical Particles 2 Formalism of Correlation Functions 2. 1 Introduction .................................. 2.2 The general FSI correlator .......................... 3 The Quantum Correlator 3.1 Correlations from Coulomb interaction .................... 3.2 Correlations from strong interactions ..................... 3.2.1 Basic partial wave analysis ...................... 3.2.2 Handling spin flipping ........................ 3.2.3 Isospin rotation under a resonance .................. 3.2.4 Constructing the wavefimction inside of the interaction range . . . . 4 The Classical Coulomb Correlator 4.1 Derivation of the classical correlation weight ................ 4.2 Numerical Evaluation ............................. 5 Numerical Results and Discussion 5.1 Choosing particle pairs for correlation analyses ............... 5.2 12K“ correlations ............................... 5.3 p7r+ correlations ................................ 5.4 E‘rr+ correlations ............................... 6 Conclusion and Outlook II Modeling Balance Functions in Heavy Ion Collisions 7 The Balance Function Formalism 7. 1 Introduction .................................. 7.2 Balance functions as a QGP signal ...................... 8 Modeling Particle Multiplicities 8.1 Recursion relations for partition functions .................. V vii viii 5 5 7 l4 l4 l6 l7 I9 20 22 26 26 31 35 35 37 39 42 43 45 46 47 48 52 52 8.2 Optimization issues .............................. 54 8.3 Monte-Carlo generation of particles ..................... 56 9 Results and Conclusion 60 9.1 Constructing the balance function ....... ~ ............... 60 9.2 Conclusion .................................. 61 APPENDICES 63 A Recursion Relations for Coulomb Waves 64 A.l Transformation to natural coordinates .................... 64 A2 Evaluation of Eq. (3.36) ........................... 65 A21 Case I: L sé 0 ............................ 66 A.2.2 Case 2: L = 0 ............................ 68 B Particles used in the Thermal Model 69 BIBLIOGRAPHY 73 vi List of Tables Particles used in the thermal model ........................... 69 vii 4.1 5.1 5.2 5.3 5.4 5.5 7.1 7.2 7.3 9.1 List of Figures Rotation of the source frame for Monte Carlo integration ............ pK+ correlations for a Gaussian source; Rlong = Rside = 4 fm, Rout = 8 fm Angular pK+ correlation; Rlong == Rside = 4 fm, Rout = 8 fm ......... Angular dependence of the p7r+ correlator near the A++ resonance; Rlong = RSIdC = 4 fm, Rout = 8 fm .......................... Shadowing by strong interaction near the resonance momentum; Rlong = Rside : R0“! = 4 firl; Zofiset : 40 fl“ .................... E'7r+ correlations at resonant momentum; Rlong = Rside = 4 frn, Rout = 8 fm Balance functions obtained from hadronic cascade models ........... Narrowing of the pion balance function with increasing centrality of the colli- sion as measured at STAR .......................... Balance functions for charged pions as measured at STAR in. comparison to theoretical data from hadronic cascade models ................ Thermally modeled balance function of identified pions in comparison to STAR measurements ............................. viii 32 38 39 40 41 42 51 Chapter 1 Introduction When speaking of the nucleus of an atom, we are used to thinking of it as a quantum me- chanical composite object consisting of finite-size particles, protons and neutrons. How- ever, one should realize that this picture is only a good description of a special phenomenol- ogy that we observe at lower energies. Neutrons and protons are not true elementary par- ticles as we know from early high-energy physics experiments. Instead, they show the behaviour of a composite system consisting of three quarks which interact strongly with each other by the exchange of gluons. Unlike any other composite system in physics, the constituents of protons, neutrons or any other bound quark system cannot be observed in- dividually as any attempt of separation is immediately followed by the creation of a new composite system. This quark confinement phenomenon comes from the special long-range nature of the strong force and its explanation requires non-perturbative models. Lattice QCD provides a numerical treatment from first principles by evaluating the partition function through 1 randomly sampling the many-body quantum state. By performing these calculations at finite temperature, lattice QCD can give the equation of state and phase structure of the QCD vacuum. These calculations suggest that a deconfined state of quarks and gluons may exist at temperatures above 150 — 200 MeV, or at energy densities of 1 — 3 GeV/fm3, respectively [33, 34]. At these temperatures, a nucleus would undergo a phase transition and form a quark-gluon-plasma (QGP) in which the major degrees of freedom are quarks and gluons rather than neutrons and protons. However, as this result is a purely theoretical prediction, it is now of significant interest to obtain experimental confirmation. The latest generation of relativistic heavy ion colliders was designed to operate clearly above the critical energy density for the formation of a QGP in order to search for experi- mental evidence for this phenomenon. However, since such a quark-gluon plasma interme- diate state is short-lived and quarks are reconfined after a time as short as a few frn/c, the most challenging part of the search for the quark-gluon plasma in heavy ion collisions is to find observables that can give a clear indication whether or not an intermediate QGP state was formed during the collision. Many such observables (“QGP signals”) have been suggested over the last two decades [15]: J / ti) suppression due to QGP screening of the c — 6 potential [42, 17] Strangeness enhancement [36, 53] Charge correlations [13, 4] Dileptons [2, 59, 8, 55] . Flow [22, 21, 57] 0 Jet quenching [27, 25, ll, 10, 12, 3, 5, 7] o HBT measurements of lifetime [44, 56, 54] In the first chapters of this thesis, we will investigate the last topic concerning statistical correlations in momentum space of the particles that are emitted in the collision process. Whilst identical boson correlations have been used exhaustively in Hanbury-Brown Twiss (HBT) interferometry to measure the size and lifetime of sources from which these parti- cles were emitted [48, 61], relatively little attention has been paid to correlations of non- identical particles [41, 58, 39]. These correlations are caused by final-state Coulomb and strong interactions rather than through identical-particle interference. It will be shown that such correlations can give important information on source parameters especially those that are associated with the lifetime of the source. These are especially relevant since a phase transition, and the associated latent heat, should lead to a reduced pressure and an extended lifetime for emission. The second part of this thesis will discuss charge balance functions as suggested in [13]. These fitnctions exploit the constraint of local charge conservation in a heavy ion collision and are a promising candidate for a QGP signal. From analysis of relative rapidities of par- ticles and their balancing partners, one can infer the timing of hadronization. In this thesis, we will develop an implementation of a thermal model for generating particle multiplicities under charge conservation constraints. The predictions of the thermal model will finally be compared against the shape of the balance function that has been measured at STAR [4]. Part I Final-state Correlations of non-identical Particles Chapter 2 Formalism of Correlation Functions 2.1 Introduction Statistical correlations in momentum space of particle pairs were initially used in astron- omy by Hanbury-Brown and Twiss to measure the size of distant stellar objects [2 8]. Later, Hanbury-Brown Twiss (HBT) interferometry also became a standard method in heavy ion physics [44, 29, 61] for investigating the size and lifetime of the emitting source, i.e. the distribution of space-time coordinates for which the final-state measured particles experi- enced their last interaction. The basic idea of HBT measurements is that particles emitted by a chaotic source are initially statistically uncorrelated. On their way from the source to the detector, however, identical bosons are afi‘ected by quantum mechanical Bose-Einstein symmetrization. This results in an enhancement in the coincidence probability for two par- ticles with similar momenta. By fitting the measured coincidence probabilities of identical bosons to model-generated data, source sizes can be estimated. 5 In relativistic heavy ion physics, source sizes provide important information on the physics of the intermediate state. If a quark-gluon plasma is formed and the phase transition involves a significant amount of latent heat, it is expected that the lifetime of the emitting source would increase significantly [44, 54]. If the emission probability for a particle with velocity v at spacetime point (t, x) is S (t, x, v), the correlation fimction is determined through the object _ fd4xld4x281(x1,v)82(z2,v)63(r - (x1 - x2 — v(t1 - t2))) g(r,v) = fd4zrlfd42251(:131,v)52($2,v) (2.1) which defines the probability that particles of type A and B with the same velocity v would be separated by r after emission. A long-lived source would lead to a shape of g which is significantly longer in the direction of v, also called the “outwards” direction. For rela- tivistic heavy ion physics it is the convention to define v in a frame moving along the beam axis such that the component of p 2 pa + p, along the beam is zero. Thus, v is manifestly perpendicular to the beam axis. The shape of g in the other two directions — the “longi- tudinal” direction along the beam axis and the “sidewards” direction that is perpendicular to the two others would not be affected by the lifetime. Correlation analyses measure the characteristic radii of g in these three directions. Assuming that S (:12, t) has comparable outward and sideward sizes in the lab frame, its lifetime can then be estimated by 1372 z Rama — 212/8) — 122 (2.2) side where Rom and Rside are the experimentally inferred dimensions of g in the two-particle 6 rest frame. Remarkably, HBT correlation analyses at RHIC have shown nearly equal outward and sideward source sizes [6] which indicates extremely sudden emission. This strongly con- tradicts theoretical models involving a first-order phase transition. As this result comes quite unexpected, alternative methods for confirming these findings were explored. One of these methods is the use of non-identical particle correlations that will be discussed here. It was noticed early in HBT interferometry that aside from Bose-Einstein symmetriza- tion, significant additional correlations arise due to Coulomb and strong final-state inter- actions between the particle pair. These additional correlations have been investigated ex- haustively [9, 60, 14] with the aim of correcting for them. Non-identical particles do not have any (anti-)symmetrization of their final states. Their correlations are dominated by Coulomb and strong FSI effects which for a long time were considered to be of little value in measuring the shape of g(r). In the following chapters, it will be shown that strong and Coulomb induced correlations provide equivalent information for the lifetime and spatial extent of the source and can be used to independently verify HBT results. 2.2 The general FSI correlator The two-particle correlation fimction is a measure for the enhancement in probability P2(pa, pb) to detect two particles from the same event at momenta pa and pi» so if P1(pa) and P1(pb) are used to designate the corresponding one-particle probabilities, one defines the correlation function as P2(pa, pb) P1(Pa)P1(Pb) 7 C(Pa, 1):») E (23) For non-identical particles, we will now show that this enhancement is approximately given by folding the square of the nonrelativistic wavefimction that solves the Schrédinger equation for the relative motion with the source emission function g(r), C(qw) = [(1371) g(ro,\t)|u”(l'o,<1)l2 (14) Here, q is the relative momentum of the two-particles in their center-of-mass frame. For particles without interaction, the wavefunction would be a plane wave and the correlator would be unity since g(r) is normalized to unity. Interactions between the particles distort the wavefunction and give higher weight to certain regions of relative momentum which results in the observable structure of the correlator. The derivation starts with considering the quantum-mechanical two-particle T -matrix element T(.r1,;1:2): If one could draw a complete Feynman diagram in coordinate space of the entire collision process, one could factorize it into an horribly complicated part T where all particles interact in arbitrary ways, and into another well-understood part U in which only the two particles that are finally detected at t —+ 00 with momenta pa and pb continue interacting with each other. We will label the two vertices where the diagram is factorized 1'1 and .272. The total two-particle probability is the sum of the squared amplitudes of all dia- grams which have a such a factorisable two-particle part with an asymptotic outgoing state 8 Ipm pb), so 2 P2(pavpb) : Z i/firld4x2Z(irla$2)U($la$2apaapb) (25) f Just as in HBT interferometry [45], one can assume fully chaotic emission behaviour of the source. This is equivalent to requiring that the two-particle T-matrix should simply be the product of two one-particle T-matrices: ’1',(.2:1, $2) *5 7}a(1’1)ij(-732) (2-6) It is technically convenient to change from the complex-valued amplitude T into a Wigner density S(.r, K) defined by 1 . , Sir, K) 2 Z (27,). / d4y 8mm, + gm — g) (2.7) Wigner densities are real but not necessarily positive definite. They contain all phase-space emission characteristics of the source and can be used like a probability density for particle emission at point a: with momentum K as long as one of these two parameters is integrated over. The probability for one-particle emission with momentum p therefore is mp) = fan: SW) (2.8) P0=EK The two-particle emission probability can also be expressed in terms of the Wigner density S. One can separate the asymptotic center-of-mass motion of the two particles ma 171 + mng ) u,(;r1 — 1'2, k) (2.9) U(l‘lil"27paipb) E exp (2K ma + mb where K = pa + 1);, and k = (mbpa — Tnapb)/(m.a + mb) are the total and the relative four-momenta of the two particles. Note that this separation is essentially non-relativistic. For relativistic particle pairs it still holds as long as the relative particle momenta after a boost to the center-of-mass frame are non-relativistic. The particle pairs investigated in correlation analyses have usually sufficiently low momenta to satisfy this requirement. The next step is to define the Wigner transformed Wk of the relative FSI interaction matrix element I'Vk(q, :17) E fd‘y eiq'yu'(.r + y/2, k)u(;r: — y/2, k) (2.10) (2W By performing Fourier transforms of 5', one can re-express the T- and U -matrix elements in terms of the corresponding Wigner functions. Equation (2.5) then reads , maK P(paapb) = [d1qd4231d41'2 Ls ((131, —— + q) ma + mb . K x 8(12—1’L— —- q) W.(q,a:2 — x.) (2.11) ma + mb The densities S have a clear physical interpretation as emission @seudo-)probability den- sities of the source. If the product of S 1 and 52 is largely independent of q, one can make the so-called smoothness approximation and neglect the entire q dependence in 31 SQ. For 10 systems with a broad 7' matrix element this assumption is safe, and for instantaneous ther- mal models it is even exact since the product of the two functions gives oc exp(E1 + E2) which is independent of q. If the q dependence can be entirely neglected, the dq integration gives a delta function in y, and one obtains P(Pa,Pb) = fdd’l‘r ($11725 (1131, aft—:%) 5 (1172» $5752) ”*($2-$1ak)U($2—$1ak) (2.12) If the emission of the two particles would be simultaneous, further simplification would be possible since the u‘u matrix element would then become the wavefunction of the corre- sponding field. After boosting into the K = 0 frame, a." = A(K):c, one can make the ansatz that the particle that was emitted earlier will propagate along the classical trajectory of free motion with velocity v until the second particle is emitted. The matrix element u‘u is then taken at equal time t’ = t’z. This gives the wavefimction [112(k’, x’2 — x’1 — v(t’2 — t’1))|2 of the underlying field. Even though the fields of pions, protons, kaons and 3 particles are relativistic spin 0 and spin 1/2 fields, the relative motion in the center-of-mass frame can be qualitatively described by a Schriidinger field as long as the relative momenta in the K = 0 frame are in the non-relativistic regime. The corresponding wavefunction will then be a solution of the Schtidinger equation that describes the relative motion of the particle pair. As the Schrodinger equation is well-investigated for a large class of problems, the charm of this method is that all the wealth of information on how to treat such interactions in nonrelativistic quantum mechanics becomes available for use in final-state correlation analyses. 11 The two-particle probability in the c.o.m. frame finally reads Pam = f .14.; f .14.; / d31‘05(2"1,v)S(1r.'2,v)63(ro—x'1—x'2—v(t’l—t’2))) wean? (2.13) Here, q is the relative momentum in the c.o.m. frame. Performing the d41t1 and (141:2 integration and dividing by the corresponding one- particle probabilities (2.8) one obtains for the final-state correlator the result We» = f can. g(v,1‘o)lit/’(romr)!2 (2.14) I did I da.:_,s(x’1, v)S(:1:’2, v)63(r0 _ 361 _. x12 _ v(t’1 _ t’2) fdfi18($a,V)fd%&S($é’v) (2.15) g(VJ‘O) This form of the correlator was given first by Koonin [3 8] in 1977 and has since then been re-derived in various ways [9, 45]. Note that g(v, r0) contains the emission Wigner density 8 only in a time integrated form. Unique extraction of S by correlation measurements is therefore not possible. For- tunately, for using correlations as a QGP signal this is also not necessary since delayed freeze-out can be detected as an enhancement of g in the outwards direction due to the source propagation properties that we discussed in the previous section. For practical correlation analyses, one has to choose a model for g since the true emis- sion Wigner densities S are unknown. A common choice is to use a Gaussian that has an enhancement in its outwards size, 2 2 1 rgut rside + Tlong 'v,r = .. , ex —————‘——'—' 2'16 g( 0) (47r)'3/ZRourRi p ( 4Rgut(v) 4R1“) ( ) 12 Correlations from final-state Coulomb and strong interactions can then be simulated by folding this source function with the appropriate wavefunction. We will demonstrate this for Coulomb, classical Coulomb and strong interaction in the following chapters. 13 Chapter 3 The Quantum Correlator Treating Coulomb and strong final-state interactions is straightforward once the wavefunc- tion formalism (2.14) has been established. One just needs to solve the Schrodinger equa- tion for the relative motion of the particle pair and choose a model for the source func- tion. For the Coulomb case, solutions will be obtained in parabolic coordinates, whilst the partial-wave analysis for including strong interactions will be performed in spherical coordinates. 3.1 Correlations from Coulomb interaction The effects of Coulomb interaction can be calculated quite easily as it is possible to con- struct an analytic solution for the continuous energy spectrum of the Schréidinger equation with Coulomb interaction, 2) 212(r) = EMI‘) (3-1) This is the equation for the reduced system that is obtained after separating center-of-mass and relative motion using the the reduced mass parameter II = 77117722 / (m1 + 77l2). The next step is to choose a coordinate system in which separation of solutions becomes possible. Clearly, given the spherical symmetry of the problem, a separation approach in spherical coordinates seems to be most promising. Indeed, separation in spherical co- ordinates is possible and leads to an infinite series of partial Coulomb waves which are simultaneous eigenstates of energy and angular momentum. This partial wave decompo- sition will be useful later; for constructing the total Coulomb wavefimction numerically, however, separation in parabolic coordinates { = r + z, r} = r — z, tp = arctan (2) (3.2) ,1? is more convenient as it immediately leads to a closed solution. We will just briefly review the solution as it can be found in quantum textbooks [40, 43]. The inverse transformation for parabolic coordinates is found to be 1 ”flaw. y=\/E—Usin = \g l1,0)|T) + \/g I1,1>|i> (3-18) $11 $10 m — 3|v>|l>+ Eli) ( ) l9 NIH [\DIi—I H V II Solving for the initial state gives 2 3 1 1 l 1 After phase-shifting the J = 3/ 2 part by 62‘6"J=3/2 and the J = 1/2 part by 62‘6'-J=1/2 one can use equations (3.18) and (3.19) to write (3.20) again in terms of the spherical harmonics of the partial wave states: 2 _ 1 , 2:0 216 _« 216 __ 6 IH=I.m=0 M (38 I'J—J/z + 36 ("l-U2 l=l,m=0 IT) 2 . 2 - + (_\:/3_:€216,.J=3/2 _ éewéua/z) Yl=l.m=l ll) (32]) 3.2.3 Isospin rotation under a resonance If two particles interact through a resonance, they can convert into different species as only the sum of the particles’ isospin is conserved. As an example, the two-particle product state of E’ and a 7r+ can go through the 30* resonance and give either a E’7r+ or a 3%“ pair. This effect can be calculated quantitatively by simply making an isospin decomposition of the incoming state. The E‘vr contributes to two different isospin states, which are l1=—m,= —>= \/:|:-7r +)+\/;|— 30°) (3.22) and its orthogonal state, l1=§ml=2 —=> —\:}|E_7r +>+\/;2;|EO7r 0) (3.23) From this immediately follows ._._ 1 .—.- \/§ _ I- 7r+> _ (5|- 71.+> | 3 I‘07r0>) 1:3/2 2 ,__ + x/i .. + (3 I: 7r )_ 3— |:07r0>) (3.24) I=l/2 In the vicinity of the 30‘ resonance, the isospin dependent phase-shift for I = 1 / 2; J = 3 / 2 becomes quite large which results in a non-zero conversion amplitude from between IE‘7r+) and IEOWO). The overall expression for the phase shifts becomes quite long as the 3 particles also have spin 1/2. By combining Eq. (3.24) with eq. (3.21), one arrives at the result l‘f’lprevious-incoming) E |E_W+;l = 1,771, = 0> IT) M 2 . 4 . s— W: 1, = o) m 6 + 6 3 9 9 9 2 .. 2 2 + IE"7r+;l = 1,771,, = 1) H) (£62ml,l=3/2.J=3/2 + f __62151.1=1/2.J=:s/2 9 _ £e2i6l.1=3/2.J=1/2 __ .2—[262i6tlzl/2le/2 9 9 + l50n°;z = 1, m1 = 0) IT) (LB/5e2'6u=w=w — QTfieM-euz-hm + £62i6l.1=3/2.J=t/2 _ £621511=u21=u2 9 9 2 - 2 »- + IEOTrOJ = 1,mz = 1111) (-62'6’~’-3/2'J-3/2 — -62'°’-’=1/2»J=3/2 _ 362i51.1=3/2.J=1/2 + 382i5l.1=1/2.J=1/2) (3.25) 21 One can clearly see that if all phase-shifts are zero the original incoming state is obtained. Furthermore, if the phase shifts do not depend on I or J, the interaction is becomes elastic and the other three components are not mixed in. Resonance phase shifts are modeled using the Breit-Wigner formula. The relativistic expression is H tand, — m, (3.26) 3 (1 MR H 2 —.-—H, q} M 0 H0 = MRI“. Here, (10 is the resonant momentum and I‘ is the width of the resonance. The simpler non- relativistic Breit-Wigner formula, 1q2l+1 F0 “"6’ 2 26,31“ Mo — M ’ (3.27) can be used if M / M R ~ 1. 3.2.4 Constructing the wavefunction inside of the interaction range For relative distances smaller than the range of the strong interaction ~ 1 fm, the wave- function is no longer a physical quantity since the hadrons no longer interact like two point particles. This is especially true for resonances where there is significant microscopic rear- rangement of the quarks, e.g. 7r + 7r +—> p. To avoid artifacts that originate from integrating over this unphysical region with high weight, one has to replace the wavefunction by a smooth function which has to be consistent with the phase-shifts of the outgoing waves. 22 For T < 6 one can choose Mg, 1);? = mam, r)? + We. q) (3.28) where 11/,va is the Coulomb wave and W is constant for all r < e. W cannot be chosen arbitrarily since the wavefunction and the change in the density of states are related to each other [31]. Consistency with these parameters requires that the following equation is fulfilled: (IN (21 + 1) d6) 47rq2 / 3 , 2 . 2 __ = __ = d - ’l _ . ,9 A dq Z 7r dq (27r13 r (1141‘): r)| Idc(q,l‘)l ) 47rq2 47re3 = ’V (W1 3 l (6") oo 2w 1 + f d. r2 f (125/ (10059 (|«.;:.,(q.,r)|2 — Iwc(q,r)I2) c O —l (3.29) Afier inserting the expansions (3.10) and (3. 13) and applying the orthogonality relation for Legendre polynomials, 1 2 [dx f)l(1)l)m('r) — W611", (330) -1 one finally obtains dN 2q263 2 0° ~ - 2 2 _ = ' W _ ' I) a . _ F 1 dq 37. W+7r§lj<2l+1>fdr lamp» 7 on) I (qr m (3.31) 23 - 1 .. «.6, E F1+—(e2'°’—1)(F)—2’G,) (3.32) Combining (3.31) and the phase-shift relation in (3.29), one can solve for W by evaluating phase-shift derivatives and integrals of the type 1102261) a /d7'lwz((1r,m6z)12 (3.33) where (it, stands for either one of (,5, and F}. In order to avoid a numerical evaluation of these integrals, one can rewrite them using the recursion relations of Coulomb wave- functions: Since "([7 is an eigenstate of the radial part H r of the Hamiltonian, it obeys the one-dimensional radial Schrodinger equation: ’22 (92 [(1 ‘1' 1) ZIZQH ([2 111"], E ———‘ . ' .’ = E ’ : _'/ 3.34 L: 2p (91:2 + 2’1,er + r ([1 (If: 2” U1 ( ) The trick for solving the integral (3.33) for a given state 212) with outgoing momentum q is to take a second state 1/2,’ that has infinitesimally higher momentum q’ = q + dq. Then, (q'2 — (12)»? / dn/Jfll; = 2m; / dr [(H.w,")w, — w;*H,.2/),] (3.35) In the limit dq —> 0, all r.h.s. terms except for those involving the differential Operator cancel. Afier integrating by parts, a a“) ,,_ 0"”! am] . (3.36) 2011(E.Q.51) = ER [(37—071— E297 C 24 Evaluation at the upper limit r —* 00 is not necessary as it is canceled in the subtraction (3.31) of the two integral expressions for phase-shifted and non-phaseshifted waves. The lower limit can be evaluated numerically using recursion relations for Coulomb waves [1]. The explicit steps are quite lengthy and are shown in appendix A of this thesis. Now that the integral 1, is known, one can calculate W using I) and derivatives of . d6 . . . . the phase shifts, — The strong interaction correction radius 5 was set to 1 f m for all , z dq ' simulations. The dependence of the correlator on the choice of E is not significant as long as e is significantly smaller than any characteristic dimension of the source. 25 Chapter 4 The Classical Coulomb Correlator In the previous chapter, it was shown that Coulomb forces influence the correlation fiinction through the distortion of the two-particle quantum wave function. However, if the relative momentum of the two particles is large, the De Broglie wavelength becomes small, and a classical description should be reasonable. In this chapter, we present the calculations for Coulomb-induced correlations using a purely classical treatment and compare with the analogous quantum treatment. 4.1 Derivation of the classical correlation weight To determine the correlation weight, we can consider a source which emits particles of mass ,u and momentum qo at r0 according to a probability S (r0, qo). If the particles are then deflected from their straight-line trajectories by a Coulomb potential, V = k/r, their 26 final momentum distribution is dN (Pg 3 M . (4.1) : d37‘os(l‘0, (lo) d3q Solving this problem represents a solution to the correlation problem as the motion of two particles under a central force reduces to a one-particle problem. The correlation function is the ratio of the spectra with and without the presence of the Coulomb potential. _ dN/d3qO C(Q)-— dN/d3q- (4.2) If S is independent of qo, or depends only on a constant of the motion, e.g., the energy, the q dependence is S disappears and one can express the correlation in a simple expression. d3qo d3—q q?) dqo sin 90 (1190 (1% q dq sin 9 d6 dqb C(q) = /d3’rog(l‘o) = [(137090‘0) , (4.3) where g(ro) is a normalized distribution that represents the spatial shape in S. By compar- ing Eq. (4.3) with Eq. (2.14), one can see that the classical analog of the squared quantum relative wavefunction is |d3q0 /d3q|. To play the role of the squared relative wave function, the Jacobian |d3q0/d3q| must be expressed in terms of the asymptotic momentum q and the initial spatial separation r0. Here, we will follow the derivation of [3 5] and derive an equivalent, though more compact, expression. 27 The solution involves finding the three components of qo. This can be accomplished using the three constants of motion. The first constant is energy, (10 k (12 2 2W — — = — => = — — 4.4 2p + To 2p 1001 q 11'01 ( ) 2pk (13qu - 1 M q, qqu (4.5) Other constants of motion are the eccentricity vector [24] .quL+:=q_d:_Xi>+: (4.6) pk 1' pk r and angular momentum L. To simplify the algebra, we choose a coordinate system where the particle is initially located on the positive 2 axis, and the momentum is confined to the 122 plane. qo sin 60 0 (10 = 0 1‘0 1' O (4.7) (10 COS 60 To q sin 6 sin 6 q = 0 r /r 3° 0 (4.8) q cos 6 cos 6 Note that 6 is the angle between final momentum q and the initial relative position vector to. Equating the eccentricity vector, (4.6), for the time of emission and for very large times yields two conditions pk + 7‘0 q?) sin2 60 = pk cos 6 + roq qo sin 6 sin 60 , (4.9) m (13 cos 60 sin 60 = r0 q qo cos 6 sin 60 — pk sin 6 , (4.10) 28 where the angular momentum is qoro sin 60. Solving for sin 60 and cos 60, respectively, we obtain sin 60 = 2q_ {sin 6 i [sin26 — 27(1— cos 6)] 1/2}, (4.11) (10 k ' 6 c0860 = lCOSO — £781;— qo 7"qu 811160 cos 6 7 81116 = (_ _ , (4.12) Although these expressions are redundant, it is convenient to retain both of them. Here, the abbreviation 7 = 2pk/q2r0 has been introduced for convenience. Note that there are two possible solutions which correspond to the fact that for any given final momentum q there are two different trajectories which pass through r0. Inserting (4.11) in (4.12) gives cos6 7sin6(/1—7 1‘7 (1—7) [sin6:l:\/sin26—27(1-cos6)] cos6 (1+cos6):|:(/(1+cos6)2—27(1+cos6) cos 60 = _ . 4.13 1 — 7 2M1 - 7 ( ) Finally, we differentiate with respect to cos to obtain dcos6o _ I :1: 1 1+cos6—7 . (4.14) dcos6 — 2(/1— 7 2\/1—7 \/(1+cos6—7)‘2 _72 If the two charges are of the same sign, the potential is repulsive and 7 is positive. The argument of the square root in Eq. (4.14) can then become negative which restricts the limits of 6, (27 — 1) < cos 6 g 1. This constraint represents the physical fact that a particle of a given energy whose initial position is at 6 = 0 can not be emitted with asymptotic 29 momentum at 6 = 7r because of the deflection of the Coulomb force. We enforce this restriction with a step function, (1 cos 60 dcos6 i 1 1:1: 1+cos6—7 261—57 (/(1+ cos6 — 7):? — 72 x 9(1 + c036 — 27) (4.15) ‘6 99 The “+” solution is always positive and the solution is always negative, so absolute values can be taken easily, and by adding the two solutions, one obtains dcos60 _ 9(1+cos6—27) 1+cos6—7 (416) dcos6 tot (/1—7 \/(1+cos6—7)2—72 ° Although there is a singularity at cos 6 = —1 + 27, the expression remains integrable. In fact, to see that this result is physical, one can check that the integral (4.16) over cos 6 yields f dcos 60 = 2. Making the substitution u = 1 + cos6 — 7 and replace the step fiJnction with the integration limits, 1 [dcos6 —l (1 cos 60 (1 cos 6 _ 2 (2-7 2 72 _ 1 f _ 1 _ 7 2.6—7 1 : / d’UJ 1-7 0 ._ 7,2 dv (A? 2—7 _ / 1 u du tot \/1 ’7 \/u 7 0 II to (4.17) Combining this with equation (4.5) and the relation (1qu = dab which comes from an- 30 gular momentum conservation, we have as our final result for the J acobian: 11':qu q?) dqo (1 cos 60 dqbo d3q q2 dq d cos 6 (19$ 1 + cos 6 — 7 = \/(1+cos6—7)2 _72 @(1-1-0036—27) (4.18) 4.2 Numerical Evaluation With the inspiration from the previous section, we can now apply the result of Eq. (4.18) by considering a simple Gaussian example for the source. 2 2 1 r3," 'rside + 1long 1' = , ex - . — (4.19 9” (4n)3/2Roij p( 4R3... 433 ) The resulting integral is rather complicated and will be solved numerically using Monte Carlo techniques. Complications arise from the singularity in Eq. (4.18). This singularity can be eliminated by performing a transformation of variables as was done in Eq. (4.17) when checking the integrability of the Jacobian, w = \/(1+cos6—7)2—72, 0p (4.25) This final expression can now be evaluated numerically. There are no problems with limited floating point accuracy as the values over which is W averaged have similar order 33 of magnitude. A more detailed underflow test gave a relative contribution of numerical errors to the final result of < 10’8. Explicit data obtained from Monte Carlo simulation runs will be presented together with quantum mechanical results in the next chapter. 34 Chapter 5 Numerical Results and Discussion 5.] Choosing particle pairs for correlation analyses In principle, F SI correlation analyses can be performed with all observable particles. How- ever, there are some issues that have to be dealt with that are less significant in identical particle HBT interferometry. Truly chaotic emission of the source is important for the suc- cess of any correlation analysis. For final-state correlations, this becomes a crucial issue as effects from Coulomb and strong interactions can be as low as 1%. Competing correlations must therefore be analyzed carefully. There are mainly two mechanisms through which particles can acquire unwanted cor- relations before they enter the final state: 0 If two particles share a quark-antiquark pair, there is a good chance that they are sta- tistically correlated due to momentum conservation in pair creation. As an example, detecting a proton at high relative momentum enhances the probability that particles 35 containing 11 or J quarks are also detected in the same event at high relative momenta. Therefore, the use of particle pairs sharing a quark-antiquark pair for F SI correlation analyses is discouraged if the FSI peak is below l %. 0 Correlations from jets and collective flow can also cause problems in a F SI correla- tion analysis. Careful construction of the correlator from events that have the same reaction plane should allow to minimize the effects of collective flow. Correlations from these two effects should be of the order of 0.5 %. In general, experimental analyses should be restricted to ranges of Qinv E 2q for which competing correlations are known to be significantly lower than the expected signal. For Coulomb correlations, which fall off as 1 /Q12nv’ this typically limits Qinv to below 200 MeV/c. Further, it is important to remember that the wavefunction formalism is not covariant, so particle momenta in the center-of-mass frame should always be significantly lower than the corresponding particle masses. A truly covariant treatment of final-state interactions would be desirable. For small 62, relativistic corrections can be crudely accounted for by altering the definition of the reduced mass, E1E2 = _ 5.1 u E1+E2 () Fortunately, for the examples discussed in this study, the most important source information is hidden in the part of the correlator for which the invariant relative momentum Qinv is still non-relativistic. 36 5.2 pK+ correlations Proton-kaon particle pairs are promising candidates for correlation analyses due to the fact that fewer particles originate from decays of long-lived resonances than in pion-based anal- yses. In the context of HBT interferometry, the advantage of using kaons instead of pions has already been discussed some time ago [26]. Further, pK+ correlations do not have correlations from charge conservation in pair creation as uud and us have no common quark-antiquark pair. The statistics for constructing the pK+ correlation function should be available with current RHIC data sets. Correlations of pK + particle pairs are dominated by Coulomb effects. Strong effects were included using phase shifis for l = O and l = 1 from experimental tables [18], but as it can be seen from Figure 5.1 they affect the correlator only moderately. In consistency with the findings of [14], the classical description of Coulomb effectsis remarkably close to the quantum mechanical solution for large relative momenta Qinv- For measurements of the source shape, the dependence of the correlator on the angle between q and the outwards direction has to be analyzed. As already shown in eq. (4.5), the angle-averaged dependence of the correlator goes as 2d 2k k q°q°= 1——”—z1—E— (5.2) (12d!) 927‘0 (127‘ Plotting anv(C (q) — 1) instead of C (q) should largely compensate this Qim, dependence and emphasize the physical structure of the angular correlation plot. All of the following angular correlation plots will therefore employ this scaling. 37 IIrIITfIIIIIIIIfiII LO - 023 - A> - -5(16 9, T 0 0.4 __g — Classical (no strong) .— ' Cl Quantum (no strong) -‘ 0'2 _ 0 Quantum '— 0.0 111111111111L11L14 O 50 100 150 200 Qinv (MeV/c) Figure 5.1: pK+ correlations for a Gaussian source; Rlong = Rside = 4 fin, Rout = 8 fm The pK+ angular correlation plot (fig. 5.2) shows again the good agreement with the classical approximation for large relative momenta. The minima at cos 6 = :tl are caused by the fact that the momenta of particles moving along the relative position axis are de- flected away from this direction by repulsive Coulomb interaction. The depth of the dips at cos6 = i1 is proportional to Ram/Rifle for large values of Qinv- Careful investigation of the classical correlator shows that for large q, where the correlation is confined to a small region where q is parallel to re, the probability for the particles being separated by a given distance scales under a 90 degree rotation of the 3 3 source as Rout/ Rside whilst the classical phase—space enhancement scales as Esme/Rout. In fact, the simulation for Rout/Rside = 3 gives an angular correlator that is nine times deeper for cos 6 = 1 than for cos 6 = 0. This favourable scaling behaviour of the angular dependence of the classical correlator shows that Coulomb induced correlations have strong 38 9\ Q -200 > (D 3 —400 '3‘ o A 8 —200 O :V—L 01400 I l l l I I I I Qinv=3O MeV/c _ 5.3 p7r+ correlations "U 1:1 0 a 1:1 1:1 :1 E O O O 0 O O O O L — Class. (no strong) ‘ __c1 Quant. (no strong) fl 0 Quantum l l l l L 1 l l r I l I fl l I _ Q. =1SO MeV/c .. mv 0 o o o o o o o "— U‘ l l l l l l L l 0.5 1.0 cos(6) Figure 5.2: Angular pK+ correlation; RIO" = Rside = 4 frn, Rom = 8 fm g diagnostic power for resolving asymmetries of the source. For large Qinv’ the quantum results approach the classical result and thus have similar power to discern the shape of So far, the strong interaction has contributed only minor corrections to a Coulomb- dominated correlator. We will now see that in the vicinity of a resonance, the strong inter- action becomes the dominant factor and Coulomb effects become almost negligible. The 39 impact of the A++ resonance on correlations of 1271+ particle pairs will be investigated here as an example of correlation analyses in the vicinity of a resonance. Phase shifts for the l = 0 and l = 1 channels were again taken from experimental tables [18]. Near the resonant momentum, Qim, = 450 MeV/c, the angular dependence of the cor- relator becomes quite large as one can see from Fig. 5.3. For resonant scattering, the ‘ ‘1 q (>8 : 1 I r 1 l I Z 8 a > - 8 8 8 _ g : 8 8 D : v 1000 [- 0g -—_ 2:1 I 83 i l I. _‘ ,~% 0 _ 3 4 g, _1000 g_ . Qim=450 MeV/c Sj :4 I 040 Cmfiowb E E I 1 O _2000 1 1 1 1 1 1 1 1 1- 0.0 0.5 1.0 cos(6) Figure 5.3: Angular dependence of the 1271*” correlator near the A++ resonance; Rlong = Rside = 4 fm, Rout = 8 frn angular dependence of the correlator arises due to shadowing in the forward direction. One can verify that shadowing is indeed the reason for the dip at cos 6 = i1 by running the correlation simulation for a greatly exaggerated scenario where the source is displaced by Zoffset = 40 frn in the direction of Rom. If shadowing is the real cause for the dip in the angular dependence graph, one can expect that for the displaced source this would result in a very narrow dip at cos 6 = -—1 where the particle first travels in the direction of its partner 40 and is then scattered. All other areas of the correlator should at most show small enhance- ments over unity that are proportional to da/dQ as all scattered particles are deflected in some other direction and give a contribution to the correlator in this direction. Coulomb interaction has to be neglected since its long-range repulsive nature would modify the result significantly and therefore complicate the investigation of strong interaction effects. é é _ Qinv =450 MeV/c, no Coulomb Qinf 1C _ _ § 13500 ['1 r“ r :130000805588855pm- - v—"l _ O O 8 (:3 C] A O D <13 12500— 0 0 — > O B .E F o 9 12000— D — U . ,_.— + 0 E, *-' :1 a. 7: at Q. = 295 MeV/c N >115m_ _— + "IV 0'1 0.33 0 c. it , Coulomb corrected 4 l L l 1 l l L l l 10000 0.2 0.4 0.6 0.8 1 cos 9 H Figure 5.5: :‘7r+ correlations at resonant momentum; Rlong = Rside = 4 fm, Rout = 8 fm 42 Chapter 6 Conclusion and Outlook Two-particle correlation functions provide tremendous information on the physical nature of the freeze-out state of a heavy ion collision. The three most important parameters that can be extracted from the correlations are the lifetime AT of the source as extracted from the ratio Rout/ Rsidea the average emission time < T > that is obtained from Rlong, and the overall volume of the emitting region. If a quark-gluon plasma is formed in the collision process and if the phase transition is of first order, it would have a large amount of latent heat. This would result in low pressure and a long lifetime [44, 54, 56] of the order ~ 10 — 20 fm/c. Particle emission should be continuous over the entire lifetime of the plasma. The first results of HBT measurements at the Relativistic Heavy Ion Collider (RHIC) [6] were a real surprise for the heavy ion physics community as they showed nearly equal outwards and sidewards source sizes. This places an upper limit on the lifetime of the source of AT < 10 fm/c which is significantly less than expected. Moreover, the average 43 emission times (7') were also unexpectedly small. The fact that none of the current dy- namical models is able to match both the observed particle momentum spectra and these small HBT radii is referred to as the “HBT puzzle” [47, 56, 37]. There are two possible interpretations of these facts: either HBT measurements are unreliable for some unknown reason, or our picture of a relativistic heavy ion collision is incorrect. The first interpretation can be probed by measuring source sizes through alternative methods. In this thesis, it was shown that non-identical particle correlations from final- state Coulomb and strong interactions carry important information on the source that allow the extraction of its size and lifetime and thus provide an alternative method. An article that has been submitted for publication [51] will summarize these results. The HBT puzzle is challenging our models of the quark-gluon plasma. The fact that our models overpredict HBT radii could be caused by the mere non-existence of a QGP. Hydrodynamical models are based on minimum-increase-in-entropy propagation, yet they overpredict the measured HBT radii and the measured per-particle entropy. This would indicate that the true initial state had much lower entropy and could not have been a quark- gluon plasma. On the other hand, other observables like the charge balance functions that will be discussed in the following part of this thesis seem to indicate that a QGP might indeed have been formed in RHIC Au+Au collisions. Intensive work will be required in order to finally resolve all these issues and to give a final answer on the existence and nature of the quark-gluon plasma. Part II Modeling Balance Functions in Heavy Ion Collisions 45 Chapter 7 The Balance Function Formalism It was already mentioned in the first part of this thesis that the constraints of charge con- servation in pair creation cause significant statistical correlations in momentum space be- tween particles and their balancing antiparticles. Whilst these correlations are harmful in the context of two-particle F SI source measurements which rely on chaotic emission, they may carry important information on the timescale of pair creation. Charge balance func- tions were introduced as a new observable [13] that allows to quantify such correlations. The strong sensitivity of these functions with respect to the timing of hadronization makes them good candidates for a QGP signal. The current chapter will give a brief overview of the balance function approach whilst the following chapter will address specific issues that are relevant for Monte-Carlo simula- tions of balance functions. 46 7.1 Introduction The canonical picture of a relativistic heavy ion collision that forms a QGP assumes that the plasma phase would be quite long-lived as its internal pressure would be low due to the high amount of latent heat stored in the plasma. As the plasma expands and cools, it returns to ordinary hadronic matter through a phase transition that might be of first order [33, 34]. This phase transition is accompanied by a massive increase in the number of quark-antiquark pairs as the entropy stored in the plasma is converted to mesons. Since gluons carry no charge, and hadrons contain at least two quarks, most of the electric charge, strangeness and baryon number is therefore produced at a late stage of the collision if a QGP intermediate state exists. The timing of hadronization is the most obvious difference between a collision that involves a QGP and a collision that involves just purely hadronic intermediate states and that is characterized by early hadronization ~ 1 f m / c. If charges are created early, they can separate further in rapidity due to the collective expansion of the collision process. A simple toy model can show how drastic this effect should be: In a Bjorken expansion model, the rapidity gradient along the beam axis is approximately given by 1 / T. If quarks separate by 0.5 f m at a time of 0.5 f m/ c, their final difference in rapidity should therefore be of the order of magnitude of 1. On the other hand, if late hadronization due to a QGP intermediate state occurs at 7' ~ 5 f m, the difference in rapidity of a quark-antiquark pair initially created at relative distance of 0.5 f m would be just ~ 0.1. The purpose of balance functions is to identify a meson’s balancing charge and to mea- 47 sure its difference in rapidity on a statistical basis. One defines the balance function as B(PZIP1)E%{1V+—(P11P2)_N++(P11P2) N-+(P11P2)_jV-—(P17P2)O} (7.1) N+(P1) N-(Pl) Here, N+_(P1, P2) is the number of events where a positive particle of a given species was observed in the bin P1 and its balancing particle was observed in the bin P2. There are many possible binning criteria, but we will focus on relative rapidity in this thesis. Also note that the counting for Ni does not necessarily have to refer to electric charge. The balanced “charge” could for just as well be strangeness, baryon number, or isospin [3. One can construct balance functions using all positive and all negative, all strange and all antistrange particles, or with particle and antiparticles of a given species, e.g. n+7r‘. The normalization of the balance function would be unity if the colliding nuclei had no surplus charges and the detector had an acceptance of unity, i.e. all particles were detected. For a lower detector acceptance, the balance fiinction will sum to the fraction of particles carrying the globally conserved charge that is actually detected. Furthermore, if the balance function is constructed using a subset of particles, e. g. K +K ‘ , the normalization would be reduced since the balancing charge might be carried by another particle. 7.2 Balance functions as a QGP signal There is an especially interesting class of balance functions that uses the specific binning that P1 should be the detection of any particle in the experimental acceptance, and P2 should be the detection of a particle with relative rapidity Ay. The properties of this specific 48 charge balance function have been investigated in several papers [13, 46, 32, 49] and have also been measured at RHIC by STAR [4]. The expectation is that if a QGP is formed in a relativistic heavy ion collision, the peak of the balance function at Ay = 0 should be much narrower than the peak of a bal- ance function of a pp collision. Moreover, the peak should also become narrower with increasing centrality of the Au+Au collision the higher temperatures would result in an enhancement to the lifetime of the QGP and an additional delay in emission time. For a purely hadronic collision scenario, emission times should not change and the balance func- tion should not become narrower with increasing collision centrality. Simulations using hadronic cascade models have confirmed this expectation. Figure 7.1 shows the balance function as obtained from the hadronic cascade models HIJING and GROMIT which sim- ulate immediate hadronization according to pp phenomenology. I I I I n I I I I I I I 04 33 1:1 0 HIJING+GROMIT1‘ [08080 :1 HIJING . ’9. . 081:1: . g 02 L oog — _ 00 - _ Una . Figure 7.]: Balance functions obtained from hadronic cascade models 49 A first measurement of balance functions at RHIC has been performed by the STAR col- laboration [4]. It was found that the balance function becomes narrower by approximately 20% for central collisions (see fig. 7.2) whilst for peripheral collisions it approaches the width that is predicted by hadronic models. 0.7 0.68 0.66 9, 0.64 $ 0.62 0'6 0 Data 0 58 E HIJING-GEANT 0 0.2 0.4 0.6 0.8 1 b/bmax Figure 7.2: Narrowing of the pion balance function with increasing centrality of the colli- sion as measured at STAR The theoretical modeling of balance functions for any kind of emission scenario is quite complex since one has to simulate particle multiplicities and momenta, and disturbing effects like the unbalanced charges of the two nuclei [49], the limited detector acceptance and effects from HBT correlations [32] have to be accounted for. A problem in balance function simulations is that the observed particle multiplicities have to be modeled in a way that the constraint of charge conservation is observed. Models that just create particle-antiparticle pairs neglect the possibility that the balancing charge of a particle may actually be spread over several other particles due to charge mixing and 50 - I - Data Central 3) - -B- Data Peripheral ‘ b). 0-3 -I -HIJING-GEANT 1 E} 0.2 '1. , v 1 m 0.15 ‘ ; j I : 0.1 i‘ x" —f 0.05 In. O . . . L 1 . . i . i . . . . 1 . . 1“ 0 0.5 A1 1.5 2 y Figure 7.3: Balance functions for charged pions as measured at STAR in comparison to theoretical data from hadronic cascade models chemical equilibration at an early stage of the collision. Therefore, we have developed a new thermal model for generating particle multiplicities which accounts for charge mixing whilst charge conservation is still strictly enforced. We will discuss the theory, implemen- tation details and present preliminary results of this model in the following chapters. 51 Chapter 8 Modeling Particle Multiplicities Particle multiplicities generated according to a thermal model using a canonical ensemble have been shown to reflect observed multiplicities quite well [16]. Local charge conserva- tion can be modeled by providing for each particle a balancing antiparticle, so for instance every 7r+ would be accompanied by production of a 71‘. However, this simple model does not account for the effects of charge mixing in the process of chemical equilibration. As an example, the a K + could be balanced by a p and a A. To include such effects, a new algo- rithm was developed in order to simulate charge-conserving particle creation according to a canonical ensemble. 8.1 Recursion relations for partition functions Let us consider a system which consists of A particles that have a total charge (baryon number, strangeness, isospin 13 etc) of Q. Assuming that there are P different species of particles with charge (7,. and a single particle partition fiinction wk, the total partition 52 function for A particles is _, (want: Z = . 8.1 AlQ) Z; a [£11 "k! ( ) Z nqu = Q, k Z 71ka = A k The number a), must be a positive quantity for all particles. For calculating the total par- tition function Z = 20 Zn, the assignment of 0.), becomes arbitrary. We will choose to assign a.)c = 1 to all particle types that are sufficiently long-lived to reach the detector, and (1),. = 2 to extremely short-lived particles. With this choice, 2k a), roughly reflects the number of particles that are created, though some unstable species might have some 3-body decay channels. Expression (8.1) is almost useless for practical purposes since solving the combinatoric problem of finding all possible 71,, has exponential complexity. However, by using recursion relations and caching all results in memory, the computational effort required for solving the problem can be drastically reduced [19, 23, 20, 52]. The recursion relations are easily obtained by inserting a factor of 1 and rewriting the expression as F "1: 21(0) = Z (2 my.) H (0:2! (82) 11kI—l F m, wklap 0.1.; (U. 2 Z A Z (7111— 1)! H 71:! (8.3) k, 27116;. = Q. If :1. k A ,,I anak = A . ¢ 1: F 71k _ 0.21.1011 wk _ Z A Z - H ”k! (8.4) k Emil: = Q — (17:, k=1 I: 2 mm, = A - (1;, k = Z wkjk’ 2.1—add). '— (17.) (8'5) kl This simple recursion relation can be used to build a partition function table in memory for a relatively high number of particles (A ~ 50). However, memory and CPU time consump- tion can be tremendous if poorly optimized code is used. Since algorithmic optimization is crucial for generating partition functions for larger A on workstation-class hardware, the specific optimization techniques employed in the final C++ implementation are discussed more in detail in the next section. 8.2 Optimization issues One of the key issues in creating an efficient implementation of the recursion relations is the reduction of the number of partition functions actually calculated. Quite often, there is no combination of charges so that 21. 111.17), = Q — q}: and 2k nkak = A — up. In these cases, the contribution of Z A..,, k, (Q — (17,") in 8.5 is zero. Since the number of zero values by far exceeds the number of non—zero values, an efficient implementation must avoid calculating and storing zero-value contributions. 54 Therefore, the first step toward an optimized implementation involves changing from the original recursive algorithm to an algorithm that starts at the bottom (A = 0) and works its way up to the top (A = Amax) by only calculating the non-zero partition functions. Its basic working principle can be described by the following pseudo-code: Allocate and assign Z[A=O](Q=O) = 1; Run for A from 0 to A;MAX For (every non-zero value ZIA](Q) in memory) For (every particle species k) { Allocate Z[A+a[k]](Q + q[k]) in memory if necessary; Make the respective contribution to Z[A+a[k]](Q + q[k]); In comparison to an algorithm that calculates all values of Z for the full charge range that could possibly be combined from A particles, the requirements in CPU time are reduced by a factor ~ 50. I Results should be stored in memory in a way that allows iterating the non-zero values of Z g(Q) efficiently whilst providing fast access to the individual values. The C++ imple- mentation uses an array of STL template classes of type map< charge , double> which map a user-defined charge class to a double that holds the non-zero partition fiinction value. Iterating through all non-zero values using STL container iterators is eflicient and easy. As no zero values are stored in memory, total memory consumption is reduced by several orders of magnitude compared to an implementation that relies on arrays. Data access for a single item is of complexity 0(log n) which is acceptable even for intensive Monte Carlo simulations. 55 The fact that particle and antiparticle partition functions are always identical allows further optimization by skipping the calculation of Z (—Q) and simply returning Z (Q). However, contributions to partition fiinctions that cross the axis of charge symmetry can lead to double-counting and finally wrong results if the optimization is not done carefully. The implementation of charge-anticharge symmetry brings a performance boost of ~ 70% and cuts memory consumption by another 40%. The C++ implementation that was developed in the context of this thesis takes care of three charge conservation constraints: Isospin 13, baryon number and strangeness. The par- ticle species that were included for the simulation are listed in appendix B. The calculation of the partition function up to A = 55 using the optimized C++ implementation took about 20 minutes and consumed 950 MB of memory on an Intel-based 2 GHz system. It can be clearly seen that without the optimization techniques described above which all together reduced memory consumption by a factor ~ 10, the calculation of the partition function for such high values of A would not be possible on a workstation-class computer. 8.3 Monte-Carlo generation of particles Particle multiplicities are obtained from the constructed partition functions using the fol- lowing algorithm: First, the total particle number A is chosen according to Z A(Q = 0). This is physically correct as the total partition function of a charge-neutral system with an arbitrary number of particles is simply Z( (8.6) 101 ll 8 II a? <01 ll 8 where Z A gives the weight of the A-particle state. When the particle number A has been determined, one can Monte-Carlo the individual particles: aka/’1; _. A ZA‘ak(Q _ (173‘) 1. A particle of species It is chosen proportional to the weight 2. A is replaced by A — oh, and Q is replaced by Q — (17,, then step 1 is repeated until A = 0 is reached. It will now be shown formally that this algorithm chooses every particle configuration with a probability that is proportional to its contribution to the partition function Z A(Q). Again, we denote the charge vectors for the different particle species as 171.. Note that these charge vectors are not always linearly independent. A set of particles will be referred to as {1.7, . . . ,kN} where the i-th particle is of type k,, 1 S 2' S N. Note that this def- inition is sensitive to the ordering of the particles. For a given value of total charge Q and associated number A we will call the set {k1, k2, . . . , 1.70} a path to (Q, A) and write N Q = Z 7,, and A = Zak, (8.7) Note that for fixed (Q, A) there is a large number of possible paths, and not all paths necessarily have the same number of particles N. Let us now define the path contribution Z for a given path Di] ({kl, . . . , kN} ; Q) as Z (8.9) all paths _. Dflfki, - - -,kN};Q) It is possible to factorize out the contribution of the last particle of this path: N 21(0) = Z Z Hgk‘a": (8.10) I N all paths {:1 Z aki 1251'({k1,.-- ,kN};Q) ’ N = z: 2 H5” A (W k’ N 1:1 i=1 “*1 all paths 251:3,(111,....k~_.};0 — 41') a w , ~ ~ _. = z .A. 2..-...10 — q.) (8121 kl We see that 2,; obeys the same recursion relations as Z A. Since Z A is well-defined by the recursion relations and the initial value Z0(0) = 1, the new object 2,; must be equal to Z A. Thus, the partition function ZA(Q) is really the sum over the path contributions 1:3(Df;I ) of all paths to (Q, A). Now we calculate the probability that the algorithm will select a specific path Dfi({k1, . . . , kN} ,Q) that has a path contribution of Z(D). Starting at (Q = O, A), the probability to select the particle species 117.; for the last ( N -th) particle is _ who“ ZA-akN (Q — 6N) ’P(kN)- A 2.40:0) . (8.13) 58 Similarly, the probability to obtain the particle species kN for the last and kN_1 for the next particle is .0 — 6.1V (N—l)(Q — (RN — (jkjv_l) A A - (1,, 2.40 = 0) 21-.., (0‘ — (it) (8.14) ) ZA—akN —ak One can see that the partition function Z A-..” (Q — (71w) cancels. Iteratively one obtains for the probability to choose a given path Dg’({1~,, . . . , 1910} ,Q’) the value Z('D§({k1,- - - ,k1\’}ié)) a , (8.15) ZA(Q = 0) P(D§({k1,...,kN}7Q)) : which is proportional to the physical weight that this path contributes to the partition func- tion Z,.((Q = 0). Thus, the algorithm chooses paths according to the definition of Z in Eq. (8.8) which is consistent with the contribution of the configuration to the partition function. 59 Chapter 9 Results and Conclusion 9.1 Constructing the balance function In this section, we will briefly sketch the other steps that are necessary for constructing balance functions from the thermal model that was discussed in the previous chapter. After generating particle multiplicities according to the described Monte Carlo procedure at a chemical freeze-out temperature of T j, = 175 MeV, the particle momenta in the source frame are assigned according to a thermal distribution at lower kinetic temperature T)c = 120 MeV which is the hadronic decoupling temperature measured at RHIC. The thermally generated particle multiplicities contain a significant amount of very unstable resonances that decay within a short period of time. This is accounted for by an algorithm that performs decays according to a list of known particle decay channels. Decay products are assigned momenta such that total momentum is conserved and the momentum correlations resulting from decays are reflected accurately. 60 Finally, a center-of-mass frame for the point-like source is chosen according to a blast- wave description at Tblast = 120 MeV, and all particles from a given source are boosted into that frame. Our specific choice of a blast-wave model is such that sources are assumed to be uniformly distributed in the transverse direction up to a maximum radius R, dN _ constant forr < R (9 l) r dr _ 0 for r > R ' The transverse rapidities are assumed to increase linearly with r, T ll/trans = E11 trans,max (9-2) Longitudinal rapidities are chosen randomly between y = +2 and y = —2. After boosting the source to the frame according to the blast-wave model, particle pair interactions are neglected for simplicity. Interdomain correlations [50], i.e. correlations arising from long-range Coulomb interaction between particles from different sources have also been omitted in these simulations. The acceptance of the STAR detector has been modeled for the (pseudo)rapidity range of —1 < 17 < 1. The thermally generated balance firnction of identified pions resulting from this simu- lation procedure shows good agreement to the actual STAR data as one can see in fig. 9.1. 9.2 Conclusion Balance functions were introduced in order to obtain information on the timing of particle generation in relativistic heavy ion collisions. Although the formalism is still very young, 61 IITTIIIrIIIITIIT 1:] STAR, peripheral ' A STAR, mid-central — 0 STAR, central . — BLAST WAVE Figure 9.1: Thermally modeled balance function of identified pions in comparison to STAR measurements balance functions seem to be already quite well-understood given the good agreement of theoretical predictions and actual STAR balance data. However, there are still aspects that have to be investigated deeper before the balance function can actually be called a a “QGP signal”. Aspects that may affect the broadening of the balance functions like interdomain interactions have to be examined carefully in order to fully understand the width balance function in a nearly model-independent fashion. Once all these issues have been resolved, it should be possible to clearly state whether or not the observed narrowing of the balance function with increasing collision centrality can be seen as a signal for the formation of a QGP. At the moment, balance functions can be seen as promising observables for the observation of a quark-gluon plasma. 62 APPENDICES 63 Appendix A Recursion Relations for Coulomb Waves A.1 Transformation to natural coordinates Let D := {(r, q) E IRQI q ¢ 0} and let w : D —+ R2 be the coordinate transformation defined by 2 Let u : R2 —> C be some complex function, and 11 : M —+ (C the function induced on M by the coordinate transformation, 11 := u o 11}. From the chain rule, (117 = du odw (811 811 0 meQ/qz 5W7) 64 _ (12 _?fi%+,.?_u — (180’ (12671 01) area = (.—,., .—.-> A.2 Evaluation of Eq. (3.36) Let now 11,, : R —> C be a solution to the “Coulomb wave” differential equation 217 L(L +1) (“)2 p p2 w =0 (Al) The solutions 11,, (Coulomb wavefunctions) then fulfill the following usefirl recurrence re- lations: 2 d L L d“; = \/ L2 + 712111,-1 — (? + 77) 11,, (A2) 2 (L+1)4—‘ = (“21) +77)UL—\/(L+1)2+772UL+1 (A-3) We want to rewrite Eq. (3.36), . (‘9 -, 3 - -, 8 6 ~ (1 (>1) in coordinates (n, p) such that it contains only partial derivatives of 111, with respect to 77 and terms in 11,, or 111,1. This can be achieved by combining the results of the first section with the recurrence relations (A2) and (A.3): 65 In the new coordinates, equation (A.4) becomes -_ 2333314 :1. Bi- 211. J _ 17 6p r1162 pap 817 —u* mezg~ 77 patrL _ 912 L '1] (9p 77162 (9p (917 For more convenient handling, we will treat the two terms in this expression separately: . (9112 611.1, BuL QUL 82111, 02111, ) = — — — — I — — A5 J 8p (p (9p 371 ) :1" ( 61) + p 3P2 778/1811 J ( ) A if A.2.1 Case 1: L # 0 Term A: (911‘ (9111, 811 L A = __l; __ _ l__ 6p (p 0p an) ‘2 = 1 (ML? + 712 111-1 - (147+ 77) HZ) [136/112 + 712 ILL-1 L (L2 ) ) 6111,] — —+n m. —n—.— p 071 2 = i [(L2 + 712) ”Liar—1 - v L2 + 712 (L? + 71) (“11114—1 + "Z-rub) L2 L2 2 * 7] /—— I. 811;, L2 .6111, J P Term B: (911. L (921171 82 11L (911. L 3211.] 0211,, B = u: ++ _ — :11: ——+ ,:—-,, L ( 010 p 0122 "011017 L 'p ' n . . . ()1 This time we use recurrence relation (A3) to express a 2 NHL—1 _ L: . LL Li +L L 8p (p 77 d ZUL 7 811. _ 77 L2 811 _Z (W 81771 L2 2 L") (— ) 01 —UL)] * 1 2 2 L2 L L + 77 '11L_1— — + 77 UL p 1 L2 [ L2 + 7721: ((7 + 77) “IL—1 — V L2 + 712111.) L 1 . L2 L2 p . .7) , (.r—iz . - (7 +1) u.) + ,4] _ L2 2 — _ — — -— - L( +7] 877 + L2+772UL1 (p+77) 37? 1%)} 1 . , L2 . L («L2 + 772 “LUL-i — (y + 77) ULUL) + L2 2 L3 . 6 117 + n) + 714.24 - (1.2 + "21“”:1 77 . 011L_1 77 L2 871,, —— i/L2 + 772 11" , + 11" 11.1,-1 — (— + 77) u' ,— — 11" 71;, L( L ()7) \/L2 +772 L P L ()7? L S: 5 PM Terms A and B have similar structure which makes it easier to see that some terms are identical and cancel in the subtraction. The final expression for j = A — B for L 31$ 0 is 1 L2 — 417+) «— > i +L2 [(142 + 712) (11'Z_1UL_1 + 11211,) 67 . . L2 . —\/L‘2 +13 (7 + 77) (111211114 + 21L_1u1,)] _2 ‘ /L2 + 72 ”U“ % _ u" 8111)-} _ —L—U‘ 'U +171." L I L4 077 L 07] \/L-2+772 1L [Fl L L L - —11211L (A.6) p A.2.2 Case 2: L = 0 In the case L = 0, the steps are similar, but the other set of recursion relations has to be used. One then obtains j = (V1 + 772115111 — (7:- + 77) :10110) +p (1 +112>