Jet Tomography of Harmonic Fluctuations in the Initial Condition of Heavy Ion Collisions

In this paper we study the jet response (particularly azimuthal anisotropy) as a hard probe of the harmonic fluctuations in the initial condition of central heavy ion collisions. By implementing the fluctuations via cumulant expansion for various harmonics quantified by $\epsilon_n$ and using the geometric model for jet energy loss, we compute the response $\chi^h_n=v_n/\epsilon_n$. Combining these results with the known hydrodynamic response of the bulk matter expansion in the literature, we show that the hard-soft azimuthal correlation arising from their respective responses to the common geometric fluctuations reveals a robust and narrow near-side peak that may provide the dominant contribution to the"hard-ridge"observed in experimental data.


I. INTRODUCTION
The structure and properties of QCD matter under extremely hot and/or dense conditions are of fundamental interest and provide unique environments for studying the strongest force of Nature. The hot deconfined QCD matter, the so-called quark-gluon plasma (QGP), was part of the history for cosmic evolution after the Big Bang and has now been created via relativistic heavy ion collisions (the "Little Bang") and explored in laboratory experiments at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). In such collisions, highly energetic jets born from initial hard scattering provide natural "tomography" of the created hot QCD matter. Jet quenching due to energy loss along the jet path through the medium encodes essential information about the dynamics of jet-medium interaction and the medium properties as well, which shall be inferrable from experimental observables such as the highp t hadron suppression and azimuthal anisotropy, as well as the di-hadron correlations (for reviews see e.g. [1]).
While the jet quenching has been experimentally established as a very robust phenomenon at RHIC and now at LHC, the microscopic mechanism of jet energy loss is not yet fully understood. A conventional observable for quantifying jet quenching is the nuclear modification factor R AA which compares the particle production in the AA collision to the naive expectation from simply scaling up single NN cross section by the binary NN collision number. A measured R AA significantly smaller than unity for the high-p t hadrons implies strong in-medium energy loss for the jets: indeed for central collisions at both RHIC and LHC we've seen R AA ≤ 0.2. The R AA provides direct information on the average opaqueness of the created hot medium, and it is customary in various jet quenching models to use R AA in the most central collisions to normalize their respective parameters for the average jetmedium interaction strength. More sensitive are the geometric features of jet quenching observables that are par-ticularly useful in discriminating different models of energy loss. These include the A-dependence (when changing colliding systems), the b-dependence (when changing the collision impact parameter or centrality class), and the φ-dependence (when changing the probe jet's azimuthal orientation with respect to the reaction plane). For any given model with its parameters fixed in the most central collisions, the above geometric dependence and the correlations among different observables then provide crucial tests of the model [2][3][4][5][6][7][8][9][10].
Let's elaborate a bit on the φ-dependence of R AA . In non-central collisions, the medium "thickness" as seen by a penetrating jet depends on the azimuthal angle φ of the jet with respect to the reaction plane, therefore leading to the reaction-plane dependence of high-p t hadron suppression i.e. R AA (φ) [2]. The dominant anisotropy in R AA (φ) (for non-central collisions) can be attributed to the second harmonic term cos(2φ − 2Ψ R.P. ) with its coefficient being the elliptic "flow" parameter for highp t hadrons, V hard 2 which is a non-collective component though [11]. Despite the success of many models in describing R AA and its centrality dependence, most models significantly under-predicted the V hard 2 and failed the test by geometric data [3] [9]. The lack of a simultaneous description for R AA and V hard 2 in a single model was not resolved till a new insight suggested in [4]. Motivated by the "magnetic scenario" for sQGP [12], the authors of [4] pointed out that the energy loss of a jet may not simply scale with the local medium density as most models have assumed, but actually has nontrivial dependence on matter density (or temperature). It was particularly shown that the geometric data R AA and V hard 2 versus centrality can be successfully described together by including a jet quenching component with strong enhancement in the near-T c matter by a factor 3 ∼ 5 compared with higher-T QGP. Such an enhancement of jetmedium interaction may originate from non-perturbative structures created by the (color-)electric jet passing a plasma of (color-)magnetic monopoles that dominate the near-T c matter [12] [13]. A natural prediction of this scenario is that the effective jet-medium interaction would be rapidly reduced when going from RHIC to LHC energies, or the hot matter created at LHC will be more "transparent" (apart from the trivial density factor) to a penetrating jet at LHC as compared with that at RHIC. Interestingly, recent quantitative computations and comparison with LHC data at √ s = 2.76TeV [14][15] [8] indeed suggest that 1) fixing the jet-medium interaction at RHIC and applying it directly to LHC would lead to "over-quenching" as compared with data; 2) reducing the jet-medium interaction by about a factor of 2 would allow a good description of the LHC data -such reduction is remarkably rapid provided only a 30% increase in temperature from RHIC ( √ s = 0.2TeV) to current LHC ( √ s = 2.76TeV)! The future LHC heavy ion data at √ s = 5.5TeV will be essential in establishing how much and how fast the jet-medium interaction decreases with temperature.
In this paper we explore a new and interesting geometric aspect of jet quenching: the hard probe of the geometric fluctuations (in terms of various harmonics in azimuthal angle) in the initial condition of heavy ion collisions. It has recently been shown in measurements [16,[18][19][20]36] and demonstrated in various modelings [21][22][23][24][25][26][27][28][29][30][31][32][33] that there are very strong fluctuations in the initial matter profile from event to event. Such fluctuations contain various higher order harmonics in azimuthal angle ∼ cos(nφ) rather than the naive expectation of 2ndharmonic dominance in averaged geometry. In the collective expansion of the bulk matter, these fluctuations lead to observed harmonic flows up to about n = 6 and also explain the soft di-hadron azimuthal correlations (the "soft-ridge"). In particular, it was shown in [22] that there is primarily a linear response of the low-p t hadrons' harmonic flows to the initial harmonic fluctuations. A natural question to ask is, therefore, how a penetrating jet responds to the strongly fluctuating matter density from event to event. There have been some discussions in the literature [5][7] [34] and also data from both RHIC [35] [36] and LHC [37] [38], but a systematic and clear picture has not been achieved. Intuitively one would expect the event-by-event azimuthal distribution of highp t hadrons shall also reflect such harmonic fluctuations in the initial condition and it is our purpose to systematically quantify such response. In addition, both the quenching of hard jet and the expansion of soft matter are commonly correlated to the same underlying particular matter profile in a given event and such hard-soft correlation will survive the event average. We will examine how such correlation contributes to the di-hadron azimuthal correlation with a hard trigger and a softer associate hadron, in particular a possible explanation of the "hard-ridge" [39] for which the origin has so far not been fully understood [40][41] [42]. To highlight the role of pure fluctuations in leading to azimuthal anisotropy, we limit our study in the present paper only for the most central collisions, and to be specific we do calculations for RHIC √ s = 200GeV collisions. Unlike in non-central collisions with strong anisotropy already in the average geometry and dominated by the 2nd harmonics, in the perfectly central collisions (b = 0) the average background geometry is isotropic and the anisotropy from various harmonic fluctuations will be best manifested.
The paper is organized as follows. In section II we will describe our model setup for the present study. The simulation results for jet response to the harmonic fluctuations in central collisions, characterized by v n /ǫ n for n = 1, 2, ..., 6, will be presented in section III. Such results will be used in section IV to compute the hardsoft di-hadron correlations due to harmonic fluctuations where we will find a robust and narrow near-side peak in the azimuthal angle ∆φ dependence as a possible explanation of the so-called "hard-ridge". The summary and discussions will be given in section V.

