STAR Collaboration

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) Collins and Perry (1975); Chin (1978); Kapusta (1979); Anishetty et al. (1980). 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 viscosity-to-entropy density ratio ($\eta$/s) smaller than any other fluid known in nature Arsene et al. (2005); Back et al. (2005); Adams et al. (2005a); Adcox et al. (2005); Muller (2007); Zajc (2008); Gale et al. (2013a); Chatrchyan et al. (2014); Abelev et al. (2015). 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 Prakash et al. (1993); Arnold et al. (2003); Csernai et al. (2006). 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_{n}^{2}\{2\}=\langle\cos n(\phi_{a}-\phi_{b})\rangle$ (where $\phi_{a,b}$ are azimuthal angles of particle momenta) Adams et al. (2005b); Adare et al. (2011); Aad et al. (2014a); Abelev et al. (2014); Chatrchyan et al. (2014) 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 Gale et al. (2013b). 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 Novak et al. (2014); Pratt et al. (2015); Bernhard et al. (2015). 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 three-particle 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 Bhalerao et al. (2013)

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 $\phi$ 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 $\sqrt{s_{NN}}=$ 200 GeV. In a longer companion paper Adamczyk et al. (2017) we present our measurements at lower energies ($\sqrt{s_{NN}}=$ 62.4-7.7 GeV). The $C_{m,n,m+n}$ are related to event-plane correlations like those measured in Pb+Pb collisions at 2.76 TeV Aad et al. (2014b); Jia (2014); Adam et al. (2016). If $v_{n}$ and $\Psi_{n}$ denote ¹¹1 $$v_{n}\,e^{in\Psi_{n}}=\frac{\int p_{T}\,dp_{T}\,d\phi\,e^{in\phi}\,\frac{dN_{\mathrm{ch}}}{d\eta\,p_{T}dp_{T}\,d\phi}}{\int p_{T}dp_{T}\,d\phi\,\frac{dN_{\mathrm{ch}}}{d\eta\,p_{T}dp_{T}\,d\phi}}.$$ (2) where $\frac{dN_{\mathrm{ch}}}{d\eta\,p_{T}dp_{T}\,d\phi}$ is the single particle distribution. anisotropic flow coefficients and their associated event planes Voloshin et al. (2008), for $m,n>1$, $C_{m,n,m+n}$ can be approximated as ${\langle v_{m}v_{n}v_{m+n}\cos(m\Psi_{m}+n\Psi_{n}-(m+n)\Psi_{m+n})\rangle}$. 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 Abelev et al. (2009a, 2010); Adamczyk et al. (2013) and the effects of global momentum conservation may be important Borghini et al. (2000); Jia (2012).

