{'fraction_non_alphanumeric': 0.035638703123677305, 'fraction_numerical': 0.016574959790061797, 'mean_word_length': 4.44136342699217, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Lectures on the Theory of High Energy A+A at RHIC

In these lectures aspects of current theory used to interpret the preliminary data on ultra-relativistic nuclear collisions at RHIC energies in terms of the physical properties of QCD matter at extreme densities are introduced. Topics include: What are the physics questions at SPS and RHIC?, Geometrical vs Dynamical features of A+A, The interplay of computable Hard pQCD vs phenomenological Soft dynamics, Baryon number transport and Junctions, How can we compute and get experimental control over the initial conditions?, How to reconcile apparent hydrodynamic behavior with partonic/hadronic transport theory. I use the Preliminary RHIC data available up to June 1,2001 to illustrate these topics. Since the main new observable at RHIC relative to SPS is jet quenching, I elaborate on this ``tomographic'' probe of ultra-dense matter. The possible discovery of jet quenching at RHIC by STAR and PHENIX is highlighted.


Introduction
Finally, after 20 years of preparation [1], a new chapter in nuclear/particle physics commenced on June 12, 2000 with the measurement of the first Au + Au collisions at √ s = 56 AGeV (GeV per nucleon pair) in the Relativistic Heavy Ion Collider (RHIC) at the Brookhaven National Lab (BNL). Soon thereafter collisions at √ s = 130 AGeV were also measured. The first results were reported at Quark Matter 2001[2] from the four major experiments, STAR [3], PHENIX [4], PHOBOS [5], and BRAHMS [6]. A small army of ∼ 1000 experimentalists measured the flavor, rapidity and transverse momentum distributions of the approximately 4000 charged particles produced in each central (head on) collision at 130 AGeV. In the summer of 2001, it is anticipated that RHIC will reach its design energy, and p + p and Au + Au collisions at √ s = 200 AGeV, will come under experimental scrutiny. These lectures provide a very condensed introduction to current theoretical work aimed to provide a consistent interpretation of observables measured in such reactions in terms of the properties of dense QCD matter. The color slides of the original lectures can be found on my WWW site [7]. These lectures note are designed to supplement those slides and update them with the preliminary RHIC data available as of June 1, 2001.
The theoretical work on the new physics that may exist in QCD matter at extreme densities began in the mid 1970's with the realization that the asymptotic freedom property of QCD implies the existence of a new phase of strongly interacting matter called the Quark-Gluon Plasma (QGP) [8]- [13]. Unlike familiar nuclear or more generally hadronic matter consisting of composite "elementary" particles (π, K, ρ, p, ∆, Λ, · · ·) in which quarks and gluons are permanently confined, the QGP phase at very high temperature and/or baryon chemical potential (T, µ B ≫ Λ QCD ∼ 200 MeV = 1/fm) is one where the interactions between quarks and gluons become relatively weak and short range The color electric (Debye) screening mass µ D (T, µ B ) increases linearly with T or µ B modulated by a slowly varying factor of the running coupling, g(T, µ B ) (see lectures of Rebhan). The thermodynamic properties of this deconfined and chirally symmetric (∼massless) phase of matter are thus expected in perturbation theory to reduce approximately to an ideal Stefan-Boltzmann gas of quarks and gluons. For the Standard Model with 3 colors and N f flavors of "light" quarks relative to T, µ B (SU c (3) ⊗ SU f (N f )), the Stefan Boltzmann constant for µ B = 0 is taking the helicity, color, flavor, and antiquark degrees of freedom into account. In reality the severe infrared singularities of perturbative QCD (pQCD) lead to large non-perturbative corrections to the ideal gas equation of state for all temperatures and chemical potentials accessible experimentally even beyond the future Large Hadronic Collider. Only numerical lattice QCD (lQCD) methods [14] (see lectures of F. Karsch) can provide reliable predictions for the thermodynamic properties of the QGP phase of matter. Effective models and resumed many-body techniques (see Alfors, T. Rebhan, J.P. Blaizot, and E. Shuryak) are however needed to interpret the lQCD "data" and provide physical insight especially at finite chemical potential. However, it is sobering to recall that pQCD thermodynamic expansion of the pressure in powers of g shows no sign of convergence [13] even before the Linde infrared catastrophe at O(g 6 ), and non-perturbative corrections to the pQCD Debye mass, µ D , remain about a factor of ∼ 3 up to T ∼ 200 GeV [16]. The full theoretical understanding of the structure of the non-abelian plasma phase of QCD therefore remains a fundamental open problem in physics because it involves strongly correlated, nonperturbative and possibly turbulent dynamical features [15]. One of the essential and intriguing aspects of the QCD many-body problem is that the physical vacuum is an extraordinarily complex coherent many body medium. The gluon and quark condensates lower the energy density of the physical vacuum by and amount B ∼ Λ 4 QCD ∼ 200 MeV/fm 3 . Drilling a perturbative vacuum bubble of volume V in this condensate costs an enormous energy BV . The QGP, if formed in V , must counteract the physical vacuum pressure, B to prevent it from collapsing. This is only possible when the temperature exceeds T c ≈ (3B/K SB ) 1/4 ∼ 150 MeV. The driving force behind the experimental effort at CERN and BNL over the past 20 years has been to try to create the extreme conditions necessary to produce and diagnose this new form of matter in the laboratory. Over the past 15  In these lectures I will focus on the most recent developments in that search that has just begun at RHIC. I must emphasize from the onset that most of the data shown here are of PRELIMINARY nature and could change as better control over the systematic errors is achieved in the next few year. Nevertheless, the new data are so exciting that it is worth trying a first pass to give an overview and possible interpretation.

Geometry and Dynamics in A+A
The main obstacle in interpreting data on collisions of finite nuclei (at any energy) is that the matter created undergoes quantum (perhaps semi-classical) many-body dynamics that may be approximated by thermodynamics only over a limited (low p T ) kinematic range. Experimentalists do not have the luxury of lattice or perturbative QCD theorists of tapping into the infinite gedanken volume or reservoir with a fixed temperature and pressure. Nuclear collisions produce dense matter in a highly dynamical environment, and the matter produced expands anisotropically near the speed of light. It is far from clear whether local thermal and chemical equilibrium concepts apply, and even so, over what domain of the 8 dimensional (x µ , p µ ) phase space can they be used. Before a collision, the partons of the two colliding nuclei are locked into a coherent field configuration. The dense virtual cloud of gluons and quarks may be described in the colinearly factorized QCD approximation by A times the known structure functions, f a/p (x, Q), of nucleons when the resolution scale is high enough Q > 1 − 2 GeV. However, many body initial state interactions could lead to strong modifications of this naive parton picture (see lectures of L. McLerran). The nuclear QCD fields continue to interact after the nuclear valence quark pancakes pass through each other. The interaction spans a space-time hyperbola over a proper time √ t 2 − z 2 ∼ 30 fm/c = 10 −22 sec. Then a "miracle" happens! The field quanta hadronize in way that is unfortunately not well understood. The dense final hadronic debris can further interact as it expands toward the detector elements. From CERES/SPS data [17] there is evidence that the in-medium mass-width (spectral function) of vector mesons may change drastically [18].
Hard probes (jets, leptons, photons, heavy quarks) are of special interest because they provide "tomographic" tools with which can map out this evolution experimentally. Hard probes are effective "external" tomographic probes because they are produced with a pQCD computable initial distribution on a much shorter time scale, ∼ 1/m ⊥ , then the plasma formation time, ∼ 1/3T . Modification of their known initial distributions therefore provides information on the medium through which they propagate in analogy to conventional X-ray or positron tomography used in medicine [19]. The primary advantage of RHIC over lower energy machines (AGS, SPS) is that hard pQCD probes are produced at RHIC orders of magnitude more abundantly over a significantly larger kinematic range. This greatly improves their tomographic resolution power. Figure 1 shows the rapid growth of high p T Au + Au → π 0 + X predicted by pQCD from SPS to RHIC and LHC. As discussed for example in [20,21], it is useful to decompose the nuclear geometry dependence of invariant hadron distributions produced in A + B → h + X at impact parameter b into a phenomenological "soft" and perturbative QCD calculable "hard" components where N part (b) is the number of nucleon participants and N coll (b) = σ pp in T AB (b) is the number of binary N N collisions at impact parameter b. The nuclear geometry of hard collisions is expressed in terms of the Glauber profile density per unity area Fig.2). The hard part scales with the number of binary collisions ∝ A 4/3 because their probability is small built up from all possible independent parton scattering processes. The soft part scales with only N part ∝ A 1 because their probability is large and therefore "shadowed". Fig. 2. Illustration of key aspects of the relation between the geometry of nuclear collisions and the participant and collision number at a fixed impact parameter. The observables (see [3,4,5,6]) used to constrain the geometry experimentally are also illustrated.
The (textbook) computable lowest order pQCD differential cross section for inclusive p + p → h + X invariant cross section is given by where x a = p a /P A , x b = p b /P B are the initial momentum fractions carried by the interacting partons, z c = p h /p c is the momentum fraction carried by the final observable hadron, f α/p (x α , Q 2 α ) is the proton structure function for parton of flavor α, and D h/c (z c , Q 2 c ) is the fragmentation function for the parton of flavor c into h. The UA1 data on pp hadron production with p T > 2 GeV can be well reproduced with this pQCD model expression.
The soft (p T < p 0 ∼ 2 GeV/c) nonperturbative contribution to the hadron yields can only be modeled phenomenological. The Dual Parton Model [22,23] and the LUND string model [24,25] are the most extensive and successful low p T multiparticle phenomenologies. The basic pQCD matrix elements have been encoded into a Monte Carlo code, PYTHIA [26]. A variant of soft string phenomenology tuned to pp, pp data, with the hard part taken from PYTHIA, a hadronization scheme taken from the LUND JETSET hadronization, and a eikonal nuclear multiple collision geometry were combined into the Monte Carlo A+B collision generator in HIJING [27]. HIJING has been used over the past decade to predict many observables at RHIC [20,27,28]. The separate soft and hard components in HIJING with a fixed A, √ s independent scale p 0 = 2 GeV/c are illustrated in Fig.3  The physics in A + A reactions that must be understood in order to be able to interpret observables in terms of the properties of dense QCD matter requires extending the above class of event generator models to include 1. Initial Conditions: The formation physics responsible for creating an incoherent gas of gluons and quarks from the initial virtual nuclear fields, 2. Parton Transport: The (x µ , p µ ) phase space evolution of that parton gas toward equilibrium, 3. Hadronization: The dynamical mechanisms that convert the parton degrees of freedom in the confining physical vacuum into the observable hadronic ones, 4. Hadron Transport: The final state interactions of the expanding dense hadronic matter prior to "freeze-out".
Each problem is fascinating in its own right but only bits and pieces are understood or phenomenologically mapped out up to now. There exists unfortunately no complete computable dynamical theory (like Magneto-Hydrodynamics for QED plasmas) that consistently takes into account all four elements. QCD is believed to be THE theory, but it is still not computable except at high p T where perturbative or classical methods may apply. There exists instead several different dynamical "scenarios" to describe A + A that attempt to patch together different approximation techniques and phenomenological models to address all the physics issues in turn. Two generic approaches to A + A can be classified by whether the Initial Conditions are (1) computed (via pQCD or classical Yang Mills (cYM)) and subsequent evolution followed by a dynamical scenario for 2-4, or (2) the initial conditions are fit by extrapolating final observables backwards via a suitable dynamical scenario. At lower energies (AGS,SPS) only the second approach is available since the momentum scales are simply too low to apply either pQCD or cYM. At collider energies RHIC and beyond, the copious production of mini-jets [20,27,34,35,36] with p T > p 0 ∼ 2GeV shown in Fig.1 makes it possible for the first time to pursue the first approach via pQCD Eq.(2). At very high energies classical Yang Mills theory [37,38,39,40,41,42] provides a general method to compute the Formation Physics which reduces to pQCD at high p T . Whether RHIC or LHC energies are high enough is an open question.
The second approach, trying to "fit" the initial conditions by extrapolating the final distributions backwards with a suitable dynamical model has been traditionally based on relativistic hydrodynamics [43,44,45,46]. The approximate longitudinal boost invariant boundary conditions at ultra-relativistic energies simplify hydrodynamic equations greatly as pointed out by Bjorken [48]. For µ B = 0 the hydrodynamic equations are, where ǫ(x), P (x) are the proper energy density and pressure and u µ (x) is the four velocity field of the fluid. The central assumption is that thermal and chemical equilibrium are maintained locally in spite of the possible large gradients in the fluid variables. The great advantage of hydrodynamics is that it provides a covariant dynamics depending only on the equation of state P (T (x)) that is directly related to the lQCD predictions. When a specific space-time freeze-out hypersurface is assumed together with the assumption, the Cooper-Frye prescription [49,50,51,52,54], the computed four fluid velocity field can be used to predict the final anisotropic flow pattern of hadrons. Since this process is assumed to be reversible, the final distributions together with an assumed freeze-out hypersurface can be used to compute the initial conditions on any desired initial hypersurface. The disadvantage of this approach as emphasized in [55] is that both the initial and final freeze-out hypersurfaces must be guessed. Also finite mean free path physics is outside the scope of ideal hydrodynamics, and transport theory solutions [55] do not support "sharp" freeze-out hypersurfaces. Thus the inversion of data in this way to deduce the initial conditions is not unique. The neglect of dissipative effects such a viscocity also makes it impossible to relate central A + A to peripheral and light ion data, especially p + p. Finally, the assumption of homogeneous or slowly spatially varying initial conditions is questionable because of copious mini-jet production [56]. In spite of all the above theoretical problems, initial conditions for RHIC have been successfully constructed that lead via ideal hydrodynamics and idealized Cooper-Frye freeze-out to distributions that reproduce amazingly well many of the low p T observables at RHIC [57,58,59,60] (see next section).
In order to bring the freeze-out assumption under better theoretical control covariant, nonequilibrium transport theory [61] must be solved. Until recently, only simplified 1+1D Bjorken transport theory was soluble in the linearized relaxation time approximation (see [62] and refs therein). This is due to the great numerical complexity of the 3+1D nonlinear Boltzmann equations [55,63]: where W is the square of the 2 → 2 scattering matrix element, the integrals are shorthands for where g is the number of internal degrees of freedom, while f j ≡ f (x, p j ) is the parton phase space distribution. The initial conditions are specified by a source function S(x, p) that corresponds to the assumed initial conditions. Yang Pang's parton subdivision technique [63,64] and the speed of current workstations have finally made it possible to solve eq.(3) numerically. (Codes can be obtained from the OSCAR Web site [65]). The solutions [66,67] prove that elastic parton scattering with pQCD rates is insufficient at RHIC to keep the plasma in local equilibrium due to the extreme rapid longitudinal "Hubble" expansion of the system [62]. Unfortunately, there exists no practical algorithm at this time to solve the more nonlinear inelastic transport equations involving gg → ng processes. Therefore, if hydrodynamics applies to A+ A at RHIC, then most likely strong nonperturbative mechanisms must be assumed to exist on faith or hypothesis (see E. Shuryak lectures). This is an important open theoretical problem.
I would also like to call attention to a new class of hydrodynamic model [68] that side-step the final freeze-out problem by assuming that local equilibrium is maintained only up to an intermediate hyper-surface, just after hadronization on a T = T c − ǫ isotherm. Using that intermediate freeze-out as the initial conditions of a hadronic transport theory, the subsequent evolution of the hadronic gas toward a dynamical freeze-out is then determined by known hadronic cross sections via URQMD [69,70].