II. MODEL SETUP
In this section we introduce our model setup. First we will show how we implement various harmonic fluctuations on top of the isotropic fireball in a perfectly central collision by using a modified cumulant expansion based on that of [26]. In particular we will discuss the fluctuations in both the participant density and the binary collision density profiles and their mutual relations. Second we will briefly discuss the geometric model we use for calculating the jet energy loss and its azimuthal angle dependence for a given matter distribution.

A. Parametrization of Harmonic Fluctuations in Central Collisions
In each central collision event, the participant density (ρ p ) and the binary collision density (ρ c ) are related to the thickness functions of the two colliding nuclei (T a,b ): Here σ(= 42 mb) is the N-N inelastic scattering cross section at √ s = 200GeV, A is the number of nucleons in each nucleus, and It should be emphasized that in general one shall use T a ( r − b/2) and T b ( r + b/2), and the above is only to be used for the perfectly central collisions with b = 0. Suppose event by event, T a,b fluctuate around an isotropic background T 0 , i.e. T a,b = T 0 a,b + δT a,b . As the result, by using the linear approximation, the fluctuations of ρ c and ρ p are directly related, as shown by T 0 is obtained by the (optical) Glauber model. To parameterize the fluctuations, we focus on δρ p and later use the above equations to compute the corresponding δρ c . We will use a modified cumulant expansion based on the methods in [24,26]. Generally, we can make a Fourier analysis for any profile The difficulty is to parameterize the function form for the coefficients ρ p,n (r). The small-k expansion method [24,26] starts with the transformation (for 2-D transverse plane) ρ p ( k) ≡ d rρ p ( r) exp (i k r), and examines instead the Fourier series ρ p,0 (k) + n 1 2ρ p,n (k) cos n(φ k − Ψ k n ) . One can then make a small-k expansion of the coefficients ρ p,n (k) like e.g.
is irrelevant to the azimuthal anisotropy and will not be addressed here. It is easy to see that [24] for the n − th harmonic (n 2), the leading term in ρ p,n (k)'s k expansion is ρ p,n,n (= 1/2 n < r n cos [n(φ − Ψ n )] >, < ·· > means averaging over normalized ρ p ( r)). The case for n = 1 deserves special treatment: it starts with ρ p,3,1 (= 3/8 < r 3 cos (φ − Ψ 1 ) >) as one can set ρ p,1,1 = 0 by choosing the "right origin" of coordinates. Practically one then has to truncate the series, e.g. keeping only leading terms (assuming the relevant k is small). Due to the truncation, upon transforming back to the r-space one needs to regulate the large r part. Different from previous approaches, we will simply use a Gaussian factor with a length parameter Σ that roughly reflects the scale of fluctuations. In the end, we have the following parameterizations for various orders of harmonics: Here ǫ 1 ≡ − < r 3 cos (φ − Ψ 1 ) > / < r 3 >, ǫ n 2 ≡ − < r n cos n(φ − Ψ n ) > / < r n >, < r n > 1/n = 3.88, 4.19, 4.44, 4.63, 4.80, 4.95 fm for n = 1, 2, 3, 4, 5, 6. The root-mean-radius-square √ < r 2 > = 4.19 fm is defined as σ and will be used as a "ruler" for Σ. Therefore by specifying the values of ǫ n (as well as Ψ n ) for an event, one then fixes the ρ p,n (r) above and obtains a particular density profile in Eq.(3) with harmonic fluctuations.
One technical complication is that for such parametrization at large radius, the total density can become negative.
Hence certain regularization scheme has to be implemented as discussed in details by [26]. We simply set the density beyond that critical radius as zero. For the condition of our interest, i.e. ǫ ≈ 0.1, this critical radius is around 7 fm and our regularization is physical. The regularization will require a re-calibration of the eccentricities: all the actual ǫ n values need to be evaluated directly from the generated and regularized profile which would generally differ from the input ǫ n parameters in Eq.(4), see e.g. [26] [23].