Measurements of $C_{m,n,m+n}$ provide unique information about the geometry of the collision overlap region and its fluctuations. Reference  Teaney and Yan (2011) 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}$ Adamczyk et al. (2016). 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 $\eta$ towards the Au going direction. For this reason, one expects that the harmonic planes can point in different directions for positive or negative $\eta$. Similar effects have been investigated in models and discussed in terms of torqued fireballs Bozek et al. (2011), twists Jia and Huo (2014), or reaction-plane decorrelations Pang et al. (2016). Studying the $\Delta\eta$ 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 $\varepsilon_{n}$ ²²2 $$\varepsilon_{n}\,e^{in\Phi_{n}}=-\frac{\int r\,dr\,d\phi\,r^{n}e^{in\phi}\,E(r,\phi)}{\int r\,dr\,d\phi\,r^{n}\,E(r,\phi)}$$ (3) where $E(r,\phi)$ is the distribution of initial energy density. and participant planes $\Phi_{n}$  Teaney and Yan (2011); Qiu and Heinz (2012); Teaney and Yan (2012, 2014); Bhalerao et al. (2013); Yan and Ollitrault (2015); Qian et al. (2016); Qian and Heinz (2016); McDonald et al. (2016); Noronha-Hostler et al. (2016). 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 $\varepsilon_{m<n}$, and thereby with $v_{m<n}$. In this way, $C_{m,n,m+n}$ becomes more sensitive to $\eta/s$ as previously demonstrated by phenomenological studies at LHC energies Qiu and Heinz (2012); Teaney and Yan (2014); Bhalerao et al. (2013); Niemi et al. (2016). Correlations of event planes and flow harmonics measured by the ATLAS and ALICE collaborations for $m,n\!\geq\!2$ Aad et al. (2014a); Jia (2014); Adam et al. (2016) have been compared to hydrodynamic simulations to constrain the temperature dependence of viscosity $\eta/s~(T)$ Niemi et al. (2016). However since LHC measurements are sensitive to the $\eta/s$ at higher temperatures, full constraint on $\eta/s~(T)$ is better achieved with measurements of observables like $C_{m,n,m+n}$ at RHIC Niemi et al. (2011); Gale et al. (2013a); Niemi et al. (2016); Denicol et al. (2016).

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}=\langle v_{m}v_{n}v_{m+n}\cos(m\Psi_{m}+n\Psi_{n}-(m+n)\Psi_{m+n})\rangle=$ $\langle v_{m}\rangle\langle v_{n}\rangle\langle v_{m+n}\rangle\langle\cos(m\Psi_{m}+n\Psi_{n}-(m+n)\Psi_{m+n})\rangle$, 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)  Luzum and Ollitrault (2013) of the measurement, which again, makes for a more direct comparison to theory. A more practical advantage is as follows: unlike LHC, since $v_{n}^{2}\{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 $\sqrt{v_{n}^{2}\{2\}}$ to express it purely as an event plane correlation without losing experimental significance. The magnitude of $v_{6}^{2}\{2\}$ is negligible at RHIC, $v_{5}^{2}\{2\}$ measurements suffer from large systematics, and $v_{1}^{2}\{2\}<0$ except for central events at $\sqrt{s_{NN}}=200$ GeV Adamczyk et al. (2017).

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 Ackermann et al. (2003) at RHIC. We detect charged particles within the range $|\eta|\!<\!1$ and for transverse momentum of $p_{T}\!>\!0.2$ GeV/$c$ using the STAR Time Projection Chamber Anderson et al. (2003) situated inside a 0.5 Tesla solenoidal magnetic field. We use track-by-track weights Bilandzic et al. (2011, 2014) 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 $|\Delta\eta_{ab}|=|\eta_{a}-\eta_{b}|$, $|\Delta\eta_{ac}|=|\eta_{a}-\eta_{c}|$, and $|\Delta\eta_{bc}|=|\eta_{b}-\eta_{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 $\Delta\eta\approx 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 $|\Delta\eta_{ab}|=|\eta_{a}-\eta_{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 $\Delta\eta$ dependence. For such quantifications we fit the $\Delta\eta$ dependence of $C_{m,n,m+n}$ with the combination of a short-range and a long-range Gaussian distributions as described in Ref Tribedy (2017); Adamczyk et al. (2016). 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 Wang and Gyulassy (1991). For each of our centrality intervals ($0-5\%,5-10\%,10-20\%,...,70-80\%$), we use a Monte Carlo Glauber model Abelev et al. (2009b); Miller et al. (2007) to estimate the average number of participating nucleons $N_{\rm part}$ for plotting our results ³³3See Ref. Abelev et al. (2009b) for details like centrality resolution, values of impact parameter, $N_{\rm part}$ etc..

Results : Figure 2 (a,b) shows the $\Delta\eta$ dependence of $C_{1,2,3}=\langle\cos(\phi_{a}+2\phi_{b}-3\phi_{c})\rangle$ and $C_{2,2,4}=\langle\cos(2\phi_{a}+2\phi_{b}-4\phi_{c})\rangle$. The $\Delta\eta$ dependence of $C_{1,1,2}=\langle\cos(\phi_{a}+\phi_{b}-2\phi_{c})\rangle$ was presented previously Abelev et al. (2010); Adamczyk et al. (2013) and other harmonic combinations will be presented in Ref. Adamczyk et al. (2017). The top panel of Fig. 2 shows $C_{1,2,3}$ as a function of $|\Delta\eta_{ab}|$ and $|\Delta\eta_{ac}|$. We observe a strong $|\Delta\eta_{ac}|$ dependence but a weak $|\Delta\eta_{ab}|$ dependence. The observation for $|\Delta\eta_{bc}|$ is similar to $|\Delta\eta_{ab}|$, so we omit it from the figure for clarity. For $|\Delta\eta_{ac}|\approx 0$, $C_{1,2,3}$ is positive, but as $|\Delta\eta_{ac}|$ increases, $C_{1,2,3}$ decreases and becomes negative. We study the centrality dependence of this effect and find that $C_{1,2,3}$ has the strongest dependence on $|\Delta\eta_{ac}|$ in mid-central events (20-30$\%$); in central (0-5$\%$) and peripheral events (70-80$\%$), $C_{1,2,3}$ shows weaker dependence on $|\Delta\eta_{ac}|$ (see Ref. Adamczyk et al. (2017)). This is consistent with expectations of the breaking of longitudinal invariance through forward-backward rapidity dependence as previously discussed. No such dependence is observed for $|\Delta\eta_{ab}|$ since although the third harmonic plane may rotate significantly in the forward and backward directions, the second harmonic plane should remain invariant due 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 Borghini et al. (2000). However, such contributions have been argued to be independent of $\Delta\eta$ in leading order Borghini et al. (2000); Pratt (2010); Bzdak et al. (2011). One, therefore, can not explain the strong variation of $C_{1,2,3}$ with $|\Delta\eta_{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 Lin et al. (2005) calculations from Ref. Sun and Ko (2017) 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 $\Delta\eta$ dependence of $C_{1,2,3}$ above $\Delta\eta>0.5$; at smaller $\Delta\eta<0.5$ AMPT under predicts the data.

In Fig. 2 (b) we present the $\Delta\eta$ dependence of $C_{2,2,4}$. We find much weaker $\Delta\eta$ 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 $|\Delta\eta_{ac}|$ than it is for $|\Delta\eta_{ab}|$. Once again, the HIJING predictions (not shown in this figure) are much smaller and consistent with zero. The AMPT predictions from Ref Sun and Ko (2017) do a very good job in describing the magnitude of the correlation, it however, seem to slightly under predict the slope of the $\Delta\eta$ dependence.

We find that all the correlators exhibit a significant $\Delta\eta$ dependence except $C_{2,2,4}$ and $C_{2,3,5}$ which vary by only $20\%$ Adamczyk et al. (2017). The variation of $C_{m,n,m+n}$ with $\Delta\eta$ 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 $\Delta\eta$ provide the best opportunity for comparison of $\Delta\eta$-integrated quantities with hydrodynamic models.

In Fig. 3 we show centrality dependence of $\Delta\eta$-integrated $C_{m,n,m+n}$. We multiply the quantity $C_{m,n,m+n}$ by $N_{\rm part}^{2}$ 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}\approx{\langle\varepsilon_{m}\varepsilon_{n}\varepsilon_{m+n}\cos(m\Phi_{m}+n\Phi_{n}-(m+n)\Phi_{m+n})\rangle}$ from purely initial state geometry using a Monte-Carlo Glauber model Schenke et al. (2014). We find that the Glauber model predicts negative values for all combinations of $C_{m,n,m+n}$ ⁴⁴4Our calculations are consistent with the estimation of plane correlations performed in Ref. Teaney and Yan (2014). Since only a fraction of the initial state geometry is converted to final state anisotropy, i.e., $v_{n}\lesssim 0.1\times\varepsilon_{n}$  Teaney and Yan (2012), one therefore expects ${\langle v_{m}v_{n}v_{m+n}\cos(m\Psi_{m}+n\Psi_{n}-(m+n)\Psi_{m+n})\rangle}\lesssim 10^{-3}\times{\langle\varepsilon_{m}\varepsilon_{n}\varepsilon_{m+n}\cos(m\Phi_{m}+n\Phi_{n}-(m+n)\Phi_{m+n})\rangle}$, we therefore scale the Glauber model calculations by factors of $\sim 10^{-3}\!-\!10^{-4}$ to make a consistent data to model comparison Teaney and Yan (2012).

We compare our results with three different boost-invariant 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 $\eta/s=1/4\pi$ with MC-Glauber initial conditions by Teaney and Yan Teaney and Yan (2011, 2014), 2) hydrodynamic simulations MUSIC with boost invariant IP-Glasma initial conditions Schenke et al. (2012, 2010) that include a constant $\eta/s=0.06$ and a temperature dependent bulk viscosity $\zeta/s~(T)$ Ryu et al. (2015) and UrQMD afterburner Bass et al. (1998), 3) the perturbative-QCD$+$saturation$+$hydro based “EKRT” model Niemi et al. (2016) that uses two different parameterizations of the viscosity with constant $\eta/s=0.2$ and temperature dependent $\eta/s~(T)$ with a minimum of $\left(\eta/s~(T)\right)_{\rm min}=1.5/4\pi$ at a corresponding transition temperature between a QGP and hadronic phase of $T_{c}=150$ MeV and 4) viscous hydrodynamic model v-USPhydro Noronha-Hostler et al. (2014, 2013) with event-by-event TRENTO initial conditions  Moreland et al. (2015) tuned to IP-Glasma Schenke et al. (2012), that uses $\eta/s=0.05$, a freeze-out temperature of $T_{FO}=150$ MeV Alba et al. (2017a) and the most recent 2+1 flavors equation of state from the Wuppertal Budapest collaboration Borsanyi et al. (2014) combined to all known hadronic resonances from the PDG16+ Alba et al. (2017b).

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 Teaney and Yan (2011, 2014). 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}=\langle\cos(\phi_{a}+\phi_{b}-2\phi_{c})\rangle$. 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. Teaney and Yan (2011) and initial state geometry. It may also indicate a significant contribution from momentum conservation  Pratt (2010); Bzdak et al. (2011). 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 ($\sim-v_{2}/N$, $N$ being the multiplicity) Pratt (2010); Bzdak et al. (2011). 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 ⁵⁵5In the mid-central events we find $C_{1,2,3}$ to be positive at low transverse momentum ($p_{T1}<1$GeV) Adamczyk et al. (2017). of the first harmonic anisotropy correlated with the triangularity as was first predicted in Ref. Teaney and Yan (2011). 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_{\rm 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 Longacre (2016).

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}\approx\langle v_{2}^{2}v_{4}\cos(4(\Psi_{2}-\Psi_{4}))\rangle$ 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 $\varepsilon_{4}$ and $\varepsilon_{2}$, as predicted by several theoretical calculations Qiu and Heinz (2012); Teaney and Yan (2014); Bhalerao et al. (2013); Yan and Ollitrault (2015); Qian et al. (2016). This result is also qualitatively consistent with the measurement by the ATLAS collaboration at LHC Aad et al. (2014a); Jia (2014).

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 $\eta/s$ significantly underestimate $C_{2,2,4}$. The predictions using EKRT with a temperature dependent $\eta/s$ are much closer to the data; the same using constant $\eta/s$ under predict data by about $20\%$. A similar trend is also observed for $C_{2,3,5}$ shown in Fig. 3 (d). Although all hydrodynamic models shown in this figure predict correct qualitative trends of the centrality dependence, they all significantly underestimate the magnitude of $C_{2,3,5}$. Such discrepancy for EKRT has been argued Eskola et al. (2017) to be related to large off-equilibrium correlations which depend on the details of the parameterization $\eta/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 $\eta/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) Ryu et al. (2015). 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 Alba et al. (2017a), it effectively introduces a different viscous effect by choosing a lower freeze-out temperature $T_{FO}=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 Alba et al. (2017b). 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 $\Delta\eta$ and the centrality dependence of $C_{m,n,m+n}$ in Ref. Sun and Ko (2017). Such data-model comparisons demonstrate that the longitudinal structure of the initial state, global momentum conservation and multi-phase transport can capture the underlying dynamics that drives anisotropic flow and mixed-harmonic correlations Sun and Ko (2017).