Global Constraints on Initial Conditions
The first results from RHIC from PHOBOS [5], shown in Fig.4, demonstrate that the energy dependence of the scaled charged particle (pseudo)rapidity density, dN ch /dη/N part , is different from p + p and p +p systematics. Measured pseudorapidity density normalized per participant pair for central Au+Au collisions (PHOBOS [5,29]). Systematic errors are shown as shaded areas. Also shown are results of Pb+Pb data (CERN SPS), HIJING [20] simulations and a parameterization of pp data. Normalized pseudorapidity density (PHOBOS and PHENIX) [33] dN ch /dη| |η|<1 /(0.5 × Npart) as a function of the number of participants. Predictions based on HIJING (thick solid) and EKRT [36] (thin) are shown.
Approximately 50% more particles are produced at mid rapidity per participating baryon in central Au+Au collisions then in p+p at the same energy per baryon. The curve shows that the two component HIJING model predicted well this result. However, as shown in Fig.5, another model EKRT [36], was also found to predict the same multiplicity as HIJING for central collisions. In ref. [20], we proposed that the centrality dependence of this observable could differentiate between these competing models of the initial conditions. The new data of PHENIX and PHOBOS [33] verified this prediction. While neither model accounts quantitatively for the data, but the two component HIJING model with its combined A 1 and A 4/3 dependence better describes the rate of increase of the scaled multiplicity with participant number. The observed increase of the scaled multiplicity with energy relative to p + P and with participant number dependence is supports the prediction of copious mini-jet production at RHIC. This is one of the necessary, though insufficient, conditions to form a dense gluon plasma in A + A.
The difference between HIJING and EKRT is that in the latter it is assume that all the produced entropy (multiplicity) arises at RHIC energies from hard pQCD processes. EKRT assume that there is no significant soft component, i.e. dN sof t ≪ dN hard in eq.(1). However, the hard component is cutoff at scale p 0 that is allowed to vary with both A and hence N part and with energy √ s based on the following assumption: independent and hence T AB (b) proportional number of gluons with p T > p 0 are produced only in "resolvable" domains of finite area π/p 2 0 . There are p 2 0 R 2 such domains in the transverse plane in a central nuclear collision. This so called "final state saturation" model is then specified by For β = 1 assumed in EKRT, the solution for the saturation scale is 19 . This predicts dN g /dy ∝ A 0.93 in spite of the apparent proportionality of hard processes to A 4/2 . The flat (dN g /dy)/N part ∼ A ∼0 independence of the scaled multiplicity is a general feature of saturating QCD models of the initial conditions (see also lectures of Mclerran). Such a flat behavior is however ruled out by the present data at RHIC.
An alternate (so-called initial state saturation) model was proposed by KN [71] based on the nonlinear QCD evolution equations of [72]. In this model of nuclear initial conditions, the number of liberated gluons is proportional to the number of virtual gluons participating in the reaction on a scale p 0 . The produced number is then taken to f N part xG(x, p 0 ) in terms of the nucleon gluon structure function, where f ∼ 1.2 is a factor on the order of unity. Since the interaction probability is proportional to the running coupling α s (p 0 ), the initial state saturation condition is defined by The main difference between initial and final state saturation models is therefore due to the logarithmic dependence on p 0 introduced by the running coupling. In [71] a simple ansatz was assumed for xG(x, Q) ∝ logQ/Λ based on the linear (DGLAP) evolution equations. With this ansatz KN predicted a participant dependence surprisingly close to the observed data in Fig.5. However, the x independent ansatz of KN used for xG(x, Q) for the scale Q ∼ 1 GeV/c is a guess that cannot be supported by the available ep HERA data. At small x ∼ 0.01 and low Q ∼ 1 the pQCD factorization analysis of deep inelastic e+p reactions breaks down and xG acquires a 100% systematic error bar as shown in Fig.6. Initial state saturation is a theoretically sound model only at very high energies or nuclei with A ≫ 200, when Q > 2 GeV/c and the errors based on pQCD analysis become manageable. While it is premature to conclude which approach is least wrong (see also [83]), in my opinion, it appears that the global multiplicity data and its centrality dependence can be used as indicators that the initial gluon rapidity density at RHIC is between HIJING's 200 and EKRT's 1000. The corresponding gluon density, ρ g (τ ) = dN g /dy/(τ πR 2 ) is thus ∼ 10 − 50/fm 3 at the corresponding formation time 1/p 0 = 0.1 − 0.2 fm/c. Thus RHIC may have indeed created the densest gluon plasma ever in the laboratory. As I emphasize in a later section, fortunately there are many other observables, especially jet quenching, that provide independent checks of this possibility.
It is important to emphasize that similar results for the multiplicity in central collisions in HIJING and EKRT are purely coincidental because the models differ by a factor of five on the initial gluon density. This is compensated for by the underlying very different hadronization schemes assumed. HIJING creates a large fraction of the observed hadrons at RHIC through its soft string fragmentation scheme, while EKRT assume that entropy conservation implies that dN π ≈ dN g . The lack of a detailed hadronization theory can only be overcome phenomenologically by testing experimentally all the ramifications any particular model.
Another observable that was suggested in [20] to help differentiate models of initial conditions is the shape and scaling of the whole rapidity distribution (see Fig.7). It is seen that HIJING predicts a somewhat narrower and stronger centrality dependence than observed by PHOBOS. This may be related to the baryon stopping power at RHIC. Unfortunately no predictions are available for either the initial or final saturation models on the predicted shape of the rapidity distribution. This observable is especially sensitive in those models to the x dependence of the saturation criteria. The total integrate charge particle multiplicity is shown in Fig. 8. RHIC has produced about 4000 charged particles in Au+Au at 130 AGeV. The nonlinear enhancement near central collisions is interpreted in terms of HIJING as due to the onset of the mini-jet component.

