High Production with Respect to the Reaction Plane
in Collisions at GeV
Abstract
Measurements of the azimuthal anisotropy of high neutral pion () production in Au+Au collisions at GeV by the PHENIX experiment are presented. The data included in this paper were collected during the 2004 RHIC running period and represent approximately an order of magnitude increase in the number of analyzed events relative to previously published results. Azimuthal angle distributions of s detected in the PHENIX electromagnetic calorimeters are measured relative to the reaction plane determined eventbyevent using the forward and backward beambeam counters. Amplitudes of the second Fourier component () of the angular distributions are presented as a function of transverse momentum () for different bins in collision centrality. Measured reaction plane dependent yields are used to determine the azimuthal dependence of the suppression as a function of , . A jetquenching motivated geometric analysis is presented that attempts to simultaneously describe the centrality dependence and reaction plane angle dependence of the suppression in terms of the path lengths of hypothetical parent partons in the medium. This set of results allows for a detailed examination of the influence of geometry in the collision region, and of the interplay between collective flow and jetquenching effects along the azimuthal axis.
pacs:
21.65.Qr,25.75.q,25.75.DwPHENIX Collaboration
I Introduction
Over the past few years, experiments at the Relativistic Heavy Ion Collider (RHIC) have established that a dense partonic medium is formed in collisions at =200 GeV Adcox et al. (2005); Adams et al. (2005a); Back et al. (2005a); Arsene et al. (2005). This medium thermalizes very quickly Adcox et al. (2005); Mrowczynski (1993); Arnold et al. (2005); Rebhan et al. (2005); Dumitru and Nara (2005); Schenke et al. (2006); Scherer et al. (2008); Xu et al. (), is extremely opaque to the passage of high particles Adcox et al. (2002); Adler et al. (2002), and the strong coupling of matter in the medium produces a system for which the ratio of shear viscosity to entropy () approaches zero Adler et al. (2003a); Adams et al. (2004); Romatschke and Romatschke (2007); Dusling and Teaney (2008); Song and Heinz (2008). Much of the current focus is on the extraction of key transport and thermodynamic characteristics of the matter produced in these collisions. Measurements of high parton propagation in the medium as well as mediuminduced modification of the fragmentation parton spectrum and its products provide a critical tool for probing medium properties.
One of the most striking early results from RHIC was the observation of strongly suppressed production of high particles in central Au+Au events compared to appropriately scaled collisions Adcox et al. (2002); Adler et al. (2002). High partons are formed from hard scattering between the initial colliding partons, and these partons fragment into two or more jets of hadrons. When propagating through a dense volume of deconfined matter, these high partons are expected to scatter from color charges in the medium, losing energy through a combination of gluon bremsstrahlung radiation and collisional energy transfer to partons in the medium. These radiated gluons eventually fragment into hadrons at lower , resulting in a depletion of the observed yields of hadrons at higher .
A useful way to quantify the suppression of high hadrons is the nuclear modification factor () where the cross section is scaled by the thickness function of the two Au nuclei
PHENIX has measured a close to unity in both peripheral collisions and Ạu collisions Adare et al. (2008a); Adler et al. (2003b), consistent with the expectation that these collisions would not produce an extended, dense medium. As the collisions become more central, decreases to about 0.2, indicating a stronger parton energy loss. Furthermore, the measured is nearly constant as a function of , for up to the highest currently accessibly , GeVc Adare et al. (2008a).
These data can be well reproduced by models that calculate the energy lost by the hard scattered partons as they traverse the dense medium. The amount of energyloss depends on the density of the medium Baier (2003), so measurements of high hadron suppression provide constraints on the transport coefficient , a measure of mean transverse momentum squared transferred by the medium to a highenergy parton. However, multiple models with different physical assumptions can reproduce the measured Majumder (2007); Vitev (2008). The different models vary widely in how they include the crucial interference terms between multiplescattering centers as well as the interplay between inelastic, elastic and flavorchanging processes during the parton’s passage.
To discriminate between these models we need to increase our experimental control of the path length, since the amount of energy lost by a high parton strongly increases with the distance traveled through the medium. A quadratic dependence on the path length is predicted for a static medium if the dominant energyloss mechanism is the bremsstrahlung radiation of gluons surviving the destructive interference caused by multiple scattering Majumder (2007); Vitev (2008). For an expanding plasma the quadratic increase should be moderated to a linear dependence Gyulassy et al. (2002).
The centrality dependence of offers a probe of the pathlength dependence of partonic energy loss. However, we can better test the pathlength dependence by studying the azimuthal variation of the high suppression at a fixed centrality. Since the collision zone has a nearly elliptical shape in the transverse plane due to the noncentral overlap of the colliding Au nuclei, partons that travel along the short axis of the nuclear overlap region lose less energy and should therefore be less suppressed. The key observable is then the twodimensional modification factor , where is the angle of emission with respect to the event plane. The azimuthal dependence of the spectra can be also parameterized by a Fourier expansion, where up to second order , with being called elliptic flow coefficient. While both quantities characterize azimuthal asymmetries, historically and conceptually they have different roots. The notion of elliptic flow is primarily tied to lower phenomena (“soft physics”), the domain where particle production is proportional to the number of participating nucleons (), and positive arises from the boost to the mean in the direction where the pressure gradient is highest (along the reaction plane). Conversely, and are commonly used to describe high behavior (hard scattering, which scales with the number of binary collisions ). When deviates from unity at high , it becomes a valuable probe of the loss of energy/momentum in a particular direction. However, there is no clear separation between soft and hard regions, and both and are welldefined in the entire momentum range, so in this sense is sensitive to differential energy loss at high .
PHENIX has measured high for particles from Au+Au collisions Adler et al. (2007). The energyloss models that reproduce diverge in their predictions of the azimuthal anisotropy at high . They generally underpredict the observed azimuthal variation of , or equivalently, are unable to describe the dependence of over the full range of where one would naively expect them to be applicable Hirano and Nara (2004); Renk et al. (2007); Majumder et al. (2007). These models include the hydrodynamical evolution of the medium, and therefore the high probe loses energy in a medium that is becoming spatially isotropic with time. Several early papers noticed that the measured values were larger than what one would expect from a completely opaque almondshape collision zone Shuryak (2002); Drees et al. (2005). Other early energyloss calculations came close to reproducing the measured Dainese et al. (2005); Cole (2006), but in these the plasma expansion was not taken into account, which resulted in unrealistically strong azimuthal anisotropy. Another calculation Zhang et al. (2007) has reproduced , but in this model the Au nuclei were parameterized as hardspheres instead of using a more realistic WoodsSaxon density profile, and this mechanism artificially increases the azimuthal dependence of the energy density.
One potential resolution of the problem with energy loss calculations not reproducing the measured azimuthal dependence of yields is a recent calculation that allowed the high parton to resonantly scatter with the medium Liao and Shuryak (2008) (and references therein), increasing the energy lost by a parton at plasma densities that correspond to temperatures near the critical temperature. This produces a sharper dependence of the energyloss on the spatial variation of the medium’s energy density and hence the model is able to simultaneously reproduce both and . A critical check will be to examine whether the same parameters work for the full range of collision centralities.
In order to discriminate among all the models that attempt to reproduce , the experimental challenge is to extend the range and increase the precision of observations which can be used to test different energyloss models. In this paper we extend the range of published data on Adler et al. (2007) by a) reaching higher , and thereby moving to a region that is completely dominated by the fragmentation of hard partons and reducing the possible contribution of particles from recombination Fries et al. (2003), b) using finer bins in centrality, thus achieving less averaging of the path length, and c) reducing the statistical and systematic uncertainties to further constrain models. We present in this article measurements using data collected during the 2004 RHIC running period. These data represent a highstatistics sample of collisions (approximately 50 times that of the 2002 RHIC running period) and therefore extend our ability to measure and to much higher with better precision.
Ii Experimental Details
The data presented in this paper were taken by the PHENIX experiment Adler et al. (2003c) in 2004 (RHIC Run4), and represent the analysis of 821M minimum bias Au+Au collisions at GeV. The detectors involved in this analysis are the beambeam counters Allen et al. (2003) (BBC; triggering, centrality and reaction plane determination), the zerodegree calorimeter Adler et al. (2001) (ZDC, centrality determination) and the electromagnetic calorimeter Aphecetche et al. (2003) (EMCal, measurement).
The BBCs are two groups of 64 hexagonal quartz Čerenkov radiator counters with photomultiplier readout surrounding the beampipe 144 cm up and downstream (“North” and “South”) from the center of the nominal collision diamond, covering the pseudorapidity range and the full azimuth. Coincidence of signals in at least two photomultiplier tubes in both BBCs served as a minimum bias trigger and according to simulations it captured 92% of all inelastic collisions. The size of the total signal in the BBCs increases monotonically with collision centrality at this . The collision vertex was calculated from the difference between the fastest timing signals in the North and South BBCs, respectively, with cm resolution. Only events with cm were analyzed.
The ZDCs are small tungsten/scintillator hadron calorimeters with quartz fiber lightguides and photomultiplier readout, located between the beampipes at 18 m North and South from the collision point. They measure noninteracting “spectator” neutrons in a cone of about 2 mrad, and their signal is doublevalued as a function of centrality (it is low in very central and very peripheral collisions but large at midcentrality). The correlation of ZDC vs BBC signals resolves this ambiguity and allows a precise measurement of the true centrality for all but the most peripheral collisions.
The reaction plane (spanned by the beam direction and the impact vector of the colliding nuclei) is determined eventbyevent from the azimuthal charge distribution in the BBCs, after taking into account small nonuniformities (in the response of individual radiators, PMTs, electronics, etc.), using the assumption that over a large number of events the distribution of perevent reaction planes should be uniform. Due to the large rapidity gap between the central arm () where s are measured, and the BBCs where the reaction plane is established, we assume that the reaction plane is unbiased and free from autocorrelations. However, the relatively coarse granularity of BBCs affects the resolution. Note that in this analysis precise knowledge of the reaction plane resolution, which depends strongly on centrality, is crucial. This will be discussed in detail in the next Section.
Neutral pions are measured by reconstructing their decay photons () in the EMCal. The EMCal consists of 8 sectors at midrapidity (), covering a total of in azimuth. Six sectors are lead/scintillator (PbSc) sampling calorimeter with photomultiplier readout and cm granularity, two sectors are lead/glass (PbGl) Čerenkov counters with cm granularity and photomultiplier readout. The two detectors are 18 and 16 radiation lengths deep, respectively, both ensuring essentially full containment of electromagnetic showers in the relevant energy range. The in situ energy resolution is well reproduced by simulation both in PbSc and PbGl: the peak positions and the widths both agree with the data to better than MeV over the entire momentum range. Therefore, the error on the energy (and momentum) scale is less than 1%. Timing resolution is ps and ps for the PbSc and PbGl, respectively, allowing the rejection of neutrons and antineutrons up to a few transverse momentum, which would otherwise be a major source of neutral showers up to a few GeV energy. At sufficiently high transverse momenta, decay photons from a nearly symmetric () decay may produce showers in the calorimeter that start to merge into one reconstructed cluster. In the PbSc this effect is first visible around , at the upper end of the region considered in this paper. Due to its higher granularity and smaller Molièreradius the PbGl is immune to this “merging” problem up to 15 . The hadronic response, timing properties and other sources of systematic errors are very different for the two calorimeter types. Therefore, when extracting the integrated , which serves as absolute normalization, the PbSc and PbGl were analyzed separately and the results combined to decrease the total systematic uncertainty.
Iii Data Analysis
iii.1 Centrality
As mentioned, the minimumbias trigger in the Run4 PHENIX configuration is supplied by the BBCs, and the correlation of the charge deposited in the BBCs with energy deposited in the ZDCs provided a determination of the centrality of the collision. The elliptic flow measurement presented in this paper is measured in seven bins of the centrality range 092%, with lowest corresponding to the most central: 05%, 510%, 1020%, 2030%, 3040%, 4050%, and 5060%. In addition, the combined ranges 020%, 2040%, 4060%, and minimum bias bins are included. For the yields with respect to the reaction plane, the centralities presented are 010%, 1020%, 2030%, 3040%, 4050%, 5060%. Finally, the versus nuclear path length result excludes the most central bin due to its smaller intrinsic ellipticity (the average path length is insensitive to ).
iii.2 Reaction plane determination
The technique used to determine the reaction plane on an eventbyevent basis is the same used in previous PHENIX analyses Adler et al. (2003a, 2006); Afanasiev et al. (2007). The quartz radiators of each counter are arranged in approximately concentric circles around the beam axis. The light collected in the photomultiplier tubes (PMTs) allows for an estimate of the number of charged particles passing through the detector.
The number of charged particles at a given PMT position, , is weighted in a manner to reduce the bias of the inner rings and used to measure the orientation of the reaction plane from the formula
(1) 
where is the nominal azimuth of the radiator. Subtraction of the average centroid removes biases due to various detector effects. A final flattening technique is used to remove the residual nonuniformities in the distribution of angles.
To estimate the resolution of the reaction plane measurement, we use the subevent technique Ollitrault (). The approach consists of dividing the event up into two subevents roughly equal in size. The two individual BBCs provide a natural subevent division, so we analyze the distribution of eventbyevent differences between the reaction plane angles measured in the north and south counters, . In the presence of pure flow, this distribution takes the form Ollitrault ():
(2) 
where and the functions and are the modified Bessel functions of the first kind and modified Struve functions, respectively. The parameter describes the dispersion of the flow vector and thus determines the correction required for the reaction plane resolution. Since represents the wholeevent difference distribution and we are dealing with subevents with roughly half the multiplicity of the event, we replace in Eq. 2 and fit this function to the measured distribution to extract . The resulting value is then used to evaluate the resolution of the eventplane of order Ollitrault ():
(3) 
where the true reaction plane orientation is denoted by and the observed orientation by . Figure 1 shows the resolution correction obtained using the abovedescribed procedure as a function of centrality. Both 5% and 10% wide bins are shown for comparison.
Eq. 2 is derived under the assumption that the azimuthal distributions are free of nonflow effects. Due to the large rapidity gap between the BBCs and the central arm region, it is expected that particles observed in the BBCs have no correlation with those measured in the central arm detectors. pythia Sjostrand et al. (2006) studies have been used to confirm that jets observed in the central arm have negligible effect on the reaction plane measurement from the BBCsAdare et al. (2008b).
iii.3 Neutral pion measurement
Measurement of neutral pions has played a critical role in the study of high phenomena at RHIC, and especially by PHENIX Adcox et al. (2002); Adler et al. (2007); Adare et al. (2008a). The twoparticle decay channel provides a clean signal of identified hadrons out to the highest regions.
EMCal showers are found by clustering contiguous towers with energy above a threshold energy (10 MeV) and requiring at least 50 MeV in the tower with highest energy deposit. The impact position is calculated from the positions of the participating towers weighted by the logarithm of deposited energy. The energy of the cluster is corrected for nonperpendicular incidence – the angle being derived by assuming a straight path between the actual vertex and the calculated impact point – as well as nonlinearities Adler et al. (2007). In highmultiplicity events such as central collisions, there is an increasing probability for clusters to overlap (one tower accumulates energy from more than one particle), which can distort an energy measurement from a simple sum over contiguous towers. To mitigate this effect, the EMCal clustering algorithm also provides a quantity called ecore, which is determined by extrapolating the “core” energy represented by the central four or five towers in the cluster, assuming an electromagnetic shower profile. The energy and impact angle dependent shower profile is a model developed from and checked against beam test data. In this way, ecore provides a more realistic measurement of the shower energy, less prone to contributions from accidental overlaps (particles hitting close enough they deposit energy in the same towers) than a simple energy sum of participating towers would be. We use this ecore for the energy of reconstructed clusters in this analysis.
The invariant mass of a photon pair as measured in the EMCal is calculated from the energy of the clusters and their measured position:
(4) 
where is the opening angle between the two photons and is equal to the mass for photons from the decay of the same . Since the photons from the are not tagged, such pairs have to be formed from each photon pair in the event where the pair momentum falls in a particular bin, and some of these pairs might accidentally reproduce the mass as well (combinatorial background), particularly at lower and higher centralities (multiplicities). Since s cannot be uniquely identified, raw yields are extracted statistically, by subtracting the combinatorial background from the invariant mass distribution.
A wellknown technique to reduce the combinatorial background is to place a cut on the energy asymmetry of the pair, as defined by:
(5) 
Because the angular distribution of the pairs in the rest frame of the is uniform, the asymmetry distribution should be flat. However, due to the steeply falling photon spectrum, fake (noncorrelated) pairs which still give the proper mass are strongly peaked towards . A pair of clusters in the EMCal is considered a neutral pion candidate only if the pair’s asymmetry is less than 0.8. In addition, the two photons are required to be separated by at least 8 cm for the combination to be considered as a candidate.
There remains a nontrivial background contribution which passes these cuts: pairs of photons from different s, or, more generally, from pairs of uncorrelated clusters which pass the photon identification cuts and accidentally give an invariant mass near the true mass. This remaining combinatorial background is estimated and subtracted using the event mixing method. The procedure involves forming pairs from different events, which will by definition be uncorrelated. Each photon candidate is combined with all the photon candidates in previous events stored in memory. In order to replicate the background from uncorrelated pairs within the same event as closely as possible in the mixed events, mixing is performed within bins of centrality, vertex position, and reaction plane orientation. Since all events analyzed are minimumbias, no special steps are needed to avoid the distortions of the mixedevent background by the trigger requirement. All cuts applied to the combinations of same event pairs are also applied to mixedevent pairs. The number of events buffered determines the statistics of the eventmixed distributions, chosen as a tradeoff between desired statistical accuracy and computational resources. The data presented in this article are mixed with five previous events (in each centrality, vertex, and reaction plane bin).
For a given bin, the mixedevent mass distribution is normalized to the sameevent distribution in a region away from the mass peak. The normalization region is 0.25–0.45 GeV/ for GeV/ and 0.21–0.45 GeV/ otherwise. Fig. 2 shows an example of this subtraction process for two ranges in two centrality bins.
The scaled background distribution is then subtracted from the sameevent pair distribution. The subtracted result thus represents a sample of real s. The peak is fit to a Gaussian to determine its width and mean position. The raw yield of s is determined by integrating the counts in a window of around the mean. The width and mean are recorded and parameterized as a function of and centrality based upon this integrated, large sample. The positions and widths from this parameterization are then used when we extract the (much smaller) raw yields in bins of angle with respect to the reaction plane. The maximum variation of the yields (multiplicities) with is only about a factor of 2, and therefore the means and widths are not expected to change substantially. Furthermore, the statistics are much poorer in the bins, which would make individual fits unreliable.
There is a residual background in the invariant mass distributions even after the mixedevent distribution has been removed, especially at lower (below GeV/). This is due to correlations that event mixing cannot reproduce, like the “subevent structure” due to the presence of jets or multiple, closeby showers from an annihilating antineutron, or imperfections of the reconstruction algorithm, such as cluster merging, cluster splitting, and a host of other contributions. Much of the residual background is excluded by starting the fit at 0.09 . What is left is accounted for by including a firstorder polynomial in the fits to the (already backgroundsubtracted) invariant mass distribution, and subtracting its integral from the raw yield (see Fig. 2). In the more central events, the peak deviates slightly from gaussian on the high mass side, due to overlapping clusters. The use of ecore mitigates this effect, and the systematic uncertainty on yield extraction arising from the remaining asymmetry has been estimated to be 34% Adare et al. (2008a).
iii.4 Elliptic flow measurement
To obtain the azimuthal angle dependence of production, we measure raw yields in a given bin as a function of the angle with respect to the reaction plane orientation in six equallyspaced bins of covering the range . The yields are measured in each bin using the same procedure described in III.3 for the reactionplane inclusive measurement except that the mass fits are not performed in each bin. Instead, the peak integration window is set around the mean where the width and mean are taken from the inclusive analysis. The resulting raw angular distribution can then be fit to determine the strength of the modulation in the yield. Because the PHENIX BBCs have uniform azimuthal coverage, the measurements have uniform acceptance in when averaged over a large event sample, despite the limited azimuthal acceptance of the PHENIX electromagnetic calorimeters.
Assuming elliptic flow is the dominant source of variation in the yields, we perform a fit to the angular distributions of the form
(6) 
We use an analytic linear fitting procedure that matches the integral of Eq. 6 over each of the bins to the measured yield within the corresponding bin. In the definition of the function we account for nonzero covariances between the yields in the different bins resulting from the limited acceptance of the calorimeters. These covariances have been evaluated separately for each and centrality bin. Examples of the raw distributions and the results of the fits are shown in Fig. 3. The resulting values are then corrected upward to account for reaction plane resolution using correction factors described in Section III.2.
iii.5 () measurement
The nuclear modification factor has played a critical role in understanding energy loss mechanisms. is defined as
(7) 
where is the mean Glauber overlap function for the centrality being analyzed:
(8) 
from which the mean number of binary nucleonnucleon collisions can be calculated, .
For each bin, we can calculate the ratio
(9) 
where is the number of s observed in the given bin. Since the BBC is azimuthally symmetric the PHENIX acceptance has no dependence, there should be no azimuthal dependence to efficiency and acceptance corrections. As a result,
(10) 
Thus, we can use measured inclusive to convert to . Since the detector efficiency and acceptance corrections are already contained in , there is no need to apply them to .
Prior to calculating we correct the ratios for the finite reaction plane resolution using an approximate unfolding technique. For a pure flow distribution, we can express the influence of the resolution broadening on the measured distribution
(11) 
where according to the results of Section III.2 . Then, if the measured distribution resulted from pure elliptic flow, it could be corrected back to the true distribution by
(12) 
As shown above, the general features of the measured distributions are welldescribed by pure modulation. However, we wish to preserve in our measurements of the azimuthal dependence of the production the full shape of the measured distribution, including possible small nonelliptic contributions. For this purpose, the correction described in Eq. 12 applied to the data represents an approximation to a full unfolding procedure that becomes exact when the distribution is purely in form. We have checked for a few cases that a full unfolding procedure applied to the measured distributions using singular value decomposition regulation of the response matrix reproduces the correction in Eq. 12. From the corrected ratios, , we use Eq. 10 to obtain .
Iv Results
iv.1 Elliptic flow coefficient
The results of the measurements using the methods described in Sec. III.4 are presented in Fig. 4 as a function of for different centrality bins. The data points in the figure are plotted at the mean in bins of width for and for . The error bars shown on the data points were obtained by multiplying the raw fit errors (see Section III.4) by the same reaction plane resolution correction factor applied to the values themselves. The error bars, then, represent uncorrelated statistical errors on the measured values arising from statistical errors on the data points used in the fits (these errors would be categorized as Type A uncertainties in the framework described in Adare et al. (2008c) or uncorrelated). Systematic errors on the measurements due to the reaction plane determination procedure and from systematic uncertainties in the reaction plane resolution correction are represented in Figs. 4,5 by filled boxes, which for most data points are similar in size or smaller than the data points (these uncertainties are classified as Type B Adare et al. (2008c) or correlated).
Figure 5 shows for four centrality ranges, obtained by combining data from the centrality bins shown in Fig. 4. The corrected distributions from individual centrality bins were summed over a given centrality range and then fit to obtain the corrected values shown in Fig. 5. The reaction plane resolution correction produces correlated errors in the corrected distributions for each original centrality bin, and these correlated errors persist in the summed distribution. These correlated errors are not included in the statistical errors for the fit to the combined distribution, but their impact on the final value is estimated separately by evaluating the changes in the fit parameter that result from adding to the summed values of the correlated errors. Since this estimated uncertainty results from the statistical uncertainties on the values for the original centrality bins, we include the bounds obtained from this procedure in the statistical error on the values for the combined centrality bins. Systematic errors for the combined bins are plotted similarly to those in Fig. 4.
The results presented here nearly double the range of previous PHENIX measurements from RHIC
Run2 Adler et al. (2006). Those measurements are shown for comparison purposes in Fig. 5. Good agreement is seen between the Run4 measurements presented here and the Run2 results except in the 4060% centrality bin where the new measurements are systematically higher by . This difference is attributed to improved reaction plane resolution corrections for the 4060% centrality bin resulting from the combining of corrected distributions from smaller centrality bins. This summing procedure better handles the rapid variation of reaction plane resolution with centrality in midcentral to peripheral collisons. Furthermore, we have crosschecked the procedure using 5% bins, verifying the combined result reproduces the data analyzed in wider bins. The previous Run 2 analysis did not have a sufficiently large data sample to allow the use of separate 4050% and 5060% bins, and therefore the reaction plane resolution correction was necessarily less accurate. The measured values presented in Fig. 4 and Fig. 5 are also consistent with previously published PHENIX charged pion measurements Adler et al. (2003a, 2006).
The results in Figures 4 and 5 show a rapid increase of with increasing at low , a maximum in the range , and then at higher a decrease of with increasing . An increase in at low is wellestablished Adare et al. (2007); Adams et al. (2005b); Back et al. (2005b) and is understood to result from the collective elliptic flow of the medium generated by the initial spatial anisotropy of the collision zone. Hydrodynamical models have been successful in quantitatively describing the pion in the region . However, it has also been wellestablished that the pion deviates from the hydrodynamic prediction above 1.5 , a result that is understood to imply the contribution of hard processes, distortions of the spectrum due to recombination at freezeout, and/or effects from dispersive hadronic evolution after freezeout. Thus, a change in the variation of with near GeV/c is not unexpected. If the large values at lower are interpreted as resulting from soft, collective mechanismsm, then a decrease in for GeV/c suggested by the data in Fig. 4 would naturally reflect an increasing contribution of hard processes with smaller .
To statistically test the significance of the decrease of with , we show in Fig. 6 the results of linear fits to the high values for the 2030% and 3040% centralities. The panels on the lefthand side of Fig. 6 display a series of fits beginning at higher values, the first fit starting at the near the maximum , . The righthand panels show the 1 limits of the functions for the three fits in the series, calculated from the 1 variation of the two fit parameters (and including the covariance between them). The results of the fits indicate that the decrease of with at higher is statistically significant, though the data points for GeV/c are not sufficient by themselves to establish a trend. We can state, however, that the points for GeV/c are consistent with the linear decrease obtained including the lower points. A question we would like to answer, then, is whether the data show any indications of devation from a monotonic decrease in indicating the transition to a quenchingdominated azimuthal variation.
A complete understanding of over the measured range therefore requires the treatment of the transition from soft to hard dominated physics. According to the above discussion, in the range where is maximum, particle production is likely not dominated by hard processes and the reduction of with increasing indicates increasing hardscattering contributions (or decreasing soft contamination). Motivated by this general argument, we have attempted to describe the results in Figs. 4 and 5 by a functional form
(13) 
The first term is intended to describe a rapidly rising and saturating soft resulting from collective motion while the second term represents a rapidly falling soft/hard ratio. The additive constant in the second term represents an asymptotic that could describe a constant or slowly varying azimuthaldependent quenching. We show in Figs. 78 the optimum fits to the full set of values in the different centrality bins and the result of variation of the fit parameters taking into account the complete covariance matrix from the fits.
The fits to the data show that the measured dependence of the is qualitatively compatible with a description of the of low and intermediate region in terms of a collective flow modulation diluted by a decreasing relative soft contribution with increasing . Assuming this picture, it is then important to determine at what the contamination from the soft production no longer dominates the measured variation of yield. For most of the centrality bins, the fits suggest that decreases over most of the measured range albeit with a decreasing slope at higher . The central bins are compatible with continuing to decrease beyond the measured range although the uncertainty bands also accommodate saturating within the measured range. The more peripheral bins (3040% and 4050%) suggest that the has reached a nearly independent value by GeV/c. The 5060% centrality bin has sufficient fluctuations that little can be inferred from the dependence of in that centrality bin. In all of the centrality bins, the data are consistent with a smooth reduction of from a maximum to a nonzero value at high with that value increasing in more peripheral collisions as would be expected from quenching in an increasingly anisotropic medium. While the functional form in Eq. 13 can describe the variation of within the range of the current data and within the statistical fluctuations of the data points, it is possible that this description will fail at higher with improved statistics. In fact, a statistically significant deviation of from the form in Eq. 13 might provide the most direct evidence of the dominance of quenching effects in determining .
iv.2 Nuclear modification factor with respect to the reaction plane
The nuclear modification factor as a function of for six centrality bins is shown in Figs. 914. The closed circles represent the dependent measurements described in this paper while the open circles positioned at represent the inclusive measurement Adare et al. (2008a). In both cases statistical uncertainties (i.e. Type A) are represented by the error bars. For the inclusive measurement, the total systematic uncertainties (or Type C Adare et al. (2008c)) are shown by the boxes. The upper and lower ranges of the correlated statistical uncertainties (i.e. Type B) on the measurements resulting from the reaction plane resolution correction are shown by the (blue) solid and (red) dashed lines. For all bins except the 010% centrality bin a dotted line is plotted at for reference. We note that by construction, the average from the reaction plane dependent measurements must be equal to the inclusive .
The results in the Figs. 914 show that the inplane suppression is generally weaker and varies more rapidly with than the suppression for s produced at larger angles. As the collisions become more peripheral (for example, 5060%), the small suppression seen in the inclusive almost vanishes for s emitted close to the reaction plane. In a previous analysis Adler et al. (2007), it was observed that the inplane even exceeded unity for peripheral collisions; these data exhibit no such enhancement. However, the results presented in this article agree within systematic errors with previously reported data.
The results are combined in Fig. 15 that shows the dependence of the in each of the six bins included in this analysis. We use a semilog scale so that the reduction of the integrated in more central collisions does not confuse the interpretation of the results. For clarity, the results from the 2030%, 3040%, 4050%, and 5060% centrality bins are shown on linear plots in Fig. 16.
The results exhibit a peak near 2 , which becomes more prominent for more central collisions. The peak is strongest in the 010% bin where there is little modulation of the distributions at low or high , so the peak cannot be directly attributed to elliptic flow. The peak in near 2 is much weaker in the more peripheral (4050% and 5060%) centrality bins, particularly for s produced at larger , and the primary variation seen in these peripheral bins with increasing is a reduction in that is only weakly dependent.
For the intermediate centrality bins (1020% through 3040%) the peaking in is seen in all bins, but is much stronger in the inplane bins. For these intermediate centralities and for values above the peak in ( ), the for s produced at angles normal to the reaction plane is nearly constant with while the for s produced at small angles from the reaction plane decreases rapidly with increasing . The near constancy of the outofplane together with the rapid reduction in inplane indicates that in the intermediate centrality bins, the and inclusive decrease simultaneously with increasing such that is approximately constant. We will argue below that a correlation between and may naturally result from the underlying physics responsible for the azimuthal variation of the particle yields. However, we observe that a simultaneous reduction in integrated and suggested by the more central data would be contrary to naive energy loss expectations since smaller would imply stronger quenching in the medium which would, in turn, imply larger variation between inplane and outofplane quenching.
A similar implicit correlation between integrated and is seen in the centrality dependence of the results. These are replotted in Fig. 17 as a function of for three bins – the bins closest to and further from the reaction plane and one of the intermediate bins. For , the outofplane values are nearly independent of centrality while the inplane values decrease rapidly with increasing centrality. This result would have a natural geometric interpretation for production dominated by hard scattering and jet quenching. The length of the medium normal to the reaction plane varies only slowly with centrality except in the most peripheral collisions. Then, if the suppression is determined primarily by the path length of its parent parton in the medium, the would be nearly constant. Following the same argument, the yield for pions in the direction of the reaction plane would be much less suppressed in noncentral collisions due to the short path lengths of the parent partons in the medium. However, with increasing centrality and decreasing anisotropy of the collision zone, the inplane parton path lengths would grow to match those in the outofplane direction. Thus, if the suppression depended primarily on path length, the inplane would naturally drop to match the outofplane values reproducing the behavior of Fig. 17. In order to better see the difference between the in and outofplane behaviors, these data are also plotted on Fig. 18 with a semilog scale.
One difficulty with this geometric interpretation of the results given above is that the trend in the data that it is supposed to explain persists down to low , where the values are too large to be accounted for via perturbative or formation time based energy loss scenarios Shuryak (2002); Pantuev (2007); Drees et al. (2005); Hirano and Nara (2004); Renk et al. (2007); Majumder et al. (2007). That fact coupled with the pronounced peaking in near 2 suggests that physics other than hard scattering and jet quenching must be invoked to explain the yields at intermediate . However, the fact that the outofplane yields show less pronounced peaking near 2 , that they vary little as a function of above 3 , and that they vary little with centrality for could be interpreted to imply that the suppression at angles normal to the reaction plane more directly represents the effects of quenching of hard quarks and gluons while the yield of s produced more closely aligned with the reaction plane is enhanced by the collective motion of the system. That additional enhancement could either be due to soft hadrons being boosted to larger values by the collective elliptic flow or could result from weaker quenching for partons moving in the direction of the flow field Armesto et al. (2004); Renk et al. (2007). Simultaneous description of the inplane and outofplane behavior is a sensitive test of energy loss models.
The values presented in Figs. 78 also peak near 2 , but the locations of the maxima in are shifted to higher than the maxima in . This suggests that the two effects may not directly related, but we observe that the maxima in the distributions in Figs. 1516 shift to larger for smaller . To better illustrate the shift of the maxima in we show in Fig. 19 the values for the different bins and indicate the variation of the peak position obtained using polynomial fits to the first four bins. For the 3040% centrality bin, the maximum in for is shifted by 0.4 relative to the bin. This shift in the peak with can produce a that peaks higher in than the inclusive .
The observed shift in the peak of with illustrates an important property of collective motion of the medium. The collective motion does not simply superimpose azimuthal variation on the particles produced at a given , it provides a dependent shift and/or broadening in the transverse momentum spectrum of the produced particles. The resulting distortion will be the smallest for particles produced at angles normal to the reaction plane and will be largest for particles produced in the plane. Any collective shift of soft particles to higher will increase the measured for small relative to large values producing a simultaneous increase in both the integrated and the . With increasing , an expected decrease in the soft contamination would naturally explain the simultaneous reduction in and integrated evident in the 1020% and 2030% bins where the separation between the curves for different bins decreases while the average also decreases. We will return to investigate this correlation again below.
The 4050% and 5060% centrality bins in Fig. 15 show little of the peaking near 2 , especially in bins not aligned with the reaction plane. Nonetheless the values for the more peripheral bins reach the same large maximum values, , at intermediate as the values for more central bins where the peak in is more prominent. Thus, while the peaking in is less prominent in the more peripheral bins, the relative difference between the inplane and outofplane yields in the 4050% centrality bin is comparable to that in the 2030% centrality bin. However, it is possible that the large dependence in the more peripheral bins and the apparent persistence of that variation to high in the 4050% centrality bin more directly reflects the larger spatial anisotropy of the collision zone in more peripheral collisions. The question of whether the suppression measurements presented here can be understood on the basis of geometry and jet quenching will be more fully explored in the following section.
We have observed above that the and centrality dependence of indicate a correlation between inclusive and such that the outofplane yields vary only slowly with or centrality while the inplane yields approach the outofplane yields with increasing or increasing . Such a correlation between these two seemingly unrelated quantities merely indicates that the yields or of ’s measured inplane and outofplane more directly reflect the underlying physics responsible for the azimuthal variation than the integrated yield or and the amplitude of the modulation, . Indeed, we have argued above that at higher the centrality dependence of may reflect the geometry of jet quenching. At intermediate , the results suggest contamination of the inplane yields by soft production and a simultaneous decrease in and with increasing as the relative contribution of collective soft processes to production decreases. To more directly demonstrate the correlation that forms the basis of these arguments we show in Fig. 20 a plot of versus the inclusive for centralities from 0 to 60%. Data are displayed for . The intermediate centrality bins show a correlated increase of and consistent with the discussion above and a possible saturation of for larger values. In fact, the trends for different centrality bins appear to be in general agreement. However, their exact relationship and establishing or excluding a causal connection requires further investigation.
iv.3 Nuclear modification factor dependence on path length
The centrality of a collision fixes the geometry of the overlap region between the nuclei, and fixing the angle of emission of the particles further constrains the path length through the medium. We can use this feature to study the dependence of the nuclear modification factor on the path length traversed by the partons. We investigate the path length dependence by expanding on several methods previously described in Adler et al. (2007). We start with three estimators of the path length that are purely geometric, and one that includes the color density of the medium in its calculation:

