Event-by-Event Jet Quenching

High momentum jets and hadrons can be used as probes for the quark gluon plasma (QGP) formed in nuclear collisions at high energies. We investigate the influence of fluctuations in the fireball on jet quenching observables by comparing propagation of light quarks and gluons through averaged, smooth QGP fireballs with event-by-event jet quenching using realistic inhomogeneous fireballs. We find that the transverse momentum and impact parameter dependence of the nuclear modification factor R_AA can be fit well in an event-by-event quenching scenario within experimental errors. However the transport coefficient qhat extracted from fits to the measured nuclear modification factor R_AA in averaged fireballs underestimates the value from event-by-event calculations by up to 50%. On the other hand, after adjusting qhat to fit R_AA in the event-by-event analysis we find residual deviations in the azimuthal asymmetry v_2 and in two-particle correlations, that provide a possible faint signature for a spatial tomography of the fireball.


Introduction
Collisions of nuclei at high energies at the Relativistic Heavy Ion Collider (RHIC) and, soon, at the Large Hadron Collider (LHC) create a fireball with local energy densities well above 1 GeV/fm 3 . At those densities, quarks and gluons form a deconfined quark gluon plasma (QGP) [1]. In some of the collisions, high momentum partons in the initial nuclear wave functions scatter off each other and propagate away from the collision axis. They form large momentum jets in the final state. These jets, and the hadrons fragmenting from them, can be used as hard probes of the fireball. The interactions of the scattered partons with quark gluon plasma lead to radiative energy loss and a significant suppression of the hadron yield at high transverse momentum P T [2,3,4,5,6,7,8]. One of the key results from RHIC was the confirmation of this jet quenching effect: Hadrons with P T ≥ 5 GeV/c are suppressed by about a factor 5. In addition, an extinction of away-side jet correlations has been seen in a certain kinematic regime, further emphasizing the large opacity of quark gluon plasma [1].
Email address: rjfries@comp.tamu.edu (R. J. Fries) The study of quark gluon plasma has now moved into an era of quantitative assessments of experimental results. A simple question that we should be able to answer is that of an averaged value q of the transport coefficientq = µ 2 /λ, the average squared momentum transfer µ 2 of a high momentum parton per mean-free path λ. The averaging . . . here refers to the many possible paths of a parton (thus sampling different localq along the trajectory in a cooling and expanding fireball), an average over parton species (until we have means to reliably distinguish gluon and quark jets), and the average over many event geometries for a given centrality bin (hard processes at RHIC are rare and experimental results are event averaged).
Comparative studies using the mainstream energy loss models lead to a somewhat unsettled picture. Bass et al. [9] have reported a wide range of possible values forq 0 , the local initial value of q at the center of a central collision, ranging from 2.3 GeV 2 /fm to 18.5 GeV 2 /fm depending on the energy loss model, and on how the localq(r, τ ) is modeled as a function of the local energy density or temperature at position r and time τ . In addition, we have been cautioned by results that show thatq extracted from single and di-hadron nuclear suppression factors are not necessarily compatible, and thatq is sensitive to assumptions on pre-equilibrium quenching and the initial parton spectrum [10] as well as radiative corrections to the hard process [11].
Clearly we have to discuss and constrain uncertainties in our modeling of jet quenching very carefully in order to arrive at reliable quantitative estimates. In this work we investigate the influence of inhomogeneities and fluctuations in the fireball on jet quenching observables. As mentioned above, experimental data are averages of observables over many events, where in each event jets are created and propagate through the underlying fireball. Most calculations found in the literature on the other hand turn this process around and propagate jets through an idealized fireball which can be understood as an average over realistic fireballs. These smooth fireballs come in various degrees of sophistication, from a simple overlap of nuclear thickness functions in the transverse plane to the use of detailed maps from hydrodynamics calculations that take into account the proper expansion and cooling. However, even the latter are often based on averaged and hence idealized initial conditions. It has been realized before that event-byevent computations are crucial to understand some low-P T observables, e.g. hydrodynamic elliptic flow [12]. It is important to take into account that the overlap of two nuclei (i) is irregularly shaped, (ii) is generally not aligned with the naive geometrical reaction plane, and (iii) exhibits local fluctuations with hot spots and cooler regions. This leads to appreciable differences compared to computations using averaged, idealized fireballs.
Here we investigate whether an event-by-event computation of jet quenching differs from one using an averaged event. If the answer is yes, an interesting question arises: is it possible to study some features of the spatial structure of fireballs with hard probes, despite the averaging over events? In other words, is a true tomography feasible?