Global Barometric Observable E T /N ch
An important global barometric measure of the internal pressure in the ultradense matter produced is the average transverse energy per charged particle. PHENIX data are shown in Fig.9 compared to WA98 data from CERN. What is most amazing is that E T /N ch ≈ 0.8GeV almost independent of √ s from 20 to 130 AGeV and independent of centrality! HIJING predicts that it should rise from 0.8 to 0.9 GeV from CERN to SPS due to the enhanced mini-jet activity at RHIC. The EKRT initial state saturation model predicts a growth of this quantity in the initial state by about a factor of 3. The reason that EKRT remains viable after these data is that the assumed entropy conservation implies that a large amount of pdV work due to longitudinal expansion is performed by the plasma. In 1+1D hydrodynamics the energy per particle ǫ/ρ ≈ 2.7T decrease as the system expands and cools T ∼ 1/τ 1/3 . if the freeze-out is assumed to occur at all energies and impact parameters in A + A on a fixed decoupling isotherm, then the energy per particle will always be the same. At RHIC this global transverse energy loss from the initial state is predicted to be about a factor 3. The theoretical problem of justifying hydrodynamics and the freeze-out prescription itself discussed in the previous section comes back to haunt us here [57]. The observed NULL effect in E T /N ch is very interesting because it is so difficult to obtain in any transport theory with finite pQCD relaxation rates.