B. The Jet Energy Loss
In this study, we use the geometric model approach for computing the jet energy loss. Such models reflect the generic geometric features (e.g. the path-length dependence) that are most crucial for describing geometric data [3][4][5][6][7][8]. We assume that the final energy E f of a jet with initial energy E i after traveling an in-medium path P can be parameterized as E f = E i × f P with the suppression factor f P given by: In the above the s(l) is the entropy density of local matter at a given point on the jet path, while the κ(s) is the local jet quenching strength which as a property of matter should in principle depend on the local density s(l).
There can be different choices of the parameter m for path-length dependence (e.g. LPM-motivated quadratic or AdS/CFT-motivated cubic) and of the jet-medium interaction κ(s). In the present study, we use the near-T c enhancement model as in [4][8], assuming m = 1 and introducing a strong jet quenching component in the vicinity of T c (with density s c and span of s w ) via with ξ = 6, s c = 7/f m 3 , and s w = 2/f m 3 . (see [4] for details.) As aforementioned, the parameter κ will be fixed by R AA ≈ 0.18 in the 0 − 5% collisions at RHIC √ s = 200GeV. It shall be emphasized that the jet path P is determined by the initial jet spot coordinates on the transverse plane as well as the the azimuthal angle φ for the transverse orientation of its propagation. After averaging over all jet paths (including all possible start points properly weighed by binary collision density and all equally possible orientations) one may then obtain the R AA : where the exponent n comes from measured reference p-p spectrum (see e.g. [9] for a detailed account). The value of n depends on collision energy and throughout this paper we focus on the RHIC √ s = 200GeV collisions with n ≈ 8.1. Note that in this model as in many other geometric models, assuming the fractional energy loss will lead to the R AA independent of p t which may be justified by the approximate "flatness" seen in the data at RHIC energy. The results from such modeling apply only to the high-p t region, e.g. p t > 6 GeV at RHIC. The azimuthal angle φ-dependence can be studied by averaging over jet paths with a particular azimuthal orientation, i.e. R AA (φ) =< (f P ) n−2 > P (φ) from which various harmonic components are derivable by Fourier decomposition. In this study for each given geometry, we compute R AA (φ) as a function of azimuthal angle φ by integrating over all initial jet spots on the transverse plane weighed by the binary collision density at each spot for a given value of φ.