We start by modeling the overlap region as an ellipse defined by
(14) where the minor axis is oriented in the direction and is parallel to the reaction plane. The axes and are fixed by the intersection in the transverse plane of two hard spheres with fm. In terms of the spatial eccentricity (often used with Glauber calculations), we can express the distance from the origin to the edge of the ellipse at a given angle:
(15) Since this length starts at the origin, and does not take into account color density, the expression provides a very simple estimator with which we can evaluate the dependence of the on path length. We will refer to the hard sphere result as .

Instead of an ellipse strictly defined by the transverse profile of two hard spheres, we model the collision region as an effective ellipse whose dimensions are determined by equating the RMS radius and eccentricity to the corresponding quantities calculated from the transverse distribution of participant density based on standard Glauber calculations. This length, , is evaluated using the same expression as Eq. 15, with . Both quantities are determined using the PHENIX Glauber model Miller et al. (2007).

For a more realistic approach, we evaluate the distance along the parton’s path weighted by the participant density,
(16) where is the hardscattering position and is the angle of the jet with respect to the axis. The jet production point is sampled from a Monte Carlo using a weighted distribution and a uniform distribution. The participant density, assumed to be proportional to the color density, is calculated from the Glauber model. The density in Eq. 16 is modeled using a 1D Bjorken expansion,
Thus roughly represents LPM energy loss Aurenche et al. (2000). Note that approaches the same dependence as the standard but differs from by a factor of 2 at (additionally this form is regular at ).