Discovery of Jet Quenching
One of the predicted [28,73,74] signatures of dense matter formation is the suppression of jets and their high p T hadronic debris due to energy loss of the jet in the medium. However, the search for this effect at SPS by WA98 yielded the opposite result as shown in Fig.10. Even a modest dEdx = 0.2 GeV/fm is completely ruled out by the data [81]. The problem is that at lower energies, multiple initial state elastic scattering leads to a random walk in transverse momentum. This enhances the p T of the scattered partons so that This so-called Cronin effect has been well studied in p + A reactions up to 800 GeV. At lower energies the very steep fall of the high p T tail makes the distribution extremely sensitive to this modest p T enhancement. When convoluted through two nuclei, Wang predicted[81] that the Cronin enhancement at SPS in P b+P b should be a factor of two as verified in Fig.11. What is plotted there is the ratio of the observed invariant cross section to the scaled binary collision number, N coll (b), scaled invariant cross section in p+p. Unity corresponds to naive superposition of N coll independent elementary p + p hard processes in the absence of any nuclear effects. The ratio starts below 1 since the low p T distribution grows only with the number of participants (divided by two) and N part (0)/2N coll (0) ≈ 0.15.
In stark contrast to the SPS enhancement of high p T pions, a factor of two or more suppression of p T > 2 GeV hadrons was reported by STAR [3,84] and PHENIX [85].   12 shows that for p T < 2 GeV a similar trend of increase due to the gradual transfer from participant to binary scaling is taking place as at SPS,  Fig. 11. The nuclear modification factor for hadron spectra in central P b + P b collisions at the CERN-SPS exceeds unity at high pT due to the Cronin effect. The solid line is a pQCD calculation by Wang [83].  [3,84]. In contrast to SPS, the high pT charged hadrons are suppressed.
but for p T > 2 GeV the ratio for charged particles π ± + K ± + p ± starts to drop again and reaches ∼ 0.5 at 4 GeV/c.
The PHENIX data [85] shows an even more dramatic quenching pattern for identified π 0 in Fig.13. In this experiment, it was further verified that "peripheral" collisions are not quenched while central ones are. Fig.14 shows that the suppression factor may reach a factor of three at 3 GeV/c. In this plot the ratio is not relative to pp data extrapolated to 130 GeV, but to "peripheral" collisions where the average number of participants and binary collisions is only ≈ 20. In contrast N part ≈ 360 and N coll ≈ 857 for the central collisions. It must be emphasized that current systematic errors are still much larger than statistical, but it is clear that the combined information from two independent experiments in Figs. (12,14) imply that something new has been discovered in A + A collisions at RHIC. I believe that this is the predicted jet quenching as discussed in the next lecture.
The reason that this discovery is perhaps even more exciting than the famous J/ψ suppression effect discovered by NA50 [86,87] at the SPS is that J/ψ suppression was also seen in p + A. The cold nuclear suppression mechanism in p + A is called "normal". The enhance suppression in P b + P b is "anomalous" because it is more than if the normal p + A suppression pattern is extrapolated to A + A That this is not theoretically fool proof was pointed out by Qiu et al [88]. They showed that including radiative energy loss in cold  Fig. 13. Semi-inclusive π 0 pT distribution (1/Nint)(dN π 0 /2πpT dpT dy)) in the upper 60-80% peripheral events (solid circles) and the 10% most central events (solid squares) from PHENIX [85]. The lines are a parameterization of pp charged hadron spectra, scaled by the mean number of collisions N coll = 857, 19 resp. The bands indicate the possible range due to the systematic error on N coll .  Fig.11) is due to jet quenching in the ultra-dense matter formed at RHIC [21,27,73,28,75,83,89] .
nuclei could lead to non-linear enhancement of J/ψ suppression by decreasing their formation probability. Only a rather schematic model was presented, but it emphasizes the necessity of improving considerably the theory of the "normal" processes associated with heavy quark propagation through with nuclei. The situation is rather similar theoretically with regard to the Cronin effect. There also only rather schematic models are available to simulate the effect.
The big difference between the two cases is that for J/ψ the "normal" and "anomalous" components work in the same direction. The premium is thus high on developing an accurate theory "normal" nuclear suppression. In the jet quenching case, on the other hand, the "normal" Cronin effect works in the opposite direction to the "anomalous" new jet quenching mechanism. Of course, there are possibly other "normal" effects, such as gluon (anti?)shadowing, that may work in either direction at high p T . To map out all the "normal" physics components will require detailed systematic measurements of p + A at RHIC as done at the SPS. As a final remark, I want to emphasize the "normal" component of the dynamics is not dull run-of-the-mill background, but fundamentally interesting many body QCD physics in its own right and deserves considerable more attention.