Effects of Inhomogeneities on Quenching
Let us discuss some general expectations when we go from parton propagation through an averaged fireball to an average over propagation in many fireballs. First consider the limit of extreme quenching, qR 2 ≫ p where R is the typical size of the fireball and p the momentum of the final state parton. All observed particles then come from the surface of the fireball. It is clear that we should expect more such particles from an inhomogeneous fireball compared to a smooth fireball with equal total energy, ifq is a fixed function of the fireball density. This is due to the larger effective surface area of an inhomogeneous fireball, see e.g. Fig. 2. Hence the single and double particle nuclear modification factors, should increase for partons, and for hadrons fragmenting from them. We also expect the azimuthal anisotropy i.e. the difference of parton emission out of the reaction plane and into the reaction plane, to decrease since the relative increase in surface should be larger on the out-of-plane side. Please note that the number of collisions N coll in the denominator is an averaged number estimated for the corresponding centrality bin. We will follow this procedure and will not divide by the number of collisions on an event-by-event basis. For completeness let us also give the definition for the nuclear modification factor of the two-particle correlation per trigger I AA (P T 1 , P T 2 ) = J AA (P T 1 , P T 2 )/R AA (P T 1 ) which we will use later. Let us consider a more quantitative example. Imagine energy loss ∆E = Ch β along a parton trajectory determined by an expression of the general type where r and ψ are the point of creation and the emission angle of the parton, e ψ the unit vector along the trajectory, and τ the time elapsed since creation of the parton. β encodes the path-length dependence with linear or quadratic dependence corresponding to β = 0 or β = 1 resp., and C is a coefficient. ρ(r) encodes a local property of the fireball, akin to a density, which we do not specify further at this point. Let n(r) be the probability for a parton to emerge from point r. The relevant quantity to study is the energy loss weighted with the emission probability, n(r)h β (r, ψ). This quantity characterizes the suppression of single particle spectra for not too large quenching (p ≫ ∆E) as we can infer from the following exercise for a powerlaw parton spectrum dN init /d 2 rdψdp ∼ n(r)p −α . Expanding the expression (p + Ch β ) −α in the final parton spectrum for small energy loss we can express the final spectrum as The energy loss model here is deterministic, but we do not expect major modifications of the following arguments if h β only gives an average value of a statistical process of scattering and gluon emission. Now we consider the pair of densities n i (r), ρ i (r) event-by-event by writing them as a sum of ensemble expectation valuesn(r),ρ(r) and fluctuations δn(r), δρ(r), resp. The ensemble average of the single particle suppression is given by The first term is the result from propagating through the averaged fireball, and we have omitted terms linear in fluctuations due to δn = 0, δρ = 0. The last term contains the correction due to fluctuations We have introduced the correlation function between fluctuations of the position of hard collisions and the density of the bulk fireball. Eq. (7) is a rather general statement one can make about event-by-event fluctuations without imposing too many restrictions on the energy loss mechanism. Our result indicates that the leading deviation due to fluctuations is given by correlations between the emission point of the jet and the fireball along its trajectory.
What constraints can be put on R? We expect fluctuations to be granular with a certain length scale σ (e.g. the nucleon diameter if n is related to the density of nucleon-nucleon collisions and ρ to the density of participant nucleons). Then R must be positive for |r 2 − r 1 | σ because on average the density of hard processes and the density of the soft, "underlying" event should be positively correlated. On the other hand, R will turn negative on distance scales larger than σ because the total amount of matter in the transverse plane is conserved on average. In other words, a hot spot of the fireball has to be compensated, on average, by less material around that spot. This is also the reason for the argument of a larger effective surface that we raised for the case of surface-dominated emission: Clumping of density somewhere along the boundary of the fireball introduces "holes" elsewhere. In principle ρ also carries an explicit dependence on the time τ elapsed -suppressed in the notationsince the fireball is evolving dynamically. We expect that relaxation phenomena or hydrodynamic evolution wash out the correlation function R, although this might take several fm/c by which time a large fraction of the observable jet strength has left the fireball. We check our qualitative expectations with an example shown in Fig. 1. We provide two cuts through the correlation function of the densities of binary collisions, ρ = n ≡ n coll , calculated from the initial state simulation GLISSANDO [13].
We discuss more details about GLISSANDO in the next section. We clearly see positive correlations with a radius σ ≈ 1 fm as expected from fluctuations based on collisions of nucleons. The anticorrelation region extends all the way to the point where the nuclear overlap zone ends.
For the following it is useful to look at a simplistic parametrization for the correlation function R which should be qualitatively true for a wide class of models. Let us assume a non-spherical fireball with short and long axes X and Y , resp. In a fluctuating fireball these are only expectation values, of course. Let us further assume where ∆r = |r 2 −r 1 | and we neglect the dependence of R on the center coordinate r 1 + r 2 which should be a satisfying approximation for fireballs which are on average uniform (i.e.ρ andn are smooth) and large, (i.e. X, Y ≫ σ). λ and µ are positive numbers that characterize the correlation strength on distance scales σ and the anti-correlation strength on larger distances, resp. We note that R should go to zero if the relative distance becomes too large which is not duly captured in the ansatz above. However this should not change the following interesting result on elliptic flow. First we note that (10) where l is the length of the parton trajectory in the fireball. We clearly see that the sign of the correction is determined by the competition between an increased suppression coming from more jets being emitted in regions with a denser fireball, and the decreased suppression around those regions.
While it is hard to predict the sign of δ(nh) even after integration over emission points r without any further concrete assumptions we can make the following observation. Let us use the difference of energy loss in-and out-of-plane as a proxy for v 2 . To be more precise, v 2 should be a monotonously rising function of − d 2 r n(r) (h(r, 0) − h(r, π/2)) .
for reasonable β since X < Y . Hence, the azimuthal anisotropy v 2 tends to be diminished in event-by-event calculations for a broad variety of energy loss models. Basically there is more room for the anti-correlation in R to decrease energy loss along the longer side of the fireball than along the narrow side. This is compatible with the argument we made in the case of extremely strong quenching and surface dominated emission, and it can also be seen in Fig. 1.

