Measurement of two-particle correlations with respect to second- and third-order event planes in Au$+$Au collisions at $\sqrt{s_{_{NN}}}=200$ GeV

We present measurements of azimuthal correlations of charged hadron pairs in $\sqrt{s_{_{NN}}}=200$ GeV Au$+$Au collisions after subtracting an underlying event using a model that includes higher-order azimuthal anisotropy $v_2$, $v_3$, and $v_4$. After subtraction, the away-side ($\Delta\phi\sim\pi)$ of the highest transverse-momentum trigger ($p_T>4$ GeV/$c$) correlations is suppressed compared to that of correlations measured in $p$$+$$p$ collisions. At the lowest associated particle $p_T$, the away-side shape and yield are modified. These observations are consistent with the scenario of radiative-jet energy loss. For the lowest-$p_T$ trigger correlations, an away-side yield exists and we explore the dependence of the shape of the away-side within the context of an underlying-event model. Correlations are also studied differentially versus event-plane angle $\Psi_n$. The angular correlations show an asymmetry when selecting the sign of the trigger-particle azimuthal angle with respect to the $\Psi_2$ event plane. This asymmetry and the measured suppression of the pair yield out of plane is consistent with a path-length-dependent energy loss. No $\Psi_3$ dependence can be resolved within experimental uncertainties.

We present measurements of azimuthal correlations of charged hadron pairs in √ s N N = 200 GeV Au+Au collisions after subtracting an underlying event using a model that includes higherorder azimuthal anisotropy v2, v3, and v4. After subtraction, the away-side (∆φ ∼ π) of the highest transverse-momentum trigger (pT > 4 GeV/c) correlations is suppressed compared to that of correlations measured in p+p collisions. At the lowest associated particle pT , the away-side shape and yield are modified. These observations are consistent with the scenario of radiativejet energy loss. For the lowest-pT trigger correlations, an away-side yield exists and we explore the dependence of the shape of the away-side within the context of an underlying-event model. Correlations are also studied differentially versus event-plane angle Ψn. The angular correlations show an asymmetry when selecting the sign of the trigger-particle azimuthal angle with respect to the Ψ2 event plane. This asymmetry and the measured suppression of the pair yield out of plane is consistent with a path-length-dependent energy loss. No Ψ3 dependence can be resolved within experimental uncertainties.