Where Have all the Baryons Gone?
One of the puzzling feature of Figs.(12,14) is that pions appear to be more quenched than the sum of charged particles. Usually we assume that pions are the most abundant hadron species at high p T since both quark and gluon fragmentation functions prefer to make the lightest mesons [91,92]. Surprisingly, the preliminary PHENIX data [90] on identified high p T hadron spectra suggest from Fig.15 that baryons may be the most abundant species above p T > 2 GeV/c. One possible source of such a non-pQCD like flavor distribu- tion could be hydrodynamic transverse flow. For a longitudinal boost invariant (Bjorken) expansion with a transverse flow velocity field, v ⊥ = tanh(η r ), the general formula [93] for the differential particle number is where d = 2s + 1 is the degeneracy factor, η r = Artanh(v ⊥ (z = 0)) is the transverse fluid rapidity and (r f (ζ), τ f (ζ)) is a parameterization (counterclockwise) of the freeze-out surface (isotherm of temperature T f ).
Solutions for freeze-out surfaces with arbitrary transverse velocity fields v ⊥ (ξ) can be obtained by solving relativistic hydrodynamics. For the simplest case with v ⊥ = tanh η r a constant and an isotherm freeze-out on a proper time hypersurface τ f [94], which corresponds to a blue shifted effective temperature T f e ηr . This is the uniform rapidity, transverse boosted Bjorken sausage parameterization of nuclear collision distributions. Evidence for increased transverse flow phenomena at RHIC relative to SPS comes from low p T STAR data [99] shown in Figs. (16,17). The data can be fit up to p T < 1 GeV/c with a rather radial flow velocity v ⊥ ∼ 0.6 c that is significantly larger than the radial flow ∼ 0.4 c deduced from similar SPS spectra.
Another important experimental tool to search for collective flow effects is to study anisotropic multiparticle emission patterns [95,96,97,98]. A particularly useful measure of collective behavior in ultra-relativistic energies has turned out to be the differential second Fourier component [96] of the azimuthal distribution: where φ is measured relative to a "reaction plane" event by event as determined in [98]. Azimuthal or "elliptic" flow results from the initial spatial anisotropy of the dense matter in semi-peripheral A + A collisions. The hydrodynamic model predicts an elliptic flow pattern at RHIC [58,59] that can be approximately parameterized as v s 2 (p T ) ≈ tanh(p T /(10 ± 2 GeV)) .
Up to about p T < 1 GeV, this agrees remarkably well with STAR data. At high p T this hydrodynamic flow component breaks down because of the emergence of the hard pQCD hadrons. The transverse boosted Bjorken sausage model eq.(7) predicts that asymptotically the baryon/meson ratios p/π+ =p/π − → 2 for any flow velocity  because nucleons have 2 spin states. However, numerically this ratios exceed unity only p T > 3 − 4 GeV. Thus transverse flow alone cannot account for the anomalous baryon dominance of high p T spectra in Fig.15 as emphasized in [92].
Another observation [99] that possibly provides a hint that the answer to the puzzling result may lie in novel baryon dynamics at RHIC can be seen in Figs. (18,19). As was shown by Kharzeev [101], the energy and rapidity dependence of the inclusive baryon production at mid-rapidity can be obtained using Mueller's generalized optical theorem in the double Regge limit. Here, the exchanges of a Pomeron and a M J 0 Baryon-anti-baryon "junction" pair lead to the following form for single mid-rapidity baryon production, where C B is a constant that reflects the couplings of the Reggeon and Pomeron to the proton, f B (m 2 t ) is an unknown function of m t and s 0 ≃ 1 GeV is a Regge energy scale. The cosh(y/2) rapidity dependence and the 1/ 4 √ s energy dependence follow from the assumed intercept [102], α M J 0 (0) ≈ 1/2. In contrast to simpler diquark breaking models as in the Dual Parton Model, the multiplicity of junction also enhanced events is enhanced by a factor of 5/4   in p + p, and the strangeness content is also enhanced by a large factor. The junction mechanism for baryon number (vs valence quark number) transport predicts for the unique possibility of producing S = −3 Ω − baryons at midrapidity, as were observed at the SPS in WA97. In the Monte Carlo event generator HIJING/BB [65,100], baryon junctions are implemented in terms of Y shaped strings spanning valence quarks. The junction is a topological knot in the gluon field connecting the color flux from three quarks into a color singlet state [102]. The intriguing aspect of junctions is that the conserved baryon number resides in the gluon knot and not in the valence quarks [101]. In a nuclear collision some or all of the valence quarks may fragment into mesons. However, the gluonic junctions insure that baryon number is conserved. The understanding of the dynamics of junction exchange and pair production is still rather primitive, but the consistency of the baryon stopping power at SPS and now RHIC with HIJING/BB predictions suggest that baryon dynamics at central rapidities may be especially interesting at RHIC. See ref. [103] for a discussion of possible novel junction network physics that may lead to femto-scale buckyball and even CP odd junction network production in A + A.