Numerical Study
We want to back up some of the analytic arguments from the last section through a numerical study. The distribution of hard collisions is usually taken to be the density of binary nucleon-nucleon collisions, n(r) = n coll (r).q(r) is often assumed to be a function of the local energy density ǫ(r) in the transverse plane. Around midrapidity the initial energy density is usually modeled as a superposition of the density of collisions and the density of participant nucleons ǫ(r) = αn part (r) + γn coll (r).
Here, we produce an ensemble of realistic initial distributions through the Glauber-based event generator GLISSANDO [13]. We take n = n coll as above and for simplicity identify the initial density of the fireball as ρ(r) ∼ n coll (r) as well, as in some well-known energy loss model calculations [14]. We do not implement a time evolution, since we only look at deviations of observables from their counterparts in smooth, averaged collisions. In other words we are only sensitive to the time evolution of R(r 1 , r 2 ). However we can argue that the longitudinal expansion of the fireball will not change the transverse correlation function R except for an overall scaling factor, and transverse expansion is building up from zero at early times, being not overly relevant for most measured jets. Smooth fireballs are created by averaging over 500,000 GLISSANDO events in the corresponding centrality bin. Fig. 2 compares a typical single event around b = 3.2 fm with the averaged event of the same centrality. The highly fragmented nature of this fireball is evident. We use GLISSANDO with the default values provided [13]. All our runs have the following choices made: binary collisions, no superimposed weights, and variable-axes quantities. To make contact between a range of impact parameters b in GLIS-SANDO and experimental centrality bins we use the tables in [15].
We use the software package PPM to calculate jet quenching results. PPM is a modular code developed by us to calculate hard probes observables. Here we run it in a mode that propagates samples of hard partons on eikonal trajectories through the background fireball with different leading particle energy loss models selected. We let PPM read in GLISSANDO output for both the distribution of hard processes and as a map for the fireball for a given event. The initial momentum distribution of quark and gluon jets used follows a leading order pQCD calculation [16,17]. As we check against pion data PPM uses the option for KKP fragmentation [18] which gives reasonable results for pions. Finally PPM computes R AA , I AA , and v 2 . For leading particle energy loss we espouse two options in PPM: (i) a simple, deterministic, LPMinspired model (sLPM) in which ∆E = c sLPM h 1 where h 1 is given by Eq. (4). The parameter c sLPM measures the relative quenching strength c sLPM = q(r)/n coll (r). (ii) the energy loss model known as the Armesto-Salgado-Wiedemann (ASW) formal- ism, which is non-deterministic. Instead, it assigns a probability density for energy loss which is given as [19] P (∆E; R, ω c ) = p 0 δ(∆E) + p(∆E; R, w c ) (13) where p 0 is the probability to have no mediuminduced gluon radiation and the continuous weight p(∆E) is the probability to radiate an energy ∆E if at least one gluon is radiated. In order to find these two quantities for each trajectory in our fireball we define ρ(r) = c ASW n coll , and PPM computes the integrals h 1 (r, ψ) and h 2 (r, ψ). The probability distributions are evaluated in the multiple soft scattering approximation (the ASW-BDMPS formalism) [19], by using the relations introduced in [14] : For the Monte Carlo sampling of the distribution we choose the non-reweighting algorithm explained in Ref. [14]. As in scenario (i) the parameter c ASW gives the quenching strength per density. We fit the energy loss parameters c sLPM and c ASW by comparing PHENIX data on neutral pion suppression R AA at top RHIC energy for three different centralities: 0 − 10%, 20 − 30% and 50 − 60% [20] to PPM calculations using averaged fireballs for three corresponding impact parameter bins.
The extracted values are c sLPM = 0.055 GeV and c ASW = 1.6 GeV. We note that ASW requires a much larger relative quenching, but we do not want to focus on a comparison of different energy loss Figure 4: R AA of neutral pions for three centrality bins computed in the ASW energy loss model compared to PHENIX data [20]. For each bin we show the result for the averaged fireball (dotted lines) for c ASW = 1.6 GeV, for the eventby-event computation (dashed lines) for the same quenching strength, and the event-by-event calculation for the refitted value of c ASW = 2.8 GeV (solid lines). models here. We simply use two models to estimate the uncertainties associated with our incomplete attempts to quanitify energy loss, and we only focus on relative changes between the smooth and the event-by-event case.
Next we run PPM over samples of individual events and then take the average of our observables, keeping c sLPM and c ASW constant. For all centralities and all values of P T we observe an increase in R AA , i.e. a consistently lower energy loss event-by-event compared to results using an averaged fireball. Fig. 3 compares both scenarios and PHENIX data for central collisions (around b = 3.2 fm) using both the sLPM and ASW energy loss. The deviations grow going from central to peripheral collisions. Now we check whether the decreased suppression can be absorbed in a redefinition of the quenching strength. Indeed, at not too large transverse momentum we find that calculations of R AA using event-by-event quenching can be fit to describe the P T -and centrality dependence of RHIC data by increasing c sLPM to 0.085 GeV and c ASW to 2.8 GeV. Fig. 4 shows the results for R AA for a central, a mid-central and a peripheral bin using ASW energy loss. For each bin three curves are compared to Figure 5: The azimuthal asymmetry v 2 of neutral pions as a function of P T for impact parameters around b = 11 fm compared with data from PHENIX [21]. We show computations in the ASW model using the average event (solid line), an event-by-event calculation (dotted line) with c ASW fitted to the R AA using the average event. We also show the event-by-event case for c ASW = 2.8 GeV which fits the R AA in the event-by-event case (dashed line). data from PHENIX: calculations with (i) the average and (ii) event-by-event fireballs using the old fit values for c ASW , and (iii) the event-by-event results using the newly adjusted parameter c ASW . sLPM energy loss leads to a similar picture. At low transverse momentum the new fits can be matched perfectly to the original curves from smooth fireballs, while at high P T differences can occur, however well within experimental error bars. We conclude that the use of smooth fireballs could underestimate the extracted energy loss coefficient by as much as 50% in the ASW model compared to an event-by-event analysis, and still by as much as 25% in the sLPM model. Suppose we do not trust GLISSANDO to capture spatial details of the initial collision correctly. We can still make the following model independent statement: There is an (additional) uncertainty of up to a factor 2 on extracted values ofq coming from the unknown event-by-event geometry of the fireball.
Let us proceed to discuss the azimuthal asymmetry v 2 . As expected from our analytic arguments the value of v 2 decreases for all centrality bins and for both sLPM and ASW energy loss if event-byevent computations are compared to the average fireball with the quenching strengths c sLPM and c ASW fixed. However, we observe that readjusting the strength to fit R AA for all centrality bins does not bring v 2 to the level observed for smooth Figure 6: The two-hadron correlation I AA of neutral pion triggers between 7 and 9 GeV/c and associated charged pions as a function of associated P T for impact parameters around b = 3.2 fm. Trigger particles are counted in a window from 7 to 9 GeV/c. We show computations in the ASW model using the average event (dotted line) with an event-by-event calculation (solid line) with c ASW fitted to the R AA using the average event. We also show the event-by-event case for c ASW = 2.8 GeV which fits the R AA in the event-by-event case (dashed line). PHENIX data for π 0 -charged hadron correlations and the same trigger window are taken from Ref. [23].
fireballs. This residual effect increases with impact parameter b. Fig. 5 shows the calculated values of v 2 for ASW energy loss for impact parameters around b = 11 fm compared to PHENIX data [21]. At a transverse momentum of 4 GeV/c the residual suppression of v 2 is about 25%. This rather deepens the puzzle of v 2 calculations which are routinely underpredicting the experimentally observed values at large P T [22]. On the other hand, an interesting possibility takes shape. Looking at R AA alone did not give us any handle on the geometry of the fireball since a simple rescaling of the energy loss parameters could absorb the effect. Looking at v 2 in addition could in principle put experimental limits on inhomogeneities in the fireball. Fig. 6 shows our results for the triggered dihadron correlation function I AA in the ASW model for impact parameters around b = 3.2 fm. We see that quenching in the average event is larger than for the event-by-event scenario, analogous to the single hadron case. J AA rises by up to 25% in the event-by-event case with fixed quenching strength, but this is almost canceled by the corresponding rise in R AA such that the modified per trigger yield I AA is almost unchanged. However, when we use the quenching parameter that fits R AA for event-by-event computations to data we observe that the refitting of R AA overcompensates the effect for di-hadron quenching in a dramatic fashion. I AA from event-by-event computations is now up to 25% smaller than for the averaged event. This overcompensation could serve as another signature for inhomogeneities. It is observed for both sLPM and ASW energy loss models. We conclude that a blend of single and dihadron measurements supplemented with v 2 measurements can in principle discriminate between different scenarios for the density correlation function R. At this point the uncertainties in I AA data are still somewhat large and quantitative estimates are not yet conclusive.