III. JET RESPONSE TO HARMONIC FLUCTUATIONS
Here we report our results for the jet response to harmonic fluctuations: see Fig. 1. The results are obtained with Σ = σ = 4.19 fm, which indicates the scale of the fluctuation is roughly the same as that of the isotropic background. In our approach, the fluctuations ǫ n are assumed to be small so that the jet response to densityfluctuations is approximately linear. We will first demonstrate that indeed for n − th harmonics alone with amplitude ǫ n , the jet response to a good approximation depends on ǫ n linearly (hence we can define χ h n = v n /ǫ n ). The solid curves in each panel clearly show that for n = 1, 2, 3, the relations between v n and ǫ n are linear, while for n = 4, 5, 6 only very mild nonlinearity starts to develop near ǫ n = 0 and the response stays quite linear for the physically most relevant region 0.05 < ǫ n < 0.1. For each harmonic fluctuation, we have checked that the initial axis angle Ψ n of fluctuation coincides with the final axis angle Ψ J n of the jet distribution anisotropy. One important issue is to understand the different influences on the jet response due to the fluctuation in the participant density (which mostly concerns the bulk matter that quenches the jets) and that in the binary collision density (which concerns the profile of initial jet spots). The difference could be tricky and counter-intuitive (see e.g. some discussions in [34]). Here we clearly demonstrate the difference by showing three different responses in Fig. 1: the case with only ρ p fluctuation ("Part", green dashed curves), the case with only ρ c fluctuation ("Coll", blue dotted curves), and the physical case with both ("Tot", red solid curves). We've checked that the total response agrees well with the sum of the two individual contributions which provides additional evidence for the linear nature of the response. As one can see, the contributions to the jet response due to fluctuation in participant density and that in collision density differ significantly in the absolute magnitude and can even have opposite signs (in the case of n = 1, 2, 6) thus canceling each other to certain extent. It is therefore clear that for event-by-event studies of jet response, both fluctuations have to be fully taken into account.
Finally in the Tab. I, we show the jet response coefficients χ h n extracted from the physically most relevant region 0.05 < ǫ n < 0.1 where the linear dependence is a very well approximation for all harmonics. The results as a response spectrum χ h n vs. n are also plotted in Fig. 2. The response spectrum shows a typical decrease from low to high harmonics, and somewhat surprisingly (as compared with the soft response) there is a very strong response in the first harmonics. We also see the responses are relatively insensitivity to the parameter Σ and therefore the quantified response spectrum shape in n is robust. While it would be interesting to directly compare with data, currently the RHIC data [35] [36] for higher harmonics are only available for the low to intermediate p t region so such comparison with our results would not be appropriate. We point out though the previous comparison between the results from the same model with v 2 data at high p t was rather successful, see e.g. [4] [8].  We end this section by discussing certain check we've done. One issue is about possible "mixing" or "interfering" effect among different harmonic fluctuations when they all co-exist with varied respective main-axis for each. This should be clarified as in reality each collision event would automatically come with all these harmonic fluctuations. We've done the following test: we include on top of the isotropic background all the n = 1, 2, 3, 4, 5, 6 harmonic fluctuations simultaneously with randomly assigned ǫ n (in the linear regime though) and randomly chosen axis Ψ n (different for each), and then calculate the jet quenching result R AA (φ). We've found that 1) the main-axis for each harmonics Ψ J n as determined by maximizing the corresponding cos(φ − Ψ J n ) component in final state R AA (φ) agrees with the initial Ψ n from fluctuation with less than 1% difference; 2) the response v n /ǫ n also agrees with the values we extracted above by treating each harmonics separately. Another issue that we've checked is how such response would change with different jet energy loss models. We've also done the calculation for two other models: the model with pathlength-square dependence and constant jet-medium interaction, i.e. κ(s) = κ and m = 1 in Eq.(6); and the model with path-length-cubic dependence and constant jet-medium interaction, i.e. κ(s) = κ and m = 2 in Eq. (6). The response coefficients χ h n from different models are not drastically different, show similar decreasing trend with increasing n, and may become distinguishable when realistic and accurate data comparison can be done.