Quenching of Elliptic Flow
As seen in Figs. (17,20) strong elliptic flow was discovery at RHIC consistent with hydrodynamic prediction at low transverse momentum p T < 1 GeV. However, the preliminary data from STAR [104] shows that above p T > 2 GeV the elliptic flow saturates and the azimuthal asymmetry deviates more and more from hydrodynamic behavior as seen in Fig.21. This information  provides insight into how hydrodynamic behavior breaks down at high p T due to the finite energy loss of partons in the plasma. As shown in detail in [21] the saturation pattern at high p T depends on the energy dependence of the gluon energy loss as well as on the geometry of the plasma density at finite impact parameters. It therefore provides tomographic information about the density profile and its evolution in A + A. See sec 4.4 for more details.
3.6 Where did the slowly burning plasma log vanish?
The last major RHIC result that I highlight here is on pion interferometry. Relativistic combustion theory [105,106,93] predicts that if there were a sufficiently rapid cross over between the QGP and hadronic phases of ultra-dense matter, then a deflagration burn front may appear between two phases. The main characteristic of that burn front is its very small velocity in case the entropy density jump across it is sufficiently large and no high degree of nonequilibrium supercooling arises. Even with a smooth cross-over transition, such slowly burning plasma solutions were shown to exists as long as the width of the transition region is (∆T c /T c < 0.08). The lifetime of a Bjorken plasma log is therefore significantly enhanced τ ∼ R/v d , where v d ∼ 1/25 is the small deflagration velocity in the static 1+1D case. This characteristic time delay of the hadronization from a QGP state was suggested in [107,108] to be testable via pion interferometry. In ref. [93] the 3+1D hydrodynamic equations were solved to study this plasma "stall" phenomenon in detail.
The two pion correlation function measures the coincidence probability P (p 1 , p 2 ) of two (identical) bosons with momenta p 1 , p 2 relative to the probability of detecting uncorrelated particles from different events, If the average 4-momentum is denoted as K µ = (p µ 1 + p µ 2 )/2 and the relative 4-momentum as q µ = p µ 1 − p µ 2 , then under the assumption that the particle source is chaotic and sufficiently large, where [47] is the single inclusive momentum distribution, f (x) = (e x − 1) −1 , and u µ the fluid 4-velocity. The integrals run over the assumed freeze-out hypersurface. In general, that hypersurface is represented by a 3-parametric (4-vector) function Σ µ (ζ, η, φ), and the normal vector on the hypersurface is determined by where ǫ µαβγ is the completely antisymmetric 4 tensor. For the common isotherm freeze-out temperature T f hypersurface, the fluid velocity gerenally varies u µ = u µ (Σ). For the Bjorken cylinder geometry, it is useful to restrict consideration to particles emitted at midrapidity, K z = q z = 0. Rotational symmetry around the z-axis in central collisions makes it possible to choose the average transverse momentum as K ⊥ = (K, 0, 0), and consequently, C 2 (K, q out , q side ) is a function of three independent variables only. The so called out and side projections of the relative momenta are q out = (q out , 0, 0), q side = (0, q side , 0). As shown in [107,108] the width, 1/R side , of the correlation function in q side is a measure of the transverse decoupling or freeze-out radius, while the width 1/R out of the q out correlation function is also sensitive to the duration of hadronization, ∆τ R 2 out ≈ R 2 side + v 2 ∆τ 2 . Thus a QGP stall would manifest experimentally in R out ≫ R side . In [93] it was found that for possibly realistic parameters, R out /R side ∼ 2 − 3, could be observed if a QGP stall occurred. Preliminary PHENIX pion interferometry data [109] vs different projections of the relative momenta. Similar preliminary STAR data [3] were also shown at QM01. Unlike predictions [93] there is no hint of the expected stall or time delay of the QGP transition.

Au+Au -2π
With this "warm-up" review of pion interferometry, we are now ready for the rude awakening from ideal gedanken considerations with the first splash of "cold" RHIC pion interferometry data shown in Fig.(22). The preliminary PHENIX data show that R out ∼ R side and even more disturbing all the deduced interferometry parameters are virtually identical to values seen at the AGS and SPS. To add insult to injury, it appears that R out < R side for p T > 0.4 GeV. Preliminary STAR data [3] show the same tendency.
Of course scenarios may be invented to "explain" the data a postiori, but if these data are confirmed by further measurements, then they are indeed surprising and call into question our picture of the space-time evolution of A + A. That this problem is not restricted to idealized hydrodynamics can be seen from the results of Ref. [110]. It was shown that R out > R side is also predicted in a calculation where the entropy jump is small and pion decoupling is dynamically handled via URQMD. Among the theoretical questions that should now be further investigated is whether the pion interferometry theory based on chaotic ensembles [111] is in fact applicable to A + A. Another question that needs further study is whether the assumed ensemble of initial conditions was too restrictive and whether highly inhomogeneous and turbulent initial conditions apply [56].

Jet Quenching: Theory
Having had a brief tour of some of the interesting new data harvested from RHIC during the first round of experiments, I turn next to the more specific theoretical problem of computing the energy loss per unit length of a fast parton penetrating a finite, expanding quark-gluon plasma. As I emphasized above, high p T many body pQCD physics is a new frontier at RHIC and higher energies. This requires the development of the non-abelian analogue of the radiative energy loss theory familiar from classical E&M. The interesting new twist is that we have no external beams of quarks or gluons and the medium is very thin due to the fact that nuclei are tiny. Also the formation time physics of Landau-Pomeronchuk-Migdal (LPM) results in major destructive interference effects that must be taken into account. Work on this problem over the past five years has advanced considerably but many open problems remain.
I will only highlight only one of those direction, namely the opacity expansion reaction operator method that we developed in refs. [75]. The reader is referred to BDMS [77], Z [78], and U [79] for alternative methods and approximations.
In ref. [74] we proposed a simple model to study induced gluon radiation due to multiple elastic scattering of a high energy jet propagating in a locally color neutral amorphous plasma.