I. INTRODUCTION
Energy loss of hard-scattered partons (jet quenching [1]) resulting from the interaction of a colored parton in the quark gluon plasma (QGP) formed in relativistic heavy ion collisions [2][3][4][5] has been observed in several different ways. Suppression of single-particle and single-jet invariant yields in central A+A collisions [6][7][8][9][10][11] provide a baseline measurement of jet quenching. Measurements of correlations between two particles and/or * Deceased † PHENIX Spokesperson: akiba@rcf.rhic.bnl.gov jets give more detailed spatial information of the jet quenching process inside the medium [12,13]. The first jet suppression effect observed in these correlations was an attenuation of the away-side yields in high-transverse momenta (p T ) correlations in the most central Au+Au collisions at √ s N N = 200 GeV [14]. The centrality dependence of high-p T π 0 -hadron correlations [15] shows a monotonic attenuation of the away-side yields with increasing propagation length of partons through the medium. In addition to away-side yield suppression, direct photon-hadron correlations [16][17][18], two-particle correlations [19][20][21], and jet-hadron correlations [22][23][24][25] show that low-momentum particles correlated with highp T jets are enhanced in yield, especially at large angles with respect to the jet axis. This may be attributable to the radiation from the parent parton or other lost energy recovered by the nearby medium. Thus, two-particle angular correlations have provided much of the experimental information we have about jet energy loss [15,[26][27][28][29][30].
It is important to understand the interactions of partons with the QGP at all scales from the hard-scattering scale to the thermal scale. Below E jet ≈10 GeV, full jet reconstruction is much more difficult due the underlying event subtraction. Two-particle correlations are important because they can probe lower energies. However, observations of many of the energy loss effects mentioned above, especially for lower jet and particle momenta in dijet two-particle correlations, have been obscured due to the much larger contribution from the underlying event at these momenta.
The underlying event modulations are attributed to hydrodynamic collective flow patterns where the importance of higher-order terms was established more recently [31][32][33][34][35][36]. These patterns are thought to result from the hydrodynamic response of the QGP to fluctuating initial geometrical shapes of the interacting portions of the colliding nuclei. Many hydrodynamic models have been developed which capture these effects [37][38][39], but to date, important details of these models are still under development, and their full implementation requires involved calculations. This motivates the use of a simpler data-driven model, which will be explored in this work.
The shape of the collective flow in the transverse plane is parameterized [40][41][42][43] by a Fourier expansion with v n = cos n(φ − Ψ n ) (1) where v n is the n th -order component of single-particle correlation due to flow, φ is the azimuthal angle of emitted particles, and Ψ n is the event plane angle defined by the n th -harmonic number. For the first decade of the Relativistic Heavy Ion Collider (RHIC), only the even harmonics and frequently only the n = 2 term, were considered. The shapes of two-particle correlations after subtraction of the n = 2-only background motivated the introduction of the other harmonics, most importantly n = 3 [31,[33][34][35]44]. Under the two-source (flow + jet) model assumption [45], this underlying event is directly subtracted to obtain the jet contributions. In our previous measurements and most RHIC A+A results, the subtracted flow modulations of the underlying event were limited to contributions of v 2 and the fourth-order harmonic component with respect to the second-order event plane v 4 {Ψ 2 } [14,15,19,22,27,[45][46][47][48][49]. Only the recent STAR measurement [50] took into account contributions from v 3 and the fourth-order harmonic component uncorrelated to the second-order event-plane v 4 {uc}. At low to intermediate p T in two-particle correlations, intricate features appear such as the near-side longrange rapidity correlations called the "ridge" [46,51] and the away-side "double-humped" structures [27, 45, 47-50, 52, 53]. Across the large rapidity ranges available at the Large Hadron Collider, the rapidity-independence and hence the likely geometrical origin of most of these structures have been established. Experiments have shown that the ridge and the double-hump structures in the two-particle azimuthal correlations for |∆η| > 1 for ALICE and |∆η| > 2 for ATLAS and CMS measured in p+p, p+Pb, and Pb+Pb collisions at √ s N N = 2.76 and 5 TeV [33,35,54] are the same in shape and size at much larger rapidity differences. Both the ridge and hump are successfully explained by the higher-order harmonics. Despite this conclusion, how the jet correlations combine with the flow correlations, especially at small ∆η, to yield the total two-particle correlation has not been quantified. In particular, the correlations left after subtracting a flow-based model at small ∆η have not been analyzed in detail.
In this work, we assume a two-source model where the total pair yield is a sum of a jet-like component and an underlying-event component. The underlyingevent components is modeled using the flow harmonics v n (n = 2, 3, 4), event plane resolutions, and the most important event plane correlations between Ψ 2 and Ψ 4 . We assume that the v n measured through the event plane method are the the same as those in the correlation functions. Event-by-event v n fluctuations [55], v n − v m correlations between different order [56], and rapidity dependent event plane decorrelations [57] are not considered in this background model. Measurements of v 1 using the event plane method [58][59][60][61] generally yield v 1 ∼ 0 at η = 0 integrated over all p T . Finite values of p T differential measurements of v 1 (p T ) [60] include momentum conservation and jet (mini-jet) effects which are considered signal in two-particle correlations. We therefore excluded contributions from v 1 {Ψ 1 } and event plane correlations involving Ψ 1 from the background model. For the inclusive trigger correlations, we estimated a potential impact of v 5 modulation using an empirical relation v 5 ∼ 0.5×v 4 found in ATLAS v n measurements [33].
After subtracting the underlying event with the model, we study the structures observed at high p T where the flow backgrounds are negligible. Because the jet signalto-flow background is significantly reduced in the low to intermediate p T region, studying the correlations there provides a more stringent test of such a background model. Any features left in the residuals can be used to reveal jet energy-loss effects at low and intermediate p T . However, because of our simple model, only substantially significant correlations can be attributable to the medium effect on jets (i.e. broadening or suppression) or the medium response (i.e. yields at large angles from the jet).
We also study two-particle correlations measured differentially with respect to Ψ 2 and Ψ 3 event planes as depicted in Fig. 1. Such differential correlations probe the path-length and geometrical dependence of energy loss with more event-by-event sensitivity and also extend similar studies of high-p T correlations [15] down to lower-p T . We also use a new method of distinguishing "left/right" asymmetry in the Ψ n correlations, which can be em-FIG. 1. Schematic picture of a trigger particle selection with respect to event-planes and pairing a trigger particle with an associate particle.
ployed to provide more information on the background dominated low and intermediate p T regions.
In this article: Section II describes the detector setup of the PHENIX Experiment. Sections III A, III B, and III C describe the analysis methodology of particle selections, higher-order flow harmonics, and two-particle correlations, respectively. Section IV presents analysis results and discusses their interpretations. This section first starts with the highest p T trigger selections, p T > ∼ 4 GeV/c, and makes connections to known energy-loss effects. Next, lower trigger correlations down to 1 GeV/c are presented. Finally the event-plane dependence of the intermediate p T selections are investigated. Section V summarizes this article.

II. PHENIX DETECTOR
The PHENIX detector [62] was designed to measure charged hadrons, leptons, and photons to study the nature of the QGP formed in ultra-relativistic heavy ion collisions. Figure 2 shows the beam view and side view of the PHENIX detector including all subsystems for this data taking period.
The global detectors, which include the beam-beam counters (BBC), the zero-degree calorimeters (ZDC), and the reaction-plane detector (RXN), were used to determine event characterizing parameters such as the collision vertex, collision centrality, and event-plane directions. They are located on both the south and north side of the PHENIX detectors. The BBC is located at ±144 cm (3 < |η| < 3.9) from the beam interaction point and surrounds the beam pipe with full ∆φ = 2π azimuthal acceptance. Each BBC module comprises 64 quartzČerenkov radiators equipped with a photomultiplier tube (PMT) and measures the total charge (which is proportional to the number of particles) deposited in its acceptance. The BBC determines the beam collision time, beam collision position along the beam axis direction, and collision centrality. The ZDCs [63], located at 18 m away from the nominal interaction point, detect the energy deposited by spectator neutrons of the two colliding nuclei. The PHENIX minimum-bias trigger is provided by the combination of hit information in the ZDC and BBC, which requires at least one hit in both the ZDC modules and two hits in the BBC modules. The orientation of higher-order event planes is determined by the BBC and the RXN [64], which have different η acceptance. The RXNs are located at ±38 cm from the beam interaction point and have two rings in each module; RXN-inner and RXN-outer are installed to cover 1.5 < |η| < 2.8 and 1 < |η| < 1.5, respectively. Each ring has 12 scintillators in its azimuthal angle acceptance ∆φ = 2π.
Charged hadron tracks are reconstructed in the PHENIX central arm spectrometer (CNT), which is comprised of two separate arms, east and west. Each arm covers |η| < 0.35 and ∆φ = π/2. The PHENIX tracking system is composed of the drift chamber (DC) in addition to two layers of pad chambers (PC1 and PC3) in the east arm and three layers of pad chambers (PC1, PC2, and PC3) in the west arm. Momentum is determined by measuring the track curvature through the magnetic field by means of a Hough transform with hit information from the DC and PC1 with a momentum resolution of δp/p = 1.3% ⊕ 1.2%p [65]. Additional track position information is provided by the outer layers of the pad chambers and the electromagnetic calorimeter (EMCal).
The ring imagingČerenkov counter (RICH) and the EMCal identify and exclude electron tracks from the analysis. The RICH produces a light yield for electrons with p T >30 MeV and for pions with p T >5 GeV, meaning that a signal in the RICH can be used to separate electrons and pions below 5 GeV. Above 5 GeV where this is no longer possible, the energy deposited in the EMCal can be used for this separation. Electrons will deposit much more of their total energy than pions will, so that the ratio of deposited energy to track momentum is significantly higher for electrons than for pions.

III. ANALYSIS METHODOLOGY
The results presented are based on an analysis of 4.38 billion minimum-bias events for Au+Au collisions at √ s N N = 200 GeV recorded by the PHENIX detector at RHIC in 2007.

A. Particle Selection
Charged hadrons are selected from candidate tracks using cuts similar to previous correlation analyses [19]. One important cut to reject fake tracks, especially decays in the central magnetic field before the drift chamber, is an association cut to outer CNT detectors. The track trajectories are projected onto outer CNT detectors. The nearest hits in the PC3 and the EMCal from the projections are identified as hits for the track. The distributions of the distance in the azimuthal (φ) and beam (z beam ) directions between the hits in the PC3 and the EMCal and the extrapolated line are fitted with a double-Gaussian. One Gaussian arises from background and the other from the signal. Hadron tracks are required to be within ±2σ of the signal Gaussian mean in both the φ and z beam directions in both the PC3 and the EMCal. To veto conversion electrons, tracks with p T <5 GeV/c having one or moreČerenkov photons in the RICH are excluded from this analysis. For p T >5 GeV/c, we require E EMCal > 0.3 + 0.2c×p T GeV [48,66], where E EMCal is the cluster energy associated with the track.

Event-plane and Resolution
Each event plane Ψ n is determined event-by-event for different harmonic numbers n using the RXN and BBC detectors. The RXN detectors are used to measure the nominal values of v n while the BBC detectors provide systematic checks to the extracted v n values. The observed event-plane Ψ obs n is reconstructed as Here Q n,x and Q n,y are the flow vectors where φ i is the azimuthal angle of the i-th segment in the event-plane detector and w i is the weight proportional to multiplicity in the i-th segment. The k × n th-order resolution of n th-order event plane is defined as Res{kn, Ψ n } = cos kn(Ψ obs n − Ψ n ) and can be expressed as [42] Res{kn, where χ n = v n √ 2M , M is the multiplicity used to determine the event-plane Ψ n , and I k is the modified Bessel function of the first kind.
Because the north (N) and south (S) modules of a given event-plane detector have the same pseudorapidity coverage and see the same multiplicity and energy for symmetric nucleus-nucleus collisions, the north and south modules should have identical resolution. We obtain the event-plane resolution of an event-plane detector using the two subevent method [42].
The north+south combined event-plane resolution is determined from Eq. (5) with χ n = √ 2χ N,S n . The factor of √ 2 accounts for twice the multiplicity in north+south compared to north or south. Figure 3 shows the the north+south combined event-plane resolution for both RXN and BBC.

vn measurements
Higher-order flow harmonics v n [31,40,42] are measured by the event-plane method [42]. Charged hadron tracks with azimuthal angle φ are measured with respect to the event plane angle Ψ obs n . The flow coefficients v kn are measured as an event-average and track-average and correcting by the event plane resolution. As an example, v n in 20%-30% central collisions measured by different RXN event-planes are shown in Fig. 4 (a, d, g, j). The (blue) band indicates σ RXN . To evaluate the systematic uncertainty due to track matching, the matching cut was varied by ±0.5σ from the nominal 2σ window. We calculated the uncertainty σ mat as the average deviation between the v n with the nominal cut and the varied cut.
The variation due to the track matching cut is illustrated in Fig. 4 (b, e, h, k) by showing the v n in 20%-30% central collisions measured with tracks having a matching cut of 1.5σ, 2σ and 2.5σ. The differences between the nominal 2σ and both 1.5σ and 2.5σ are also shown and scatter around zero indicating the size of σ mat . The systematic uncertainties associated with the rapidity gap between particles and the event-plane σ rap are defined by the absolute difference between v n determined by the RXN average and v n determined by the BBC.
The v n measured with the RXN, the BBC, and their difference are shown in Fig. 4 (c, f, i, l). Except in the case of v 4 , this systematic uncertainty is much less than the uncertainty due to the RXN module variation. The small uncertainty in the rapidity gap supports the claim that, even though the rapidity gap between the central arm and the RXN is less than one unit of pseudorapidity, the contamination from nonflow correlations does not dominate the uncertainty on the extraction of v n . The total systematic uncertainties σ vn are the quadrature sum of these individual systematic uncertainties.
These total systematic uncertainties are applied symmetrically for a conservative estimate of the systematic uncertainties in the correlation measurements. In nearly all p T and centrality selections, the RXN systematic uncertainty dominates the total uncertainty. The v n results are shown in Fig. 5 and compared with previous PHENIX v n measurements Ref. [34]. They are consistent within uncertainties where they overlap. For the two-particle correlations, we use larger p T bins (four bin selections) as indicated in Table I, which lists the v n results. A test of scaling between higher-order flow harmonics and the second-order flow harmonics v n /v n/2 2 is in Table II similar to that done in Ref. [33]. These v n /v n/2 2 ratios are constant as a function of p T within systematic uncertainties, which are dominant over the statistical uncertainties in 0%-50% central collisions. A possible interpretation of this relation is the acoustic scaling described in Ref. [67].

Pair Selections
Tracks that pass the single-particle cuts are paired for correlations. In real events two tracks cannot come arbitrarily close together. The tracking algorithm would split or merge the tracks. Therefore, there is an acceptance difference in real and mixed events. These effects are estimated by calculating the distances ∆φ(rad) and ∆z beam (cm) between hits in the PC1 and the PC3, where ∆φ(rad) is the relative azimuthal angle and ∆z beam (cm) is the relative length between two track hits in both real and mixed events. The ratios of these distributions are shown in Fig. 6. The ratio is normalized to one where the detectors work well to reconstruct hadrons, and inefficient and over-efficient regions deviate from unity. The dashed lines indicate the cuts used to remove the inefficient regions: (∆φ P C1 /0.04) 2 + (∆z beam,P C1 /90) 2 < 1, (∆φ P C1 /0.08) 2 + (∆z beam,P C1 /8.0) 2 < 1, (12) (∆φ P C3 /0.07) 2 + (∆z beam,P C3 /25) 2 < 1.

Inclusive Trigger Correlations
Two-particle correlations are calculated as where ∆φ = φ a − φ t is the relative azimuthal angle between trigger and associated hadrons and N real (∆φ) and N mix (∆φ) are pair distributions in the real and mixed events, respectively. N real (∆φ) reflects the physical correlation among trigger and associated hadrons from jets and from the underlying event as well as the dihadron detector acceptance effects. N mix (∆φ) is obtained by pairing trigger and associated hadrons from randomly  selected pairs of events that have similar collision vertices and centralities so that it reflects only the dihadron acceptance effects. Taking the ratio between the real and mixed distributions corrects for the nonuniform azimuthal acceptance for dihadrons so that C(∆φ) contains only physical effects. The correlation function is normalized so that the underlying event modulates approximately around unity. Within the two-source model [45], the correlation function C(∆φ) is composed of a jet-like term J(∆φ) and an underlying-event term that includes modulations from flow F (∆φ). We use the following model for the underlying event [42] n v a n cos n∆φ.
The jet-like correlation is then obtained by subtracting F (∆φ) from C(∆φ) as The scaling factor b zyam is determined with the zero yield at minimum (ZYAM) method [45,68,69]. In the ZYAM assumption, F (∆φ) is scaled such that J(∆φ) has a minimum of exactly zero. This therefore gives the lower boundary of possible jet-like correlations. The ZYAM scaling factor b zyam is determined by fitting the correlation function C(∆φ) with Fourier series for − π 2 < ∆φ < 3π 2 and identifying the single point where this fit and F (∆φ) have the minimum contact point. The statistical uncertainty e zyam of the ∆φ bin containing the ZYAM point is used to scale F (∆φ) to estimate the systematic uncertainty due to ZYAM An example ZYAM identification is shown in Fig. 7. The jet-like correlations J(∆φ) are scaled to the pertrigger yield Y (∆φ) (17) where N t is the number of trigger hadrons, N ta is the number pairs, and ǫ a is the single-hadron tracking efficiency in the associated hadron p T range. The efficiency is estimated via detector simulations for acceptance and occupancy effects as discussed in Ref. [7,15,27,48]. The ). The columns represent centrality bins 0%-10% (a,f,k,p), 10%-20% (b,g,l,q), 20%-30% (c,h,m,r), 30%-40% (d,i,n,s), and 40%-50% (e,j,o,t).
Coefficients obtained in this analysis are shown by blue points and those measured in Ref. [34] are shown by magenta points. Shaded bands and magenta lines indicate systematic uncertainties of those measurements.
tracking efficiency of a trigger particle is canceled by the ratio of d∆φN real /N t .

Event Plane-Dependent Correlations
Event plane-dependent two-particle correlations C(∆φ, φ s ) are defined as where φ s = φ t − Ψ n and N real (∆φ, φ s ) and N mix (∆φ, φ s ) are the event plane-dependent pair ∆φ distributions in real and mixed events, respectively. Similar to inclusive correlations, event planedependent jet-like correlations J(∆φ, φ s ) are obtained by subtracting the event plane-dependent flow background term F (∆φ, φ s ) from C(∆φ, φ s ) with a ZYAM scale factor as We use the same b zyam as determined from the inclusive correlations from the same trigger, associated and centrality selection. An analytical formula for F (∆φ, φ s ) including the n = 2 event plane dependence exists [70], however, it is not easily applied with finite correlations between the n = 2 and n = 4 event planes. For this reason, a Monte Carlo simulation is employed to estimate F (∆φ, φ s ). This is described in Sec. III C 4 below.
The event plane-dependent jet-like correlations are converted into event plane-dependent per-trigger yield as (20) where N t φs is the number of trigger hadrons and N ta φs is the number of pairs in the trigger event-plane bin.

Flow Background Model Including Event Plane Dependence
With the assumption that the measured v n from the event plane method are purely from collective dynamics of the medium, flow-like azimuthal distributions of single hadrons can be generated by performing a Monte Carlo simulation inputting the experimentally measured v n , the resolution of the event planes, and the strength of correlation among different order event planes. The single-hadron azimuthal distributions due to collective flow can be described by a superposition of v n as where φ is the azimuthal angle of the emitted hadrons and Ψ true n is a true nth-order event plane defined over [−π/n, π/n]. Separate distributions using v n for each p T ranges of trigger and associated particles are used in the simulation. The trigger and associated distributions in real events share a common Ψ n , while those in mixed events do not. The experimental event plane resolution is introduced through a dispersion term ∆Ψ n where Ψ obs n = Ψ true n + ∆Ψ n . We calculate ∆Ψ n as where z n = χ n / √ 2 cos(n∆Ψ n ) and erf(z n ) is the error function [42,71]. This equation can be solved for ∆Ψ n by using the experimentally-determined χ n from the measured event plane resolutions using Eq. (5).
Because a weak correlation between Ψ true We estimate ∆Ψ 42 assuming the correlation between the two event planes follow similar functional forms as the dispersion of event planes due to the resolution. That is, we assume, where assumed to be similar to Eq. (5) where cos(4∆Ψ 42 ) = v 4 {Ψ 2 }/v 4 [72]. The functional shape of Eq. (23) is verified by event plane correlation studies using the BBCs and the RXNs following the method described in Ref. [73]. The correlation strength between Ψ true 2 and Ψ true 3 , cos 6(Ψ 2 − Ψ 3 ) , is measured to be consistent with zero within large statistical uncertainties. Potential impacts of cos 6(Ψ 2 − Ψ 3 ) to the event plane dependent correlations are estimated using the value of cos 6(Ψ 2 − Ψ 3 ) reported in Ref. [73] by the ATLAS Experiment. The impact of cos 6(Ψ 2 − Ψ 3 ) is within the systematic uncertainties described later.
The event plane-dependent background shapes are determined by generated particles in this simulation using Eq. (18). Figure 8 shows event plane-dependent correlations and backgrounds with a selection of the absolute trigger azimuthal angle relative to the event-planes |φ t − Ψ n |. The backgrounds agree with the experimental correlations except at ∆φ = 0, π where contributions from jets are expected. Figure 9 shows event plane-dependent correlations and backgrounds with a selection of trigger azimuthal angle relative to event planes φ t − Ψ n < 0. Agreement between the experimental correlations and the background except at ∆φ = 0, π is also observed here. Other event-plane dependent correlations and backgrounds with a selection of trigger azimuthal angle relative to event-planes φ t − Ψ n < 0 for different collision centralities and p t,a T selections are shown in the Appendix.

D. Unfolding of Event Plane-Dependent Correlations
In this analysis, φ s is divided into 8 bins. The width of the φ s bins is π/8 and π/12 when correlating with Ψ 2 and Ψ 3 , respectively. The event plane-dependent per-trigger yields Y (∆φ, φ s ) are smeared across neighboring event plane bins due to limited experimental resolution of the event planes. We unfold the smearing to obtain the true event-plane dependence of the correlations. Two different methods are used to check the unfolding procedure: (I) iterative Bayesian unfolding, Y itr unf , and (II) correcting the event plane-dependence of the per-trigger yield based on a Fourier analysis, Y fit unf .

Iterative Bayesian Unfolding
The Iterative Bayesian Unfolding Method presented in Ref. [74,75] is applied to this analysis with the following formulation wheren(φ s , ∆φ) is the unfolded distribution, n obs (φ s , ∆φ) is the experimentally observed distribution, n(φ s , ∆φ) is the prior distribution, P (φ s,j |φ s,i ) is the conditional probability matrix where φ s,i is measured to be φ s,j , and ε i = j P (φ s,j |φ s,i ) is the efficiency. In the iterative calculation,n(φ s , ∆φ) also serves as the prior distribution of the next loop. We perform this unfolding separately for every ∆φ bin.
We define the experimentally observed distribution as n obs (φ s , ∆φ) = 1 + Y (φ s , ∆φ) using the measured event plane-dependent per-trigger yield. The offset is to prevent a divergence in the iteration due to small yields near the ZYAM point. In the initial loop of the iteration, we define the prior distribution as n(φ s , ∆φ) = n obs (φ s , ∆φ).
The probability distribution of the relative azimuthal angle between the true event plane Ψ n and the measured event plane Ψ obs n can be translated into the difference between real and observed φ s as of Ψ3. Schematic pictures in each panel also depict these ranges of the trigger particle selections with respect to event plane Ψn.

FIG. 9.
Event plane-dependent C(∆φ) (black circles) and event planedependent model flow background (blue lines) of (2 < pT < 4)⊗ (1 < p a T < 2) GeV/c. Trigger particles are selected in (a) out- Schematic pictures in each panel also depict these ranges of the trigger particle selections with respect to event-plane Ψn.
With this probability distribution of φ obs s − φ s , the probability matrix P (φ s,j |φ s,i ) is determined by the degree of the contamination by neighboring φ s bin as where s n (n = 0) is the contamination fraction from the n-th φ s bin away from a selected φ s bin, and s 0 is the fraction of the true signal in the selected φ s bin. A study in previous identified particle v 2 measurements of the PHENIX experiment [65] using the same data sample as this analysis showed that the tracking efficiency is independent of φ s . Thus, we normalize the probability as s n = 1, i.e. ε = 1. Due to the cyclic boundary condition in the azimuthal angle direction, symmetric elements of P (φ s,j |φ s,i ) are identical i.e. s 5 = s 3 , s 6 = s 2 , and s 7 = s 1 . The matrix P (φ s,j |φ s,i ) depends only on the order of event-planes and centrality. An example of corrections based on this iterative method at − π 24 < ∆φ < 0 for (2 < p t T < 4) ⊗ (1 < p a T < 2) GeV/c in 20%-30% central collisions is shown in Fig. 10 together with an example of the Fourier analysis method introduced in Section III D 2.

Fourier Oscillation Correction of the Event Plane-Dependence of Correlations
The second method to correct the event planedependence of the per-trigger yield is a Fourier analysis. Y (∆φ, φ s ) is offset by 1 to prevent divergences in the correction due to small values due to the ZYAM subtraction of the background. A Fourier series should be able to fit the event plane-dependence of 1 + Y (∆φ, φ s ), and the fit function to the Ψ 2 -dependent case at a given ∆φ can be written as 2a n cos n(φ s + ∆φ) (29) and similarly the Ψ 3 -dependent case can be written as where a 0 is a normalization and a 2 , a 3 , and a 4 are the azimuthal anisotropies of 1 + Y (φ s , ∆φ). In the fitting functions F Y (∆φ, φ s ), the phase shift ∆φ is necessary in 1 + Y (φ s , ∆φ) because the associated yields are at φ a − Ψ n = φ s + ∆φ (see Figure 1). With the assumption that the coefficients determined from the fits are diluted by the event plane resolutions, the effects can be corrected in a manner analogous to the single-particle azimuthal anisotropy v n as performed in Ref. [15]. For the Ψ 2 -dependent case, the correction is given as and for the Ψ 3 -dependent case it is given as The correction coefficient to 1 + Y cor is then given by the ratio F Y,cor (φ s )/F Y (φ s ), which then fixes the corrected per-trigger yield as

Efficiency
Systematic uncertainties in tracking efficiency are estimated to be approximately 10% for p T < 4 GeV/c and 13% for p T > 4 GeV/c independent of centrality [7,15,27,48].

Inclusive Per-Trigger Yields
Systematic uncertainties on the yields from the matching cut and from the v n measurements are determined by the variations on the parameters discussed below. The systematic uncertainty from the matching cut σ mat after flow subtraction is derived in similar manner as in previous publications [48] as The systematic uncertainties from v n are evaluated by taking the quadrature sum of residuals from the 1-σ un-certainties on the v n for all orders of n used in the subtraction. Formally the calculation is given by where the second Y (∆φ) refers to the yields resulting from the default set of measured v n values. The total systematic uncertainties σ in in the inclusive trigger yields are given by We studied the inclusion of a v 5 term assuming v 5 = 0.5v 4 consistent with the ATLAS measurements [33]. The results were completely consistent with the quoted uncertainties. Uncertainties due to ZYAM will be discussed later.

Event Plane-Dependent Per-Trigger Yields
In addition to systematic uncertainties considered in the inclusive per-trigger yields, systematic uncertainties due to the cos 4(Ψ 4 − Ψ 2 ) correlation strength are also taken into account before the unfolding of event plane resolution effects. The value of χ 42 is determined with v 4 {Ψ 4 } and v 4 {Ψ 2 }, and the systematic uncertainties of χ 42 are propagated from those of v 4 {Ψ 2 }. Systematic uncertainties on the yield due to χ 42 are given by The systematic uncertainties before unfolding are For the event plane-dependent per-trigger yields, the systematic uncertainty due to the impact of the finite event plane resolution on the correlations has contributions from the method and the number of iterations in the Bayesian method. The uncertainty due to the method is given by where Y itr unf is the result using the iterative Baysean method and Y fit unf is the result using the Fourier fitting method. The uncertainty due to the number of interactions using for unfolding is given by the difference between the number of iterations (Nit) for n = 5 and n = 10 These unfolding uncertainties are added in quadrature to the uncertainties before unfolding In event plane-dependent per-trigger yields, we also unfolded the upper and lower boundaries of ZYAM uncertainties propagated from statistical uncertainties from the data. The systematic uncertainties associated with ZYAM are not included into the total systematic uncertainties σ tot . How they are taken into account is described later.

IV. RESULTS AND DISCUSSION
A. Inclusive Per-Trigger Yields I: High-pT Trigger Particles We first present the highest trigger p T correlations. The jet-like correlations should be dominated by 2 → 2 scattering. Pairs of particles with ∆φ = 0, the near side, are from both particles fragmenting from a single jet. Pairs of hadrons around ∆φ = π, the away side, occur when each particle fragments from back-to-back jets. In high-p T correlations the jet momentum fraction for the associated particle is approximated by Per-trigger yields with trigger particles from 4 < p t T < 10 GeV/c paired with associated particles from 0.5 < p a T < 10 GeV/c are shown in Fig. 11. The band around zero indicates the systematic uncertainty due the ZYAM assumption. The band around the data points is the systematic uncertainty from all other sources. Systematic uncertainties from the associated tracking efficiency and matching and the ZYAM normalization are fully correlated point-to-point. The underlying event subtraction is correlated point-to-point and can affect the shape. For the highest trigger p T correlations, the dominant systematic is not the underlying event subtraction.
The near-side yield is centrality independent (Fig. 11). This is consistent with measurements of the two-particle correlation that indicated the near-side yields are not modified [20,26,29,46]. The lack of centrality dependence is also consistent with the picture that triggering on high-p T particles biases the origin of the hardscattering toward the surface of the QGP such that the leading parton loses little to no energy.
The away-side peak is evident in several p a T and centrality selections. The evolution of the away shape and yield with centrality and p a T is similar to previous measurements where the underlying event subtracted is modulated by v 2 only [28,47,48]. The away-side peak becomes sharper and more pronounced as p a T increases or the centrality selection becomes more peripheral. In more central collisions and lower p a T when the away-side structure is present, it is broader than in the highest p a T and peripheral centrality selection. The trends are consistent with a picture where the associated parton opposite the trigger loses energy and scatters in the medium. At the lowest p a T a very wide plateau-like away-side structure is GeV/c, and (p)-(t) (4 < p t T < 10) ⊗ (0.5 < p a T < 1) GeV/c. The columns represent centrality bins 0%-10% (a,f,k,p), 10%-20% (b,g,l,q), 20%-30% (c,h,m,r), 30%-40% (d,i,n,s), 40%-50% (e,j,o,t). Systematic uncertainties due to track matching and the vn are shown by blue bands around the points. Uncertainties from ZYAM are shown by the purple bands around zero yield. observed with similar shape and magnitude in all centralities. Similar low-momentum and large-angle yields have been observed in prior observables [19,[22][23][24]. Figure 12 shows the comparison between the highest p t T correlations for each centrality with the same distributions measured in p+p collisions from a previous analysis [48]. In that paper the lowest p a T bin was 0.4-1.0 GeV/c compared to 0.5-1.0 GeV/c in this analysis. Therefore, the lowest p a T bin from p+p was modified by a ∆φ-dependent correction determined from pythia 6 [76]. The correction, which has negligible uncertainties compared to those from other sources, was determined from the ratio of fits to the pythia dihadron ∆φ pertrigger yield distributions with 0.5 < p a T < 1.0 GeV/c and 0.4 < p a T < 1.0 GeV/c. Previous correlation analyses that relied on v 2 -only subtraction indicated the near-side yield was enhanced in Au+Au compared to p+p, the so-called "ridge" [46,48]. Our updated underlying event model has reduced the near-side yield as expected [31]. In fact, the yields are slightly suppressed relative to p+p. The integrated awayside yields show modification relative to p+p. Figures 12(b) and (d) show the comparisons of the pertrigger yields for the lowest p a T . The away-side shapes of the Au+Au distributions are different than p+p. The large-angle enhancement of the per-trigger yield at low associated particle momentum is qualitatively consistent with direct photon-hadron and jet-hadron correlations with fully reconstructed jets [19,[22][23][24].
To explore these features quantitatively, we calculate the ratio (I AA ) of the away-side yields in Au+Au to those in p+p. Figure 13(b) shows I AA vs. centrality when integrating the away side 0 < |∆φ − π| < π/2 for 4-10 GeV/c hadrons paired with 0.5-1.0 GeV/c hadrons. I AA is unity within uncertainties indicating that yield suppression is disfavored. Figure 13(a) shows I AA for two different angular regions of integration and p a T selections. First, for p a T from 4-10 GeV/c (high z T ) and integrating 0 < |∆φ − π| < π/4, the jet peak region, I AA is less than unity indicating the pair yields are suppressed relative to those in p+p. This is consistent with previous measurements of strong suppression of high p a T [19,27]. When integrating π/4 < |∆φ − π| < π/2 for p a T from 0.5-1.0 GeV/c (low z T ), I AA ∼ 1 within systematic uncertainties. This would indicate that the yield in Au+Au is similar to p+p. However, it is more instructive to compare the I AA for a fixed p t T , which approximately fixes the jet energy. Figure 13(a) shows that the low-z T fragments at large angles from ∆φ = π are significantly enhanced compared to the suppressed level of high-z T fragments within the jet region. Both the high-z T suppression relative to p+p and the enhanced level of low-z T fragments at large angles are consistent with a radiative energy loss model where the away-side jet traverses the medium, loses energy, and the energy gets redistributed to larger angles.

B. Inclusive Per-Trigger Yields II: Intermediate-pT Trigger Particles
Given the success of our background model at reproducing prior correlation results at high p t T , we study lower p t T correlations to attempt to measure jet-like correlations at lower momentum transfer Q 2 . Per-trigger yields with trigger particles of 1 < p t T < 2 and 2 < p t T < 4 GeV/c paired with associated particles of 0.5 GeV/c < p a T < p t T in several centrality selections are shown in Fig. 14. As in Fig. 11, the ZYAM uncertainties are shown as a band around zero while v n uncertainties are combined as the band around the data points. At these p t T the jet-like signal-to-underlying-event back- Away-side yields show the well-known suppression at high zT , most pronounced for the small angle region around the usual away-side peak center |∆φ−π| < π/4. At low zT , the large angle integration region, |∆φ − π| > π/4, show an enhancement in IAA which is significantly higher than the high zT suppressed values, and generally enhanced above unity The full away-side integration region at low zT is also higher than the suppressed level with at least 1σ significance for most centrality bins.
ground is reduced making the contribution of the v n uncertainties dominant. Because the v n uncertainties are point-to-point correlated, it is important to recognize that the yields and shape change due to that correlation. For example, if the v 2 subtracted is too large, the effect on the away side is a reduced peak and an enhanced large-angle yield. If the v 2 subtracted is too small, the away-side becomes more peaked. In the discussion that follows we only make statements that have a large systematic significance.
The away-side yield and shape varies with both p a T and centrality. In these p T selections and in the most central collisions, the away-side seems to completely disappear. If our background model represents all nonjet correlations, the disappearance is presumably due to jet quenching. Compared to v 2 -only subtraction [48], the very large displaced away-side peaks are reduced primarily due to the subtraction of v 3 in the underlying event [31]. Both the flat away-side and the near-side peak shape, seems relatively centrality independent.
To better assess the systematic significance of the correlation features at the lower p a T selections, Figure 15 shows a further breakdown of the uncertainty contributions from different orders of the v n subtraction. The lines and bands break out the 1σ and 2σ variation of each of the v n components independently. Again, these The columns represent centrality bins 0%-10% (a,f,k,p,u), 10%-20% (b,g,l,q,v), 20%-30% (c,h,m,r,w), 30%-40% (d,i,n,s,x), 40%-50% (e,j,o,t,y). Systematic uncertainties are shown by (blue) bands around the points. Uncertainties from ZYAM are shown by (purple) bands around zero yield.
uncertainties are point-to-point correlated where a reduction at ∆φ = π leads to an increase at angles away from π. The v 4 dominates the uncertainty in the awayside shape. The location of any double-peak structure in the away-side is also strongly dependent on the underlying event background subtraction. The peak in the uncertainties from v 2 , v 3 and v 4 are all at slightly different places. Therefore, a robust and well-motivated background model is necessary to extract detailed shape information at these p a T . With our quoted systematic uncertainties, we cannot distinguish between a single broad away-side peak structure or a double-hump structure at the 2σ level.

C. Event Plane-Dependent Correlations I: Ψ2 Dependence
An important goal of jet quenching studies has been to determine the density and path-length dependence of energy loss [15]. For example, perturbative models of radiative jet quenching and strongly coupled jet quenching models predict a different path-length dependence for the quenching [77]. Varying the path length by selecting azimuthal orientations relative to the second-order event plane has been explored for single-particle or single-jet observables at high-p T [11,78]. Potentially more differential information can be obtained from two-particle observables coupled with the event-plane. For example, Fig. 16 illustrates the trigger (Jet A) being emitted to one side of the in-plane direction and the away-side jet (Parton B) radiating two gluons (Gluon 1) and (Gluon 2). Figure 17 shows Ψ 2 -dependent per-trigger yields of trigger particles from 2 < p t T < 4 GeV/c with associated particles from 1 < p a T < 2 GeV/c in 10%-20% and 30%-40% central collisions. In Figs. 17 (b) and (d), the trigger is selected to be either in-plane 0 < |φ s | < π/8 or out-of-plane 3π/8 < |φ s | < π/2, respectively. Similar to the inclusive per-trigger yields, we observe a broad away-side structure in both cases. The use of a common ZYAM point results in a slight over-subtraction in the out-of-plane bins. The over subtraction can be corrected by determining a ZYAM point for each φ s selection. This would result in moving all yield points up and does not affect the discussion of the shape that follows.
In Figs. 17 (a) and (c), we chose the trigger to have a particular sign of φ s . That is, we choose −π/8 < φ s < 0 and −π/2 < φ s < −3π/8 for the in-plane and out-ofplane, respectively. Choosing the sign of φ s to be neg- GeV/c. The columns represent centrality bins 30%-40% (a,c,e) and 40%-50% (b,d,f). The lines and bands further break down of the uncertainty contributions from each different order of the vn subtraction. The systematic uncertainties are point-to-point correlated. If the yield at ∆φ = π is reduced, the yield off of ∆φ = π is increased. ative results in always choosing the trigger to be "below" the event plane, if the event plane is the horizontal. When sign-selecting φ s , an asymmetry around ∆φ ∼ 0 and ∆φ ∼ π is observed. Such an asymmetry does not exist when choosing both signs of φ s . In the in-plane trigger case, there is a preference for the associated particle to be emitted toward the in-plane direction i.e. the thinner side of the overlap region. Referring back to Fig. 16, our data suggests that Gluon 1 is more likely to be measured than Gluon 2. Other Ψ 2 dependent per-trigger yields with a selection of trigger azimuthal angle relative to event-planes φ t − Ψ n < 0 for different collision centralities and p t,a T selections are shown in the Appendix. The integrated per-trigger yields from sign-selected trigger angles with respect to Ψ 2 are shown in Fig. 18 as a function of associated particle angle with respect to Ψ 2 . Figure 18 shows data for all centralities and all four orientations with respect to the event plane. We note that the use of a single ZYAM level for all event plane bins can result in negative yields for certain ∆φ.
In Figs. 18 (a)-(e), the yields have been integrated for the near side |∆φ| =< π/4. For all trigger/associated combinations there is weak to no dependence of the yield on event plane orientation. This trend persists for all centrality selections.
In Figs. 18 (f)-(j), the yields have been integrated for the away side |∆φ − π| < π/4. In the most central collisions, the largest yield is out of plane. This trend is opposite in the most peripheral collision selection where the largest yield is in plane. In 10%-40% central collisions similar yields are observed in-plane versus out of plane. For the event plane selections in the range −1.2 < φ a − Ψ 2 < −0.5, no significant yields are observed. Because each of the trends with event plane is not affected by the level of the ZYAM uncertainty, we conclude there is a trend with centrality that out-ofplane yield is reduced from central to peripheral collisions whereas in-plane yield increases from central to peripheral collisions. The trends with event plane and centrality are observed in both the highest and lowest trigger p T .
Correlations selecting the trigger within a certain azimuthal angles from the Ψ 3 plane are shown in Figure 19. The triangular shape of the third Fourier component restricts the range: 0 < Ψ 3 ≤ 2π/3. The most out-of-plane direction, which we will refer to as "out-of-plane", is at π/3 radians relative to Ψ 3 . In Figs. 19 (b) and (d), we do not select the sign for φ s . Similar to the inclusive distributions, there is a broad away-side structure that is has a small yield. In Figs. 19 (a) and (c), we select for φ s < 0. No discernible asymmetry is observed. It is possible that unfolding with the smaller Ψ 3 event plane resolution could obscure any effect. Other Ψ 3 dependent per-trigger yields with a selection of trigger azimuthal angle relative to event-planes φ t − Ψ n < 0 for different collision centralities and p t,a T selections are shown in the Appendix.
Similar to the Ψ 2 -dependent correlations, we also integrate the per-trigger yields. This is shown in Fig. 20 for each centrality and associated particle azimuthal angle with respect to Ψ 3 . Figures 20 (a)-(e) show the near side integral. Figures 20 (f)-(j) show the away side integral. In all cases, no event plane-dependent or centrality dependent trends can be observed within uncertainties.

V. SUMMARY
In summary, we reported the two-particle azimuthal dihadron correlation measurements at |∆η| < 0.7 in Au+Au collisions at √ s N N = 200 GeV after subtraction of an underlying event background model. The model includes modulations from higher-order flow coefficients v n (n = 2, 3, 4) that assumes only the expected correlation of 2nd and 4th order event planes. We tested this two-source model by studying high-p T (> 4 GeV/c) triggers. We observe suppression of highz T jet fragments as well as enchantment of low-z T jet fragments off the away-side jet axis. These results are consistent with previous dihadron and γ-hadron correlations and jet analyses [6][7][8][9][10][11].
At lower trigger p T 2 < p t T < 4 GeV/c, the near-side distribution is not enhanced compared to p+p, which traditionally is associated with the ridge. When a significant away-side yield exists, the double-hump structure that had been observed when subtracting a v 2 -only underlying event is significantly reduced. Given our model assumptions and the systematic uncertainties on v n , we cannot precisely determine if the away-side distribution is a single broadened peak or has further structure that may peak away from ∆φ = π.