IV. THE "HARD-RIDGE" IN DI-HADRON CORRELATIONS
In this section, we study the correlation between the hard and soft sector due to their respective correlations to the common initial condition with fluctuations. It would be interesting to see what major features of the hard-soft correlation structures can arise from the jet response to harmonic fluctuations in the initial condition.
Suppose for an arbitrary central collision event, we have the initial azimuthal distribution of matter density fluctuations characterized by a series of harmonics with varied (and uncorrelated) magnitude ǫ n and axis orientation Ψ n . According to our results on the jet response, we would expect the final state high-p t hadron distribution  to be of the form The bulk matter, on the other hand, will generate the harmonic flows during hydrodynamic expansion and lead to the final state low p t hadron distribution of the form Note that for a given event both the soft and the hard responses to each order of the harmonic fluctuations are commonly aligned to the same corresponding angle Ψ n from the initial condition -such a correlation between the hard and soft responses will survive the average over many events. We can therefore define a hard-soft correlation function (11) where the < > means averaging over events. Since the initial harmonic fluctuations have their orientation angles Ψ n vary randomly from event to event and uncorrelated among each other (at least so in central collisions), we obtain from Eq.(9)(10) the following In the linear response approximation we can use v h,s n = χ h,s ǫ n to further get We therefore see that the common correlations of the hard and soft responses with various harmonic fluctuations indeed lead to the hard-soft correlation that survives the event average. It it tempting to quantitatively examine the features of such hard-soft azimuthal correlation. To do that, we compute C[∆φ] by including the n = 1, 2, 3, 4, 5 harmonics in Eq. (13).
For the initial fluctuations we use results from Monte Carlo simulations of initial conditions for most central collisions, which are ǫ 2 n=1,2,3,4,5 ≈ 0.037, 0.068, 0.076, 0.084, 0.091 respectively (see e.g. [22] [24][26] [27]). For the soft response coefficients we take the results from hydrodynamic modeling, which are χ s n=1,2,3,4,5 ≈ 0.15, 0.26, 0.21, 0.14, 0.086 respectively (see e.g. [22] [24][26] [27]). The hard response coefficients are taken from our results: since such coefficients depend on the parameter Σ we will show results for both values of Σ as studied in the previous section. The so-obtained correlation function C[∆φ] is plotted in Fig.3: the blue solid curve is obtained using χ h n results from Σ = σ case while the red dashed curve is from Σ = 1.25σ case. A few remarks are in order: 1) there is a robust and narrow near-side peak around ∆φ = 0. Actually for both cases we've fitted very well the peaks by a Gaussian function with width about 0.8 radian which is very close to the experimental data. This main feature, we believe, could provide a natural expla-nation of the "hard-ridge" structure as arising from superposition of multiple harmonics (dominantly from the first three harmonics). 2) the away-side structure around ∆φ = π shows an interesting double-hump shape, with a "dip" at π and two "shoulders" around 2π/3 and 4π/3. This structure is mainly due to an interplay between the v 1 response and the v 3 response and also quite robust despite the two different choices of the parameter Σ. The precise shape of the away-side, though, is very sensitive to the ratio between the responses to the first and third harmonic fluctuations. Qualitatively these features agree well with the experimentally observed di-hadron correlations. While certainly interesting to do, a quantitative comparison between the obtained correlation and the experimental data is not ready yet, due to a number of approximations used in the modeling and also the complication in experimental methods (e.g. observable definition, trigger and associate p t selection, the ZYAM procedure for background subtraction, etc). Some of these issues will be discussed further at the end.
One might also think about possible measurements of hard-hard correlations, i.e. di-hadron correlations with both hadrons' p t to be high enough that they are most likely from jets. Such measurements might become feasible at LHC energies when enough events with more than one pair of jets could be collected. Along similar consideration as before we may expect a component in the hardhard azimuthal correlation due to geometric fluctuations of the form ∼ n=1,2,3,... 2 (χ h ) 2 < (ǫ n ) 2 > cos(n∆φ). A remnant of a near-side ridge might be visible, while the usual di-jet back-to-back correlation will be dominant.