GLV Formalism
In [76] we developed a systematic graphical method to compute medium induced gluon radiation amplitudes as shown in Fig.23. The exponential growth of the number of graphs with the number of interactions makes it very tedious to go beyond order three. In GLV [75] we overcame the combinatorial problem by developing a new algebraic operator technique to solve the inclusive radiated gluon distributions recursively. The first step is to compute the three direct (single Born) and four surviving virtual (contact double Born) diagrams shown in Fig.24.
For scattering off of n scattering centers located at depths z i in a transverse homogeneous medium of large area ((µR) 2 ≫ 1), we can write the inclusive radiated gluon spectrum, P n (k, c), as a sum over products of partial sums of amplitudes and complementary complex conjugate amplitudes. Every term in the sum contributes to the same O(g 2n ). The average value of σ=(0,0,1,0,1) → l=(0 2 + 0 2 + 1 2 + 0 2 +1 2 ) ⁄ 2 2 1 2 3 4 5 Fig. 23. Induced radiation amplitude [76] contributing to fifth order and higher order in the opacity expansion of QCD energy loss in the GW model [74]. The crosses denote static color screened Yukawa interactions on a scale µ. The blob is the initial hard jet amplitude without final state interactions. 24. Three first order (singe Born) direct and four surviving (double Born) virtual or contact amplitudes [77] from which theDn andVn components of the reaction operator in eq.(16) are derived in [75]. n is referred to as the opacity of the medium. The partial sums of diagrams at order n in such and opacity expansion can be conveniently expressed in a tensor notation and constructed by repeated operations of,1,D i , orV i corresponding to no, direct, or virtual interactions at scattering center i Here G 0 is the initial hard q + g color matrix amplitude. In the inclusive probability each class contracts with a unique complementary class with the complementary class constructed as Fig. (25) shows an example of how this formalism works at 4th order in opacity for elastic and inelastic inclusive distributions. Direct interactions enlarge rank n − 1 class elements as follows: is the so-called Bertsch-Gunion amplitude for producing a gluon with transverse momentum k in an isolated single collision with scattering center n. The momentum transfer to the jet is q n . The notation ω n = (k − q n ) 2 /2ω, for a gluon with energy ω and a n is the color matrix in the d R dimensional representation of the jet with color Casimir C R . N v = n−1 m=1 δ im,2 counts the number of virtual interactions in A i1···in−1 .
Unitarity (virtual forward scattering) corrections to the direct processes involve the sum of four double born contact diagrams in Fig.(24) that enlarge rank n − 1 classes as follows: This key operator relation between direct and virtual insertions that we discovered in [75] makes it possible to solve the problem algebraically. The tensor classification of classes of diagrams makes it possible to construct the distribution of radiated gluons in the case of n interactions, P n , recursively from lower rank (opacity) classes via a "reaction" operator Using the key identity (15), the reaction matrix simplifies tô R n = (D n − a n ) † (D n − a n ) − C A = (Ŝ n +B n ) † (Ŝ n +B n ) − C A The next major simplification occurs because bothŜ andB involve the same gluon color rotation through if cand . This reduces the color algebra to multiplicative Casimir factors Re e −iωnzn e iqn·b I n−1 I n obeys a recursion relation from which the inclusive radiation probability is found to obey the soluble recursion relation whereb = i∇ k is the transverse momentum shift operator. The initial condition for this recursion relation is the initial hard vertex radiation amplitude without final state interactions that is given by P 0 = C R H 2 = C R /k 2 ⊥ . The complete solution to the problem can therefore be expressed in closed form as This expression can be averaged over any spatial distribution of interaction centers, z i as well as any z i dependent momentum transfers q n . This form is thus ideally suitable for Monte Carlo implementation for arbitrary q i , z i medium ensemble averages.

Nonabelian Energy Loss at Finite Opacity
The first application [75] of our general solution to the energy loss problem was to calculate numerically the total radiated energy loss as a function of jet energy E, plasma depth, L, and infrared screening scale, µ. In the absence of a medium, the gluon radiation associated with a spin 1 2 parton jet is distributed as where x = k + /E + ≈ ω/E, and C R is the Casimir of the (spin 1/2) jet in the d R dimensional color representation. The differential energy distribution outside a cone defined by k 2 ⊥ > µ 2 is given by where the upper kinematic limit is k 2 .. The energy loss outside the cone in the vacuum is then given by While this overestimates the radiative energy loss in the vacuum (self-quenching), it is important to note that ∆E (0) /E ∼ 50% is typically much larger than the medium induced energy loss. Averaging over the momentum transfer q 1⊥ via the color Yukawa potential leads to a very simple first order opacity result for the x ≪ 1 gluon double differential distribution where the opacity factor L/λ g = N σ (g) el /A ⊥ arises from the sum over the N distinct targets. Note that the radiated gluon mean free path λ g = (C A /C R )λ appears rather than the jet mean free path. The upper kinematic bound on the momentum transfer q 2 max = s/4 ≃ 3Eµ, (1/µ 2 ef f = 1/µ 2 − 1/(µ 2 + q 2 max )). For SPS and RHIC energies, this finite limit cannot be ignored as we show below.
Numerical results comparing the first three orders in opacity corrections to the hard distribution Eq. (17) were presented in [75]. To illustrate the result consider a quark jet in a medium with λ g = 1 fm, a screening scale µ = 0.5 GeV and α s = 0.3 . The total radiative energy loss could be expressed as with N (∞) = 4 log(E/µ)/ṽ if the kinematic bounds were ignored as in the approximations of ref. [77]. We found that finite kinematic constraints and the form of the first order result cause N (E) to deviate considerably from the asymptotic value for all energies accessible in the RHIC range. Together with the logarithmic dependence on energy, these kinematic effects suppress greatly the energy loss at lower (SPS) energies as seen in Fig.26. This is in sharp contrast to the approximately energy independent result in BDMS-ZW where the finite kinematic bounds were neglected because only the asymptotic  [77] is shown for comparison. The energy dependence of GLV suppressing radiative energy loss of low energy jets explains why no jet quenching was observed at the SPS (see Fig.10) .
limits were considered. Another remarkable result demonstrated numerically is that the second and third order contributions to the integrated energy loss remains surprisingly small in the physical range of nuclear opacities L/λ g ∼ 5.
The rapid convergence of the opacity expansion even for realistic opacities results from the fact that the effective expansion parameter is actually the product of the opacity and the gluon formation probability Lµ 2 /2xE.The leading quadratic dependence of the energy loss on nuclear thickness discovered in BDMS [77] therefore already emerges from the dominant first order term in the opacity expansion.
At SPS energies kinematic effects suppress greatly the energy loss relative to BDMS. Our estimates provide a natural explanation for the absence of jet quenching in P b + P b at 160 AGeV observed by WA98. At RHIC energies, on the other hand, a significant nonlinear (in A) pattern of suppression of high p ⊥ hadrons relative to scaled pp data is predicted.

