Constraining the initial conditions and temperature dependent transport with three-particle correlations in Au+Au collisions

We present three-particle mixed-harmonic correlations $\la \cos (m\phi_a + n\phi_b - (m+n) \phi_c)\ra$ for harmonics $m,n=1-3$ for charged particles in $\sqrt{s_{NN}}=$200 GeV Au+Au collisions at RHIC. These measurements provide information on the three-dimensional structure of the initial collision zone and are important for constraining models of a subsequent low-viscosity quark-gluon plasma expansion phase. We investigate correlations between the first, second and third harmonics predicted as a consequence of fluctuations in the initial state. The dependence of the correlations on the pseudorapidity separation between particles show hints of a breaking of longitudinal invariance. We compare our results to a number of state-of-the art hydrodynamic calculations with different initial states and temperature dependent viscosities. These measurements provide important steps towards constraining the temperature dependent transport and the longitudinal structure of the initial state at RHIC.

We present three-particle mixed-harmonic correlations cos(mφa + nφ b − (m + n)φc) for harmonics m, n = 1 − 3 for charged particles in √ sNN =200 GeV Au+Au collisions at RHIC. These measurements provide information on the three-dimensional structure of the initial collision zone and are important for constraining models of a subsequent low-viscosity quark-gluon plasma expansion phase. We investigate correlations between the first, second and third harmonics predicted as a consequence of fluctuations in the initial state. The dependence of the correlations on the pseudorapidity separation between particles show hints of a breaking of longitudinal invariance. We compare our results to a number of state-of-the art hydrodynamic calculations with different initial states and temperature dependent viscosities. These measurements provide important steps towards constraining the temperature dependent transport and the longitudinal structure of the initial state at RHIC.
Introduction : Matter as hot and dense as the early universe microseconds after the Big Bang can be created by colliding heavy nuclei at high energies. At these temperatures, baryons and mesons melt to form a quark gluon plasma (QGP) [1][2][3][4]. Data from the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory and the Large Hadron Collider (LHC) at CERN have been arguably used to show that the QGP at these temperatures is a nearly perfect fluid with a shear viscosityto-entropy density ratio (η/s) smaller than any other fluid known in nature [5][6][7][8][9][10][11][12][13]. Theoretical calculations suggest that like many other fluids, the QGP viscosity should have a dependence on temperature with a minimum at the QGP-to-hadron transition temperature [14][15][16]. The determination of the temperature dependence of these transport properties is an open problem of fundamental importance in the study of the emerging properties of QCD matter.
Over the past years the harmonic decomposition of two-particle azimuthal correlations v 2 n {2} = cos n(φ a − φ b ) (where φ a,b are azimuthal angles of particle momenta) [12,[17][18][19][20] have already helped shed light on these topics. Hydrodynamic models with different initial conditions and transport parameters have been compared to measurements at RHIC and LHC to constrain the fluid-like property of the medium [21]. Given their large number of parameters, measurements of multiple observables over a wide energy range have been found to be essential for constraining such models [22][23][24]. So far however, the temperature dependence of transport parameters like the bulk and shear viscosity are not well constrained by the existing data.
In this letter, we report on the measurement of threeparticle correlations that provide unique ways to constrain the fluid-like properties of the QGP. These new measurements at RHIC extend beyond the conventional two-particle correlations; they help elucidate the three dimensional structure of the initial state, probe the nonlinear hydrodynamic response of the medium, and will help constrain the temperature dependence of the transport parameters. We measure three-particle azimuthal correlations using the observables [25] C m,n,m+n = cos(mφ where the inner average is taken over all sets of unique triplets, and the outer average is taken over all events weighted by the number of triplets in each event. The subscripts "m, n" in C m,n,m+n refer to the harmonic number while the subscripts "a, b, c" in φ refer to the indices of the particles. We report on the centrality dependence of C m,n,m+n with combinations of harmonics (m, n) = (1, 1), (1,2), (2,2), (2,3), (2,4) and (3,3) for inclusive charged particles in Au+Au collisions at √ s N N = 200 GeV. In a longer companion paper [26] we present our measurements at lower energies ( √ s N N = 62.4-7.7 GeV). The C m,n,m+n are related to eventplane correlations like those measured in Pb+Pb collisions at 2.76 TeV [27][28][29]. If v n and Ψ n denote [30] anisotropic flow coefficients and their associated event planes [31], for m, n > 1, C m,n,m+n can be approximated as v m v n v m+n cos(mΨ m + nΨ n − (m + n)Ψ m+n ) . Such flow based interpretation is not likely to be applicable in case of m, n = 1 for which a strong charge dependence has been observed [32][33][34] and the effects of global momentum conservation may be important [35,36]. Measurements of C m,n,m+n provide unique information about the geometry of the collision overlap region and its fluctuations. Reference [37] proposed that measurements of C 1,2,3 could detect event-by-event correlations of the first, second and third harmonic anisotropies. Although it is sometimes assumed that the axis of the third harmonic is random, Monte-Carlo Glauber simulations show correlations between the first, second, and third harmonic planes. Figure 1 (left) shows the case when a single nucleon (shown by a red dot) at the edge of a colliding nucleus fluctuates outward and impinges on the other nucleus creating a region of increased energy density. This specific in-plane fluctuation generates v 1 , which reduces v 2 and increases v 3 [38]. A similar fluctuation occurring in the out-of-plane direction is illustrated in the right panel of Fig. 1. Such correlations, if observed in terms of C 1,2,3 , will for the first time, demonstrate the presence of a v 1 driven component of v 3 arising due to initial geometry.
The fluctuation illustrated in Fig. 1 (left) when the nucleon at the edge of one nucleus impinges on the center of the other nucleus, it is similar to a central p+Au collision. In p+Au collisions, the maximum of the multiplicity distribution shifts in pseudorapidity η towards the Au going direction. For this reason, one expects that the harmonic planes can point in different directions for positive or negative η. Similar effects have been investigated in models and discussed in terms of torqued fireballs [39], twists [40], or reaction-plane decorrelations [41]. Studying the ∆η dependence of C 1,2,3 should reveal these effects if they exist, and provide new insights on the three dimensional structure of the initial state.
In general, if a medium is fully describable by hydrodynamics, nonlinear couplings between harmonics are expected to change the sign of C m,n,m+n relative to what would be expected based on the initial state eccentricities ε n [42] and participant planes Φ n [25,37,[43][44][45][46][47][48][49][50]. Observables sensitive to nonlinear hydrodynamic response are ideal probes of viscosity. Since higher harmonics are more strongly dampened by viscosity, the nonlinear coupling increases correlations of v n with other lower harmonic eccentricities ε m<n , and thereby with v m<n . In this way, C m,n,m+n becomes more sensitive to η/s as previously demonstrated by phenomenological studies at LHC energies [25,43,45,51]. Correlations of event planes and flow harmonics measured by the ATLAS and AL-ICE collaborations for m, n ≥ 2 [19,28,29] have been compared to hydrodynamic simulations to constrain the temperature dependence of viscosity η/s (T ) [51]. However since LHC measurements are sensitive to the η/s at higher temperatures, full constraint on η/s (T ) is better achieved with measurements of observables like C m,n,m+n at RHIC [11,[51][52][53].
In this work we report the three-particle correlations directly instead of event-plane correlations. Expressing three-particle correlations as event plane correlations relies on factorization, i.e., approximations like C m,n,m+n = v m v n v m+n cos(mΨ m + nΨ n − (m + n)Ψ m+n ) = v m v n v m+n cos(mΨ m + nΨ n − (m + n)Ψ m+n ) , that can complicate data-model comparison. We therefore, directly compare C m,n,m+n to theoretical predictions. Another advantage of three-particle correlations is that the measurements are well defined even without assuming the flow coefficients and harmonic planes dominate the correlation. Other effects besides reaction plane correlations, particularly important for m, n = 1, can be present in C m,n,m+n and the correctness and completeness of a model needs to be judged through direct comparison to the data. Also, when the correlations are dominated by reaction plane correlations, C m,n,m+n corresponds to a well-defined limit (the low-resolution limit) [54] of the measurement, which again, makes for a more direct comparison to theory. A more practical advantage is as follows: unlike LHC, since v 2 n {2} for n = 1 − 6 is not always a large positive quantity at RHIC, it is not always feasible to divide C m,n,m+n by v 2 n {2} to express it purely as an event plane correlation without losing experimental significance. The magnitude of v 2 6 {2} is negligible at RHIC, v 2 5 {2} measurements suffer from large systematics, and v 2 1 {2} < 0 except for central events at √ s N N = 200 GeV [26].
Experiment and Analysis : We present measurements of C m,n,m+n in 200 GeV Au+Au collisions with data collected in the year 2011 by the STAR detector [55] at RHIC. We detect charged particles within the range |η| < 1 and for transverse momentum of p T > 0.2 GeV/c using the STAR Time Projection Chamber [56] situated inside a 0.5 Tesla solenoidal magnetic field. We use trackby-track weights [57,58] to account for imperfections in the detector acceptance and momentum dependence of the detector efficiency. We correct the two-track acceptance artifacts which arise due to track-merging effects by measuring the |∆η ab | = |η a − η b |, |∆η ac | = |η a − η c |, and |∆η bc | = |η b − η c | dependence of C m,n,m+n and algebraically correcting the integrated value of C m,n,m+n for the missing pairs apparent at ∆η ≈ 0. Note that, throughout this paper, the subscripts "m, n with comma" in C m,n,m+n refer to the harmonic number while the subscripts "ab without comma" for the |∆η ab | = |η a − η b | refer to the indices of the particles. We estimate systematic uncertainties by comparing data from different time periods, from different years with different tracking algorithms, by comparing different efficiency estimates, by varying the z-vertex position of the collision, and by varying track selection criteria. We also include estimates of the effect of short-range HBT and Coulomb correlations in the systematic uncertainties based on the shape of the ∆η dependence. For such quantifications we fit the ∆η dependence of C m,n,m+n with the combination of a shortrange and a long-range Gaussian distributions as described in Ref [38,59]. Finally, in order to quantify other nonflow effects such correlations due to mini-jets, fragmentation, decay etc. we compare our data to HIJING (Version 1.383) calculations [60]. For each of our centrality intervals (0 − 5%, 5 − 10%, 10 − 20%, ..., 70 − 80%), we use a Monte Carlo Glauber model [61,62] to estimate the average number of participating nucleons N part for plotting our results [63].
to the symmetry of collision geometry.
As mentioned before, since C 1,2,3 involves the first order harmonic it may have contributions from nonflow correlations such as global momentum conservation [35]. However, such contributions have been argued to be independent of ∆η in leading order [35,64,65]. One, therefore, can not explain the strong variation of C 1,2,3 with |∆η ac | even up to 2, which is strongest in the mid-central events, to be only as an artifact of momentum conservation.
The HIJING model comparisons shown in Fig. 2 demonstrate that nonflow contributions due to mini-jets can not explain data. On the other hand the AMPT model [66] calculations from Ref. [67] that involves momentum conservation, mini-jets, as well as collectivity due to multiphase transport, and three-dimensional initial state seem to provide a better description of the ∆η dependence of C 1,2,3 above ∆η > 0.5; at smaller ∆η < 0.5 AMPT under predicts the data.
In Fig. 2 (b) we present the ∆η dependence of C 2,2,4 . We find much weaker ∆η dependence for C 2,2,4 than for C 1,2,3 ; while C 1,2,3 changes sign, C 2,2,4 only varies by 20% over the range of our measurements. This is not surprising since the second harmonic event plane dominates C 2,2,4 . The dependence of C 2,2,4 is also stronger for |∆η ac | than it is for |∆η ab |. Once again, the HIJING predictions (not shown in this figure) are much smaller and consistent with zero. The AMPT predictions from Ref [67] do a very good job in describing the magnitude of the correlation, it however, seem to slightly under predict the slope of the ∆η dependence.
We find that all the correlators exhibit a significant ∆η dependence except C 2,2,4 and C 2,3,5 which vary by only 20% [26]. The variation of C m,n,m+n with ∆η makes it difficult to compare the data to models that assume a longitudinally invariant two-dimensional (boost invariant) initial geometry. Until those simplifying assumptions are relaxed, C 2,2,4 and C 2,3,5 having the smallest relative variation on ∆η provide the best opportunity for comparison of ∆η-integrated quantities with hydrodynamic models.
In Fig. 3 we show centrality dependence of ∆ηintegrated C m,n,m+n .
We multiply the quantity C m,n,m+n by N 2 part to account for the natural dilution of correlations expected from superpositions of independent sources. We find that HIJING model predicts a magnitude of three-particle correlations that is consistent with zero for all harmonics. We also estimate the expectations for C m,n,m+n ≈ ε m ε n ε m+n cos(mΦ m + nΦ n − (m + n)Φ m+n ) from purely initial state geometry using a Monte-Carlo Glauber model [68].
We find that the Glauber model predicts negative values for all combinations of C m,n,m+n [69]. Since only a fraction of the initial state geometry is converted to final state anisotropy, i.e., v n 0.1 × ε n [44], one therefore expects v m v n v m+n cos(mΨ m + nΨ n − (m + n)Ψ m+n ) 10 −3 × ε m ε n ε m+n cos(mΦ m + nΦ n − (m + n)Φ m+n ) , we therefore scale the Glauber model calculations by factors of ∼ 10 −3 − 10 −4 to make a consistent data to model comparison [44].
We compare our results with three different boostinvariant hydrodynamic model calculations that have been constrained by the global data on azimuthal correlations available so far at RHIC and the LHC. The models include : 1) 2+1 dimensional hydrodynamic simulations with η/s = 1/4π with MC-Glauber initial conditions by Teaney and Yan [37,45], 2) hydrodynamic simulations MUSIC with boost invariant IP-Glasma initial conditions [70,71] that include a constant η/s = 0.06 and a temperature dependent bulk viscosity ζ/s (T ) [72] and UrQMD afterburner [73], 3) the perturbative-QCD+saturation+hydro based "EKRT" model [51] that uses two different parameterizations of the viscosity with constant η/s = 0.2 and temperature dependent η/s (T ) with a minimum of (η/s (T )) min = 1.5/4π at a corresponding transition temperature between a QGP and hadronic phase of T c = 150 MeV and 4) viscous hydrodynamic model v-USPhydro [74,75] with event-byevent TRENTO initial conditions [76] tuned to IP-Glasma [70], that uses η/s = 0.05, a freeze-out temperature of T F O = 150 MeV [77] and the most recent 2+1 flavors equation of state from the Wuppertal Budapest collaboration [78] combined to all known hadronic resonances from the PDG16+ [79].
Correlators involving the first order harmonic C 1,1,2 and C 1,2,3 are shown in Fig. 3 (a) and (b). In Fig. 3 (a) we compare results to the hydrodynamic predictions by Teaney and Yan [37,45]. We note that since finite multiplicity effects, such as global momentum conservation, are not included in these calculations, comparisons presented for C 1,1,2 and C 1,2,3 are not intended for the purpose of constraining transport parameters.
Any dipole anisotropy with respect to the second order harmonic plane will be exhibited in the correlator C 1,1,2 = cos(φ a +φ b −2φ c ) . The negative value of C 1,1,2 observed in Fig. 3 (a) indicates that the dipole anisotropy arising at mid-rapidity is dominantly out-of-plane as predicted by the theoretical calculations in Ref. [37] and initial state geometry. It may also indicate a significant contribution from momentum conservation [64,65]. For the correlator C 1,1,2 , it was explicitly shown that a combination of flow and momentum conservation gives rise to a negative contribution (∼ −v 2 /N , N being the multiplicity) [64,65]. The models do not include such effects; therefore it is not surprising that they significantly under predict the data.
The centrality dependence of C 1,2,3 is shown in Fig. 3  (b). We see a nonzero correlation consistent with the illustrations in Fig. 1. The large positive values of C 1,2,3 in mid-central events are indicative [80] of the first harmonic anisotropy correlated with the triangularity as was first predicted in Ref. [37]. In the model, the hydrodynamic response of the medium changes both the sign and the centrality dependence and provides very good agreement with data for C 1,2,3 over a wide range of N part except for the most central collisions. Interestingly in the most central collisions, the measurements of both C 1,1,2 and C 1,2,3 are nonzero and negative while the models predict nearly zero values for these correlators which might need further investigation [81].
We next report the measurement of the correlators C 2,2,4 and C 2,3,5 in Fig. 3 (c)-(d). The correlator C 2,2,4 ≈ v 2 2 v 4 cos(4(Ψ 2 − Ψ 4 )) measures the correlation between the second and the fourth order harmonics and the corresponding event planes. While the Glauber model results for the initial state are negative, both C 2,2,4 and C 2,3,5 exhibit strong positive values. This is consistent with the linear and nonlinear hydrodynamic response of the medium created at RHIC, in which the higher flow harmonics like v 4 is driven by both ε 4 and ε 2 , as predicted by several theoretical calculations [25,43,[45][46][47]. This result is also qualitatively consistent with the measurement by the ATLAS collaboration at LHC [19,28].
The quantitive difference between the models and the measurement at RHIC is an important observation of the current study. In Fig. 3 (c), we observe that the hydrodynamic predictions by Teaney and Yan using constant η/s significantly underestimate C 2,2,4 . The predictions using EKRT with a temperature dependent η/s are much closer to the data; the same using constant η/s under predict data by about 20%. A similar trend is also observed for C 2,3,5 shown in Fig. 3  trends of the centrality dependence, they all significantly underestimate the magnitude of C 2,3,5 . Such discrepancy for EKRT has been argued [82] to be related to large off-equilibrium correlations which depend on the details of the parameterization η/s (T ). The current data will therefore provide important constraints for the transport parameters involved in the hydrodynamic modeling at RHIC energies. In Fig. 3 (e)-(f) we present the centrality dependence of C 2,4,6 and C 3,3,6 . Once again the positive values for C 2,4,6 and C 3,3,6 , in contrast to the Glauber prediction of negative values for the initial state, indicate the importance of the nonlinear hydrodynamic response. The EKRT predictions are not available for these correlators, it will be interesting to see if such calculations can describe the data in future.
We revisit the centrality dependence of higher order correlators (n > 2) in Fig. 4. Here, we compare the data with most recent hydrodynamic model calculations. The IP-Glasma + MUSIC simulations with constant η/s, tuned to global data on v n s, qualitatively reproduce the trend; however they under predict the magnitude of the correlation. The IP-Glasma + MUSIC + UrQMD simulations, that include additional hadronic rescatterings, seems to be much closer to the data. This is indicative of the fact that a large fraction of the mixed-harmonic correlation is developed in the hadronic phase below a temperature of T = 165 MeV. The addition of hadronic transport effectively increases the viscosity at lower temperature (T < 165 MeV) [72]. This indicates that current data can constrain the temperature dependent transport at RHIC energies. In Fig. 4 our data is also compared to the TRENTO+v-USPhydro model calculations. Although this model does not include hadronic transport, as discussed in Ref [77], it effectively introduces a different viscous effect by choosing a lower freeze-out temperature T F O = 150 MeV, additional resonances and a different equation of state (speed of sound), as compared to IP-Glasma + MUSIC + UrQMD simulations. A reasonable description of C 2,3,5 , C 2,4,6 and C 3,3,6 is obtained from the TRENTO+v-USPhydro model. In the case of C 2,2,4 the data are 20% higher, which will provide further constraints for the TRENTO+v-USPhydro model [79]. It will be also interesting to see other hydro calculations by using the most recent equation of state like TRENTO+v-USPhydro model.
After the appearance of this preprint, an extensive study using the AMPT model was shown to provide a good description of both the ∆η and the centrality dependence of C m,n,m+n in Ref. [67]. Such data-model comparisons demonstrate that the longitudinal structure of the initial state, global momentum conservation and multiphase transport can capture the underlying dynamics that drives anisotropic flow and mixed-harmonic corre-lations [67].

Summary :
We presented the first measurements of the charge inclusive three-particle azimuthal correlations C m,n,m+n = cos(mφ a + nφ b − (m + n)φ c ) as a function of centrality, relative pseudorapidity and harmonic numbers m, n in √ s N N = 200 GeV Au+Au collisions. These measurements, provide additional information about the initial geometry, the nonlinear hydrodynamic response of the medium and provide good promise to constrain temperature dependence of η/s. The centrality dependence of C 1,2,3 for the first time reveals a possible coupling between directed, elliptic, and triangular harmonic flow, which arises from fluctuations in the initial geometry. The strong ∆η dependence of C 1,2,3 suggests a breaking of longitudinal invariance at odds with the assumptions in many boost invariant models.
While variations of C 1,2,3 with ∆η are large, C 2,2,4 and C 2,3,5 varies by only 20% between ∆η = 0 and 2 making them most suitable for comparison to boost-invariant hydrodynamic simulations. We therefore, compared our measurements of the centrality dependence of C m,n,m+n with a number of boost-invariant hydrodynamic models that are constrained by global data. Such comparisons indicate that three-particle correlations can provide important constraints on fluid-dynamical modeling, in particular the temperature dependent transport at RHIC.