V. SUMMARY
In summary, we've studied the jet response (particularly azimuthal anisotropy) as a hard probe of the harmonic fluctuations in the initial condition of central heavy ion collisions. By implementing the fluctuations via cumulant expansion for various harmonics quantified by ǫ n and using the geometric model for jet energy loss, we've computed the response χ h n = v n /ǫ n . Combining these results with known results for hydrodynamic response of the bulk matter expansion in the literature, we've shown that the hard-soft azimuthal correlation arising from their respective responses to the common geometric fluctuations reveals a robust and narrow nearside peak that may provide the dominant contribution to the "hard-ridge" observed in experimental data. This paper is intended to demonstrate and emphasize the main idea of hard probe for harmonic fluctuations in the initial condition, and we end by discussing a number of highly interesting aspects for further developments as well as some issues that are not fully addressed in the present paper and will be further investigated. The approach of quantifying fluctuations by cumulant expansion with one harmonic fluctuation at a time has the advantage of clearly demonstrating the response of jet energy loss to each order of harmonics, but for realistic modeling and comparison with data one needs to use Monte-Carlo generated initial matter profiles born with various harmonics and more irregularity -this has now been implemented and the results are to be reported in a forthcoming publication. Using the Monte-Carlo initial conditions will also allow extending our study to non-central collisions as well as eliminating ambiguity due to choice of parameters. In addition the "reaction plane" dependence of the hard-soft correlation will be exploited with the Monte-Carlo generated fluctuations. At present the various jet energy loss models suffer from uncertainty for the energy loss in the pre-equilibrium matter [5][7] [43], and it would be interesting to find ways of using geometric features of jet quenching to tightly constrain such uncertainty. While explanation of the near-side peak structure in the hard-soft di-hadron azimuthal correlation by the harmonic fluctuations appears to be robust, the full explanation of the observed away-side structure is much tricker due to known contributions from various other sources, like e.g. transverse expansion dynamics [44] and background effects like transverse momentum conservation and cluster correlations [45], which all require further scrutiny. All that said, we emphasize again that the hard probe of geometric fluctuations provides a new and useful tool for exploring both the initial condition and the jet energy loss mechanism in heavy ion collisions.