Finally, we modify by normalizing it by the value of the participant density at the center of the collision region, . As a result, is an estimator based on geometry alone, but also accounts for the effect of the density distribution both on the jet production point as well as the path from that point to the edge of the medium.
Figs. 2124 show the dependence on , , , and respectively. The results shown in this paper cover the range  , not only extending the measurement presented previously but allowing a much finer binning in . The statistical errors in the measurements are represented by error bars (see Section IV.2). The systematic errors shown in these data are on the values only, and are indicated by the filled boxes. The major contribution to the systematic error in the values is the determination of , and is at the 1020% level.
Both the and behavior show an interesting degree of scaling. This result is all the more unexpected because of the overly simple geometric picture they represent. Despite the simplistic picture they present, they both exhibit striking universality: the hard sphere scales well in the low region (as high as ) while scales well to higher , at least one bin in beyond (though one might argue qualitatively this trend continues even higher when the most peripheral centrality is excluded). The more precise dependence available in this data set reveals a slight deviation from the universality with that was previously reported Adler et al. (2007).
By contrast, we expect the estimator to provide a somewhat more intuitive and concrete picture, as it represents a more realistic approach to the geometry and medium. Since we expect radiative energy loss to play a greater role at high , this should be the estimator that would provide the best scaling. In fact, at the higher range, a universality does emerge, though not until . Below that value, the measured points lie on parallel, but separate, curves. When is normalized to the central density, data again exhibit a more universal dependence over a wider range than what is seen with alone, as shown in Fig. 24.
When considered together these results offer a rich picture. At low to moderate simple geometry may play a larger role in determining the final level of than conventionally thought. At higher the scaling motivated by energy loss () describes the data well. We note that there are three possible (and not necessarily exclusive) interpretations: 1) at low to moderate the combined effects of the boost due to expansion and fragmentation are sensitive primarily to the difference in lengths traveled by the partons, and only weakly dependent on other parameters 2) we need to restrict the analysis of the to to be in the range where fragmentation followed by energy loss dominates, or 3) the assumption that energy loss does not depend linearly on color density Liao and Shuryak (2008) is incorrect and leads a departure from scaling with at low to moderate .