Summary : We presented the first measurements of the charge inclusive three-particle azimuthal correlations $C_{m,n,m+n}=\langle\langle\cos(m\phi_{a}+n\phi_{b}-(m+n)\phi_{c})\rangle\rangle$ as a function of centrality, relative pseudorapidity and harmonic numbers $m,n$ in $\sqrt{s_{NN}}=$ 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 $\eta/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 $\Delta\eta$ 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 $\Delta\eta$ are large, $C_{2,2,4}$ and $C_{2,3,5}$ varies by only 20% between $\Delta\eta=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.

Acknowledgement : We thank Gabriel Denicol, Jacquelyn Noronha-Hostler, Harri Niemi, Risto Paatelainen, Björn Schenke, Chun Shen, Yifeng Sun and Li Yan for providing their model predictions and helpful discussions. We thank the RHIC Operations Group and RCF at BNL, the NERSC Center at LBNL, and the Open Science Grid consortium for providing resources and support. This work was supported in part by the Office of Nuclear Physics within the U.S. DOE Office of Science, the U.S. National Science Foundation, the Ministry of Education and Science of the Russian Federation, National Natural Science Foundation of China, Chinese Academy of Science, the Ministry of Science and Technology of China and the Chinese Ministry of Education, the National Research Foundation of Korea, GA and MSMT of the Czech Republic, Department of Atomic Energy and Department of Science and Technology of the Government of India; the National Science Centre of Poland, National Research Foundation, the Ministry of Science, Education and Sports of the Republic of Croatia, RosAtom of Russia and German Bundesministerium fur Bildung, Wissenschaft, Forschung and Technologie (BMBF) and the Helmholtz Association.