Summary
We have shown that realistic fluctuations and inhomogeneities in the fireball can have significant effects on jet quenching. We tie the deviation of single particle suppression from that in an average fireball to a path integral over the correlation function n(r 1 )ρ(r 2 ) between the fluctuations in the density of hard processes and the density of the medium. We predict that for a fixed quenching strengthq(ρ) v 2 should be diminished for a wide class of energy loss models, while the sign of the correction to R AA is less obvious and depends on details of the correlation function and the energy loss model used. We expect less suppression for event-by-event jet quenching in the limit of very strong, surface-dominated quenching.
We have verified numerically with two energy loss models that at RHIC energies single hadron suppression R AA is decreased for realistic eventby-event quenching. On the other hand v 2 is decreased as expected. The quenching strengthq as a function of the medium density ρ can be increased to describe the observed single particle suppression in event-by-event calculations. In fact, we can not distinguish, at low transverse momentum, between smooth and inhomogeneous fireballs using the P T -and centrality dependence of R AA alone if the quenching strengthq(ρ) is an adjustable parameter. The quenching strength has to be increased by up to 100% which can be interpreted as an additional uncertainty in the extraction ofq from data.
We observe that v 2 is still suppressed by up to 25%, and I AA is decreased by the same amount even after adjusting the quenching strength to fit the data on single hadron suppression. This residual signal of inhomogeneities can in principle be used for a true tomography which can measure the degree of initial fragmentation in the fireball. Of course this is only viable with di-hadron data that has significantly smaller error bars, and once theoretical uncertainties from other sources in energy loss calculations are under control.
This work was supported by CAREER Award PHY-0847538 from the U. S. National Science Foundation, RIKEN/BNL and DOE grant DE-AC02-98CH10886. E. R. thanks the Cyclotron Institute at Texas A&M for its hospitality and the National Science Foundation for support of the REU program under award PHY-0647670.