The Opacity of the QGP at RHIC
As a second a application of the GLV energy loss, in ref. [89] we computed the quenched pQCD distribution of high p T hadrons as a function of the effective static plasma opacity, L/λ g . In Figs. (27,28), the jet energy dependence of the GLV energy loss for gluons is shown. The most important feature to note is that ∆E GLV /E is approximately constant in the energy range accessible at RHIC.  In order to compute the pion spectrum, note that jet quenching reduces the energy of the jet before fragmentation. We concentrate on mid-rapidity (y cm = 0), where the jet transverse momentum before fragmentation is shifted by the energy loss as in [112], p * c (L/λ) = p c − ∆E(E, L). This shifts the z c parameter in the fragmentation function of the integrand (23) The invariant cross section of hadron production in central A+A collision is then given by [91] where upper limit of the impact parameter integral is b max = 4.7 fm for 10 % central Au+Au collisions. Here t A (b) is the usual (Glauber) thickness function. The factor z * c /z c appears because of the in-medium modification of the fragmentation function [112]. Thus, the invariant cross section (23) depends on the average opacity or collision number,n = L/λ g . The calculated spectra for pions are displayed forn = 0, 1, 2, 3, 4 in Fig. 29. Fig. 30. shows their ratios to the non-quenched spectra atn = 0. We note that in contrast to previous energy independent estimates for the energy loss, the GLV energydependent energy loss leads to constant suppression of the high p T domain in agreement with the preliminary data. The peripheral collisions are consistent with a rather small opacity in contrast to central collisions, as expected. The ratio of central to peripheral PHENIX [85] data from QM01 shown in Fig. 29. Pion production in Au+Au collision including jet quenching with opacity L/λ = 1, 2, 3, 4. Preliminary QM01 PHENIX data shown (see updated data from [85] in Fig.13)  Fig. 30 clearly reveals that jet quenching at RHIC overcomes the Cronin enhancement at zero (final state) opacity. This is in stark contrast to data at SPS energies, where WA98 found no evidence for quenching in P b + P b at 160 AGeV but a factor of two Cronin enhancement as discussed before.
Figs. 29,30 indicate that an effective static plasma opacity L/λ = 3 − 4 is sufficient to reproduce the preliminary jet quenching pattern observed at RHIC. In ref. [83] it was shown that a rather small constant dE/dx ≈ 0.25 GeV/fm was also found to be consistent with the data. However, it is important to emphasize that these effective static plasma opacities and parameters hide the underlying rapid dilution of the plasma due to expansion. The GLV formalism including the kinematic constraints at first order has been further generalized to include effects of expansion in [21]. It was found in [21] that the inclusion of longitudinal expansion modifies the static plasma results in such a way that the moderate static plasma opacity actually implies that the produced mini-jet plasma rapidity density may have reached dN g /dy ∼ 500.

Jet Tomography from Quenched Elliptic Flow
So far we have not included the dilution effect of expansion on the energy loss. The generalization of GLV to the case of expanding plasmas is [21] where ρ(z) is the plasma density at time z along the jet path at position z from the production point and where the screening scale µ(z) may also depend on time.
Consider a density evolution of the form as in [77], where α = 0 corresponds to a static uniform medium of thickness L, while α = 1 to Bjorken 1+1D longitudinal expansion transverse to the jet propagation axis.
Analytic expressions can only be obtained again for asymptotic jet energies when the kinematic boundaries can be ignored. In that case, all but the path integral can be done giving which is a linear weighed line integral over the local transport coefficient [77] (µ 2 (z)/λ(z)) ≈ 9 2 πα 2 s ρ(z) however enhanced by a log 2E/Lµ 2 (z) factor that results from the structure of the GLV integral missing in the BDMS asymptotic limit. For an expanding plasma as in (25) ∆E α (L, z 0 ) ≈ C R α s 2 µ 2 (L)L α λ(L) Hereṽ = log 2E/Lµ 2 (L) and we used that µ 2 (L)L α /λ(L) is a constant independent of L for this type of expansion. We also took the z 0 → 0 limit. We therefore recover the asymptotic BDMSZ energy loss for both static and expanding media modulated by a log E/ω(L) factor that is important at RHIC energies. Using the Bjorken relation between the gluon density and the rapidity density then gives ∆E α=1 (L) = 9C R πα 3 s 4 1 πR 2 dN g dy L log 2E Lµ 2 (L) .
In practice, it is straight forward to integrate GLV numerically including the finite kinematic constraints. For non-central collisions the GLV line integral depends of course on the azimuthal direction φ of the jet. The variation of the azimuthal energy loss with respect to φ at a given impact parameter b can be expressed in terms of with results shown in Fig.31 The  energy loss is to induce an apparent elliptic flow at high p T not related to hydrodynamic phenomena of low p T . In [21] we proposed a simple interpolation between the hydrodynamic and jet quenched p T eikonal regimes This interpolates between the hydrodynamic and the pQCD regimes because at high p T , dN h ≫ dN s . For our numerical estimates the low p T interpolation is achieved by turning off the pQCD curves with a switch function f h (p T ) = 0.5[1 + tanh(3(p T − 1.5 GeV))]. We see in Fig.32 that the magnitude and shape of the high p T elliptic flow provides a complementary probe of the initial gluon density and is also sensitive to the geometrical distribution of the plasma. The saturated v 2 increases systematically with increasing initial plasma density and thus provides an important complementary constraint on the maximum initial parton density produced in central b = 0 collisions. The consistency of the quenched elliptic flow in non-central with the central quench pattern will be very important to test when the final data become available. Hydro v 2 (p T )=Tanh(p T /12) Fig. 32. The interpolation of v2(pT) between the soft hydrodynamic [58] and hard pQCD regimes [21] is shown for different gluon rapidity densities in central b = 0 collisions. The gluon rapidity density at b = 0 is assumed to scale with the binary collision number. Solid (dashed) curves correspond to sharp cylindrical (diffuse Wood-Saxon) geometry

Summary
If confirmed by further measurements and theoretical refinements, jet quenching may have already provided the first evidence that initial parton densities on the order of 100 times nuclear matter density may have been produced at RHIC. The full analysis of flavor composition, shape, and azimuthal moments of the high p T spectra appears to be a promising diagnostic probe of the evolution of the gluon plasma produced at RHIC. However, it is too early to tell what the preliminary say about the properties of that extremely dense form of matter. There are too many pieces of the puzzle that simply do not fit well into any scenario. The beam energy and centrality independence of the transverse energy per charged particle is one of them. The anomalous baryon number transport to high transverse momenta and central rapidities is another. Finally, the puzzling beam energy independence of the preliminary pion interferometry results is a mystery. As the tera-bytes of RHIC data continue to stream in during the next few years, they will certainly pose many interesting new QCD many body problems. The new chapter on the physics of ultra-dense matter and the dynamics of ultra-relativistic nuclei is now unfolding at RHIC.