Normalized symmetric cumulants as a measure of QCD phase transition: a viscous hydrodynamic study

TL;DR

This study demonstrates that normalized symmetric cumulants of charged pions can distinguish between different QCD equations of state in viscous hydrodynamic models, aiding in identifying the QCD critical point.

Contribution

The paper introduces the use of normalized symmetric cumulants as a novel observable to differentiate QCD equations of state in hydrodynamic simulations.

Findings

Normalized symmetric cumulants distinguish between crossover and first-order phase transition EoSs.

Results are consistent across various initial conditions and shear viscosities.

These cumulants can serve as an EoS meter in experimental analyses.

Abstract

Finding the existence and the location of the QCD critical point is one of the main goals of the RHIC beam energy scan program. To make theoretical predictions and corroborate with the experimental data requires modeling the space-time evolution of the matter created in heavy-ion collisions by dynamical models such as the relativistic hydrodynamics with an appropriate Equation of State (EoS). In the present exploratory study, we use a viscous 2+1 dimensional event-by-event (e-by-e) hydrodynamic code at finite baryon densities with two different EoSs (i) Lattice QCD + HRG with a crossover transition and (ii) EoS with a first-order phase transition to studying the normalized symmetric cumulants of charged pions $v_n$ $(n=2-4)$. We show that the normalized symmetric cumulants can differentiate the two EoSs while all other conditions remain the same. The conclusion does not change for various initial conditions and shear viscosity. This indicates that these observables can be used to gain information about the QCD EoS from experimental data and can be used as an EoS meter.

Equations

$${\partial }_{\mu }{T}^{\mu \nu }$$

$${\partial }_{\mu }{J}^{\mu }$$

$$T^{\mu \nu }$$

$$J^{\mu }$$

$${\Delta }_{\alpha \beta }^{\mu v}{\tau }_{\pi }D{\pi }^{\alpha \beta }+{\pi }^{\mu v}=$$

$${\partial }_t(A)+{\partial }_x(v_{x}A)+{\partial }_y(v_{y}A)=B(t,x,y)$$

$$n$$

$$ϵ$$

$${\tilde{v}}_{\perp }=\frac{T^{\tau \perp }}{T^{\tau \tau }+P(ϵ)},$$

$${\tilde{v}}_i={\tilde{v}}_{\perp }\frac{T^{\tau i}}{T^{\tau \perp }}$$

$${ϵ}_p=\frac{\int d^{2}x\left(T^{xx}-T^{yy}\right)}{\int d^{2}x\left(T^{xx}+T^{yy}\right)}$$

$$\mathrm{SC}(m,n)\equiv {⟨⟨\cos \left(m{\phi }_1+n{\phi }_2-m{\phi }_3-n{\phi }_4\right)⟩⟩}_c=$$

$$-⟨⟨\cos \left[m\left({\phi }_1-{\phi }_2\right)\right]⟩⟩⟨⟨\cos \left[n\left({\phi }_1-{\phi }_2\right)\right]⟩⟩$$

$$=$$

$$\mathrm{NSC}(m,n)\equiv \frac{\mathrm{SC}(m,n)}{⟨v_m^2⟩⟨v_n^2⟩}.$$

$$\rho (x,y,z)=\frac{{\rho }_0}{1+\exp [(r-R)/a]},$$

$$T_{A,B}=\int dz{\rho }_{A,B}(x,y,z).$$

$$f\sim (T_A+T_B)+\alpha T_A{T}_B.$$

$$f(c{T}_A,c{T}_B)=cf(T_A,T_B),$$

$$f=T_R\left(p;T_A,T_B\right)\equiv {\left(\frac{T_A^p+T_B^p}{2}\right)}^{1/p}.$$

$${ϵ}_n{e}^{in{\Phi }_n}=-\frac{\int dxdy{r}^n{e}^{in\varphi }\epsilon \left(x,y,{\tau }_0\right)}{\int dxdy{r}^n\epsilon \left(x,y,{\tau }_0\right)}.$$

$${\Phi }_n=\frac{1}{n}\arctan \frac{\int dxdy{r}^n\sin (n\varphi )\epsilon \left(x,y,{\tau }_0\right)}{\int dxdy{r}^n\cos (n\varphi )\epsilon \left(x,y,{\tau }_0\right)}+\pi /n.$$

$$\frac{d^{3}N}{p_{T}d{p}_{T}d\varphi dy}=\frac{d^{2}N}{2\pi p_{T}d{p}_{T}dy}\sum _{n=-\infty }^{n=\infty }{v}_n(y,p_T)e^{in\left(\varphi -{\Psi }_n(y,p_T)\right)}.$$

$${\Psi }_n=\frac{1}{n}\arctan \left(\frac{⟨p_T\sin n\varphi ⟩}{⟨p_T\cos n\varphi ⟩}\right).$$

$$c(x,y)={⟨\frac{\left(x-⟨x{⟩}_{\mathrm{ev}}\right)\left(y-⟨y{⟩}_{\mathrm{ev}}\right)}{{\sigma }_x{\sigma }_y}⟩}_{\mathrm{ev}},$$

$$v_2=m{ϵ}_2+\delta ,$$

Full text

††preprint: APS/123-QED

I Introduction

It is well known that at low temperature and small baryon chemical potential, the degrees of freedom of nuclear matter are color-neutral hadrons; whereas at high temperature or large baryon chemical potential, the matter is in the form of quark-gluon plasma (QGP), in which the fundamental degrees of freedom are colored objects quarks and gluons. The nuclear matter at high baryon density and finite temperature are believed to undergo first-order phase transitions, from the hadronic phase to the QGP phase, and the first-order phase transition line terminates at a critical point Datta et al. (2013); Borsanyi et al. (2010); Gavai and Gupta (2005). This is because lattice QCD shows that the hadron to QGP transition is a crossover for vanishing baryon chemical potential at temperature $\sim$ 170 MeV Bernard et al. (2005); Bhattacharya et al. (2014); Borsanyi et al. (2010); Bazavov et al. (2012). For details on QCD phase diagram and critical points see review Fukushima and Hatsuda (2011); Baym et al. (2018)

Present theoretical models widely disagree with each other regarding the value of critical temperature and baryon chemical potential corresponding to the QCD critical point on the QCD phase diagram Schaefer et al. (2007); Kovacs and Szep (2008). Also, the existence of the QCD critical point is yet to be confirmed experimentally Aggarwal et al. (2010); Adamczyk et al. (2014a, b, 2018); Luo (2015); Mohanty (2009). It is crucial that phenomenologically motivated studies of the heavy-ion collision, such as relativistic hydrodynamics, be capable of accounting for the potential influence of such a critical point on experimental observables. Recently it was shown Pang et al. (2018) that the goal mentioned above can be achieved by using relativistic hydrodynamic model, experimental data and a state-of-the-art deep-learning technique that uses a convolutional neural network to train the system.

The present exploratory study aims to find a unique observable which connects QCD Equation of State (EoS) and the experimental data of heavy-ion collisions using one of the available dynamical models. We believe this effort will be complementary to the finding of Pang et al. (2018). It is well known that the initial energy/baryon densities follow very different space-time trajectories in the QCD phase diagram for crossover and first-order phase transitions Stephanov et al. (1998). This fact is shown in Fig. 1(a), the three different lines correspond to different entropy per baryon numbers ($s/n_{B}$), and the trajectory for the crossover will follow $\mu_{B}=0$ line. The first-order phase transition involves finite latent heat (kink in the trajectories), and the speed of sound is zero in the mixed-phase. We note that in the fluid dynamical picture, converting the initial spatial deformation to the final momentum anisotropy depends on the speed of sound (EOS) and other factors such as the viscosity of the medium. In the case of crossover, the speed of sound is never zero; the consequence of the two different EoS’s can be seen in the temporal evolution of momentum space anisotropy $\epsilon_{p}$ (defined later) shown in Fig. 1(b). A different value of $\epsilon_{p}$ corresponds to different elliptic flow. Although event averaged elliptic flow is sensitive to the EoS used, it is also known that its value is suppressed in the presence of finite shear viscosity. Therefore the event-averaged elliptic flow is not a good indicator of the EoS. Nevertheless, it is interesting to investigate what imprint of different EoSs one can find in the correlation between the initial fluctuating geometry and the final flow coefficients in the event-by-event collisions Niemi et al. (2013).

We use the relativistic hydrodynamic model which has been very successful in simulating heavy-ion collisions and explaining experimental observablesKolb et al. (2000); Nonaka et al. (2000); Hirano (2002); Aguiar et al. (2001); Chaudhuri (2006); Nonaka and Bass (2007); Song and Heinz (2008a); Baier and Romatschke (2007); Muronga (2007); Pratt and Vredevoogd (2008); Petersen et al. (2008); Molnar et al. (2010); Schenke et al. (2010); Werner et al. (2010); Holopainen et al. (2011); Roy and Chaudhuri (2012); Bozek (2012); Pang et al. (2012); Akamatsu et al. (2014); Noronha-Hostler et al. (2013); Del Zanna et al. (2013). We find the linear/Pearson correlation (defined later) of initial geometric asymmetry to the corresponding flow coefficient (particularly the second-order flow coefficient $v_{2}$) is a unique observable which can differentiate between EoS with a first-order phase transition to that with a crossover transition irrespective of the initial condition used. It has been known that the event averaged $v_{2},$ and the eccentricity of the averaged initial state, $\epsilon_{2}$ are approximately linearly correlated Niemi et al. (2013); Plumari et al. (2015); Kolb and Heinz (2003); Song and Heinz (2008b),and the initial condition used Niemi et al. (2013). Although the above finding sounds promising, we need to keep in mind that the initial eccentricities are not directly measurable in experiments; hence the correlation between the eccentricities and flow coefficients may not provide any practical information about the EoS’s. But we can use the cumulants between flow coefficients such as $\mathrm{NSC}(v_{2},v_{3})$ and $\mathrm{NSC}(v_{3},v_{4})$ Adam et al. (2016); Bilandzic et al. (2014) which can always be measured experimentally without referring to any particular model. That is what motivates us to investigate these observables described above as an EoS meter.

For the present study, we use a newly developed $2+1$-dimensional event-by-event viscous hydrodynamic code with an EoS Kolb et al. (2000) with finite baryon chemical potential/density for low $\sqrt{s_{NN}}$ collisions, and a lattice QCD+HRG EoS Huovinen and Petreczky (2010) for the higher $\sqrt{s_{NN}}$ collisions. There are ongoing efforts to construct EoS with a critical point Plumberg et al. (2018); Monnai et al. (2019). Due to the present uncertainty in the location of QCD critical point, we refrain to use sophisticated EoSs.

The conservation equations are solved numerically by using the time-honored SHarp And Smooth Transport Algorithm (SHASTA) Boris and Book (1973). In the next section, we shall discuss the details of various tests performed to find the numerical accuracy of the code. As we will show, our code passes all the test cases satisfactorily.

The paper is organized as follows: in the next section, we discuss various aspects of the newly developed hydrodynamic code. In section II we show the comparison of experimental data of charged pion invariant yield and elliptic flow in Au+Au 200 GeV collisions to the simulation result for optical Glauber initial condition to test the code. In the same section, we also discuss the correlation between different observable and the effect of EoS on them. Finally in section III we conclude and discuss some future possibilities. Throughout this article we adopt the units $\hbar=c=k_{\mathrm{B}}=1.$ The signature of the metric tensor is always taken to be $g^{\mu\nu}=\operatorname{diag}(+1,-1,-1,-1)$. Upper greek indices correspond to contravariant and lower greek indices covariant. The three vectors are denoted with Latin indices.

I.1 Conservation equations

Relativistic hydrodynamics model of high-energy heavy-ion collisions assumes that matter created in collisions reaches a state of local thermal equilibrium at time $\tau_{0}$. The evolution of the thermalised nuclear matter is governed by the conservation equation of energy-momentum tensor and net baryon current,

[TABLE]

where the the energy-momentum tensor and the net baryon current can be expressed as

[TABLE]

Here $\epsilon$ is the local energy density, $P$ is pressure, $g^{\mu\nu}$ is the metric tensor, $n$ is the net baryon density, $u^{\mu}$ is the time-like 4-velocity with $u^{\mu}u_{\mu}=1$ and $\pi^{\mu\nu}$ is the shear-stress tensor. The evolution equation for the shear-stress tensor is given as Denicol et al. (2012); Israel and Stewart (1976)

[TABLE]

where $\eta$ is the shear viscosity coefficient, $D=u^{\mu}\partial_{\mu}$ is the comoving time derivative, $\sigma^{\mu v}=\Delta_{\alpha\beta}^{\mu v}\partial^{\alpha}u^{\beta}$ is the shear tensor, $\theta=\partial_{\mu}u^{\mu}$ is the expansion rate, and $\Delta_{\alpha\beta}^{\mu v}=\left(\Delta_{\alpha}^{\mu}\Delta_{\beta}^{v}+\right.$ $\left.\Delta_{\alpha}^{v}\Delta_{\beta}^{\mu}-2/3\Delta^{\mu v}\Delta_{\alpha\beta}\right)/2$, with $\Delta^{\mu v}=g^{\mu v}-u^{\mu}u^{v}$. The transport coefficients in the above equation were obtsined in the massless limit with the relaxation time $\tau_{\pi}$ being $\tau_{\pi}=5\eta/(\epsilon+P)$.

In the present work, we will be using two kinds of EoS’s:

(i) A parameterized EoS Noronha-Hostler et al. (2019) shown in Fig. 2 (left panel) which has a cross-over transition between high temperature QGP phase obtained from lattice QCD and a hadron resonance gas below the crossover temperature. In order to have a smooth pressure profile as a function of $\varepsilon$ and $n$, the pressure is first expanded as a as a Taylor series in powers of $\mu_{B}/T$. The associated Taylor coefficients which are known form LQCD Borsanyi et al. (2014); Bazavov et al. (2014) are parameterized using polynomial of ninth order. (ii) A parameterized EoS Kolb et al. (2000) with a first order phase transition at finite net baryon number density shown in Fig. 2 (right panel). The EoS connects a non-interacting massless QGP gas at high temperature to a hadron resonance gas, of masses up to 2 GeV, at low temperatures through a first order phase transition. The bag constant $B$ is a parameter adjusted to $B^{1/4}=230~{}$MeV, to yield a critical temperature of $T_{c}=164~{}$MeV. We note that the choice of bag parameter used here is not unique, it may vary between $B^{1/4}=150-300$ MeV Baym et al. (2018). Both of these EoS will be used in the hydrodynamic simulation corresponding to $\sqrt{s_{\mathrm{NN}}}=$ 62.4 GeV collisions at a finite net baryon density.

Instead of solving the full $3+1-$ dimensional hydrodynamics equations, we consider a simplified evolution by assuming boost invariant expansion in the $z$ direction. This can be done easily by working in the Milne coordinates with metric tensor given by $g^{\mu\nu}=(1,-1,-1,-1/\tau^{2})$. The velocity 4-vector in this coordinate (physical quantities in this coordinate is written with a tilde sign), $\tilde{U}^{\mu}=\tilde{\gamma}(1,{\tilde{v}_{x}},{\tilde{v}_{y}},0)$, where $\tilde{v}_{i}=v_{i}\cosh\eta,\hskip 2.84544pt(i=x,y)$ and $\tilde{\gamma}=1/\sqrt{1-\tilde{v}_{\bot}^{2}}$ with $\tilde{v}_{\bot}^{2}={\tilde{v}_{x}}^{2}+{\tilde{v}_{y}}^{2}$ .

Conservative equations of the form Eqs. (1,2,5) can be solved accurately using flux-corrected transport (FCT) algorithms Refs. Boris et al. (1993); Powell et al. (1999); Marti and Mueller (1999) without violating the positivity of mass and energy, particularly near shocks and other discontinuities. Here, we use the SHASTA (”SHarp and Smooth Transport Algorithm”), which is a designed to solve partial differential equations of the form

[TABLE]

where $A=A(t,x,y)$ is for example $T^{00},T^{0i},..,,v_{i}$ is the $i^{\mathrm{th}}$ component of three-velocity, and $B(t,x,y)$ is a source term. The local rest frame charge and energy densities are given as

[TABLE]

while the velocity components are calculated using an one-dimensional root search algorithm, by iterating the following transcendental equation

[TABLE]

such that the components of velocity are given by

[TABLE]

where $T^{\tau\bot}=\sqrt{{(T^{\tau x})}^{2}+{(T^{\tau y})}^{2}}$. In this velocity finding algorithm we keep the accuracy $\sim 10^{-22}$.

We have tested our numerical simulation against $(1+1)-$ dimensional Riemann problem which has an analytical solution in the perfect-fluid limit. The numerical solution reproduces the analytic solution with a very high precision with sufficiently high numerical resolution. We have similarly checked and verified the applicability of our numerical simulation to $2+1-$ dimension by having the initial discontinuity in the $1+1-$ dimensional Riemann problem in the plane $y=-x$. For small viscosities $\eta/s=10^{-2}$, we have verified that the viscous version of our simulation reproduces the result of our simulation with $\eta/s=0$ but with anti-diffusion coefficient $A_{\mathrm{ad}}=0.9$, which is the numerical analog of physical viscosity. Finally, we have checked that our numerical simulation reproduces Gubser flow with initial condition fixed such that, the system expands and cools down much faster than typical heavy-ion collisions. In the numerical code we use the CORNELIUS subroutine reported in Huovinen and Petersen (2012). The subroutine uses an improved version of the original Marching Cube algorithm and extends the different distinct possible topological configuration from $15$ to $33$ and thus creates a consistent surface with lesser holes and no double counting which is required for event-by-event hydrodynamics study.

II Result

II.1 Testing for smooth Glauber

Before we show the result of event-by-event hydrodynamics simulation we show here the results from smooth optical Glauber  model for $\sqrt{s_{NN}}=200$, $62.4$ GeV collision energies respectively. For all the calculations, the spatial extension of the numerical grid is set to $36\times 36$ fm2. The spatial grid spacing is set to $\Delta x=\Delta y=0.09$ fm and the temporal spacing is set to $\Delta t=0.04$ fm. The initial time for all collisional erergies is fixed to $\tau_{0}=0.6$ fm.

For Au-Au collisions at $\sqrt{s_{NN}}=200$ GeV we use EoS-Lattice shown in Fig. 2(a) with $\mu_{B}=0$ MeV and the freeze-out energy density is taken to be $\varepsilon_{F}=0.28$ GeV/fm3 which corresponds to a temperature $T_{f}=137$ MeV. The central energy density is fixed to $\varepsilon_{0}=45$ GeV/fm3 for $10\%-15\%$ collision centrality which explains the corresponding experimental data.  The top panel of Fig. 3 shows the comparison of invariant yield  of $\pi^{-}$ for $10\%-15\%$ and $40\%-50\%$ centrality collisions. Experimental data measured by PHENIX collaboration Adler et al. (2004) are shown by different symbols and the corresponding simulation results are shown by lines. We found reasonable agreement with the experimental data except at low momentum $p_{T}<0.5$ GeV, this is because we neglected the  contribution of $\pi^{-}$ yield coming from various heavier resonance decay. We plan to incorporate the resonance  decay in a later version of the code. The comparison of simulated elliptic flow with the experimental data of STAR collaboration Adams et al. (2005) are shown in the bottom panel of Fig. 3. We can see that the simulated result nicely explain the observed asymmetry in the final momentum spectra of $\pi^{-}$.

For Au-Au collisions at $\sqrt{s_{NN}}=62.4$ GeV we use EoS-$1^{st}$ order PT shown in Fig. 2(b) and the freeze-out energy density is fixed to $\varepsilon_{F}=0.3$ GeV/fm3. The initial central energy density is set to $\varepsilon_{0}=16$ GeV/fm3 for $0\%-15\%$ collision centrality. Similarly, the initial central net baryon density is fixed to $n_{0}=0.4$ fm*-3*. We have checked that the above parameters explains the invariant yield of $\pi^{-}$Back et al. (2007) across various centralities. In Fig. 1(a) we have plotted the trajectories of different regions of the fireball during hydrodynamical evolution as a function temperature ($T$) and net baryon chemical potential ($\mu_{B}$) with initial energy densities corresponding to 100$\%$, 50$\%$, 25$\%$ of the maximum energy density $\epsilon_{0}=16.2$ Gev/fm3 till freeze-out $\epsilon_{F}$. The corresponding constant $s/n_{B}$ for the above lines are 156, 175 and 206 respectively. Momentum space anisotropy $\epsilon_{p}$ defined as

[TABLE]

where $T^{ii}$ are components of energy-momentum tensor $T^{\mu\nu}$ are plotted in Fig. 1(b) both for EoS-Lattice and first order phase transition at $b=8$ fm. This completes the testing of our code. In the next section, we shall study the effect of phase transition on the various flow correlation which is the main goal of this paper.

II.2 Normalised symmetric cumulant and event-by-event flow correlation

Flow observables obtained from multiparticle correlations, symmetric cumulants (SC), are introduced, which are defined as Adam et al. (2016); Bilandzic et al. (2014)

[TABLE]

with the condition $m\neq n$ for two positive integers $m$ and $n$. The double angular brackets indicate that the averaging procedure has been performed in two steps- first over all distinct particle quadruplets in an event, and then in the second step the single-event averages were weighted with number of combinations. The four-particle cumulant in Eq. (12) is less sensitive to non-flow correlations than any 2- or 4-particle correlator on the right-hand side taken individually Borghini et al. (2001a, b). The last equality is true only in the absence of non-flow effects. The observable in Eq. (12) is zero in the absence of flow fluctuations, or if the magnitudes of harmonics $v_{m}$ and $v_{n}$ are uncorrelated. Experimentally it is more reliable to measure the higher order moments of flow harmonics $v^{k}_{n}$ $(k\geq 2)$ with 2- and multiparticle correlation techniques, than to measure the first moments $v_{n}$ with the event plane method, due to systematic uncertainties involved in the event-by-event estimation of symmetry planes.

One can also define normalised symmetric cumulants (NSC) as,

[TABLE]

Normalized symmetric cumulants reflect the strength of the correlation between $v_{m}$ and $v_{n}$, while $\mathrm{SC}(m,n)$ has contributions from both the correlations between the two different flow harmonics and the individual harmonics. In the present work, we calculate the symmetric cumulants using the last equality in Eq. (12), due to the absence of non-flow effects in our hydrodynamic formulation.

The hydrodynamics equation of motion has to be initialized at an initial time $\tau=\tau_{0}$, by specifying the initial energy densities $\varepsilon(\tau_{0},x,y)$ or the entropy densities $s(\tau_{0},x,y)$ and the flow velocities $u^{\mu}$ of a relativistic fluid. Two of such models that we are going to use are the Glauber model Miller et al. (2007) and the TRENTo model Moreland et al. (2015).

Given a pair of projectiles labeled $A$ and $B$ collide along the beam axis $z$ and let $\rho_{A,B}(x,y,z)$ be the density of the nuclear matter that participate in inelastic collisions. $\rho(x,y,z)$ is usually given by the Woods-Saxon profile, given as

[TABLE]

where $\rho_{0},R$ and $a$ are the normalization, size of the nucleus, and the stiffness of the edge of nucleon distribution profile respectively while $x,y,z$ are the spatial coordinates. Each projectile may then be represented by its participant thickness, given by

[TABLE]

We shall assume that their exists a function $f(T_{A},T_{B})$ which converts the projectile thickness to entropy production, namely $f\propto ds/dy|_{\tau=\tau_{0}}$. In a two-component Glauber model the function is the sum of $f\sim T_{A}+T_{B}$, which is proportional to the number of wounded nucleons $N_{WN}$ and a quadratic term $f\sim T_{A}T_{B}$, which is proportional to the number of binary collisions $N_{BC}$. The complete function $f$ is given as

[TABLE]

The proportionality constants and the other relevant details can be found in Miller et al. (2007). However, in general the energy density profile in the reaction zone of nucleus-nucleus collision fluctuates from event to event due to the quantum fluctuations of the nuclear wave function. This fluctuations are attributed to the fluctuations in the positions of participating nucleons. In such a case, the previous decribed smooth Glauber distribution will then be an ensemble average of a large number of fluctuating initial distribution. This can be achieved by extending the smooth Glauber model to their Monte Carlo (MC) versions. These fluctuating initial distributions breaks the rotational and reflection symmetry of the smooth distributions.

The TRENTo model Moreland et al. (2015) for the initial condition assumes the a scale invariant form of the function $f$, i.e.,

[TABLE]

for any nonzero constant $c$. One can clearly see that the above is broken by the binary collision term $\alpha T_{A}T_{B}$ in the Glauber model. TRENTo assumes a reduced thickness function given by

[TABLE]

Various limiting cases of the function $f$, for different choices of parameter $p$ can be found in Moreland et al. (2015). In the present study we assume $p=0$, in which case $f=\sqrt{T_{A}T_{B}}$, which is the geometric mean of the $T_{A}$ and $T_{B}$. We will be using the above two models in order to see the sensitivity of various correlations to the initial conditions used for the simulations.

We are interested in the study of the difference in fluid dynamical response of the system to the initial geometry (anisotropy) for two different EoS. In the literature this initial geometry/anisotropy of the overlap zone of two colliding nucleus is quantified in terms of eccentricities $\epsilon_{n}$ Alver and Roland (2010); Alver et al. (2010):

[TABLE]

where $r^{2}=x^{2}+y^{2},\phi$ is the spatial azimuthal angle, and $\Phi_{n}$ is the participant angle given by,

[TABLE]

$\epsilon_{n}$ is basically the eccentricity of a polygon of $n^{\mathrm{th}}$ order, which can be reconstructed from a initial distribution generated from a MC Glauber or TRENTo model in a given event. Such a $n^{\mathrm{th}}$ order initial distribution generates a $n^{\mathrm{th}}$ order harmonic flow $v_{n}$ analogous to $v_{2}$.

As described before, the azimuthal momentum anisotropy is characterized in terms of the coefficients $v_{n}$ of the Fourier expansion of the single particle azimuthal distribution. We use the following definition to calculate $v_{n}$ in each event:

[TABLE]

In Eq. (21) $v_{n}$ can be calculated as the expectation $v_{n}=\langle\cos(n\phi-n\Psi_{n})\rangle$, with respect to the associated event plane angle $\Psi_{n}$ of $n$-th order harmonics and is defined as

[TABLE]

The angle brackets indicates the averaged value with respect to the particle spectrum.

Event-by-event hydrodynamics is one of the most natural way to model azimuthal momentum anisotropies $v_{n}$ Eq. (21) generated by fluctuating initial-state anisotropies $\epsilon_{n}$ Eq. (19), which are generated by the highly fluctuating initial conditions in experiments. The largest source of uncertainty in these hydrodynamic models are the initial conditions. Since, direct measurement of bulk properties of matter like EoS in experiments is not possible, we try to identify hydrodynamic responses which in one hand can be calculated in experiments and in other hand are also tolearant to uncertainities in the model parameters viz., the initial conditions, shear viscosity etc.

It has been known that the event averaged $v_{n}$, and the eccentricity of the averaged initial state, $\epsilon_{n}$ are approximately linearly related Kolb and Heinz (2003); Song and Heinz (2008b) for $n<4$ but the same may not be true for higher-order flow coefficients. It is also known that the same linear relationship holds approximately between $\epsilon_{n}$ and $v_{n}$ even for event-by-event fluctuating conditions Niemi et al. (2013). Here we first study how the event-by-event correlation between $\epsilon_{n}$ and $v_{n}$ is changed when the system undergoes either a first phase transition or a cross over. In order to quantify the linear correlation we use Pearson’s correlation coefficient which is defined as

[TABLE]

where $\sigma_{x}$ and $\sigma_{y}$ are the standard deviations of the quantities $x$ and $y$. The correlation coefficient ranges from $-1$ to $1$. A value of $1(-1)$ implies that a linear (anti-linear) correlation between $x$ and $y$. A value of [math] implies that there is no linear correlation between the variables. In a previous study, it was shown that the Pearson correlator $c(\epsilon_{2},v_{2})$ is almost insensitive to the different initial condition and the value of shear viscosity over entropy density of the fluid Niemi et al. (2013). Further, assuming an approximate linear relationship between the $\epsilon_{2}$ and $v_{2}$, we can write for event-by-event case

[TABLE]

where $m=\left\langle v_{2}\right\rangle_{\mathrm{ev}}/\left\langle\epsilon_{2}\right\rangle_{\mathrm{ev}},$ and the average error $\left\langle\delta\right\rangle_{\mathrm{ev}}=0.$ The values of $m$ indicate how efficiently the initial deformation is transformed into the final momentum anisotropy.

In Fig. (4) (a) (Top row) we show the event-by-event distribution of $v_{2}$ vs $\epsilon_{2}$ for $0\%-5\%$ Au+Au collisions at $\sqrt{s_{NN}}=62.4$ GeV. The initial energy density and $\epsilon_{2}$ are obtained from the MC-Glauber model with the contribution coming only from wounded nucleons. The result is obtained for EoS-Lattice, i.e, crossover transition. Fig. (4) (b) (Top row) shows the same, but with EoS having a first order phase transition. Using two different EoS we found $\sim 10\%$ decrease in $c(\epsilon_{2},v_{2})$ for the case of first order phase transition, which clearly indicates that $c(\epsilon_{2},v_{2})$ can be treated as a good signal of phase transition in the nuclear matter. Above results are not surprising since the speed of sound becomes zero (hence the expansion) for a certain temperature range in the first order phase transition. Fig. (4) (a,b) (Middle row and bottom row) shows the same results but for $20\%-30\%$ centrality using MC-Glauber and TRENTo model initial conditions. The first order phase transition shows a $\sim 6\%$ and $\sim 5\%$ decrease in the value of $c(\epsilon_{2},v_{2})$ respectively. Thus, we infer that although, the value of $c(\epsilon_{2},v_{2})$ for first-order phase transition is always less than that with crossover transition independent of the model used, the difference is more prominent in the central collisions than at higher centrality. Similarly, other higher order correlations e.g., $c(\epsilon_{n}$ ,$v_{n})$ (for $n=3$ or $n=4$) is found to be smaller for the case of first order phase transition.

However, note that the initial eccentricities $\epsilon_{n}$ are not accessible in real experiments (and are model dependent) and hence the $c(\epsilon_{n},v_{m})$ are not as interesting as $c(v_{n},v_{m})$ which can be calculated from the available experimental data. As pointed out earlier instead of $c(v_{n},v_{m})$, a clean experimental observable would rather be the normalised symmetric cumulants $\mathrm{NSC}{(m,n)}$.

This EoS dependence of the flow correlations can be more clearly seen from NSC$(m,n)$ Fig. (5) (left panel), where we show $\mathrm{NSC}(2,3)$, $\mathrm{NSC}(2,4)$, and $\mathrm{NSC}(3,4)$ for EoS-Lattice (solid orange circles) and EoS first order phase transition (open blue circles) with corresponding errors for $20\%-30\%$ collision centrality. We have used the MC-Glauber $\varepsilon_{WN}$ initialisation. Fig.5 (right panel) shows the same but for TRENTo model. The errors are calculated by using bootstrap method. As can be seen from the Figs. (5) that the $\mathrm{NSC}(2,3)$, $\mathrm{NSC}(2,4)$, and $\mathrm{NSC}(3,4)$ always distinguishes the two different EoSs. These observations may be attributed to very different evolutionary dynamics of the system for the two different EoS, as the speed of sound becomes zero in first-order phase transition hence the linear/non-linear coupling of $\epsilon_{n}$ - $v_{n}$ and $v_{n}$-$v_{m}$ is different in the two scenario. Although the absolute values of the $\mathrm{NSC}(m,n)$ varies for different initial condition (energy density scales with wounded nucleons or in the TRENTo model), the difference in them remains almost same for two different EoSs. We found in the mid central collisions $\mathrm{NSC}(m,n)$ is larger for the EoS with first order phase transition irrespective of the initial conditions used here.

Finally, in order to see the influence of viscosity in the observables discussed above, we did the simulation for various values of $\eta/s$. In Fig. (6), we show the difference between $\mathrm{NSC}(m,n)$ calculated for EoS- with first order phase transition and that with crossover. The initial energy density is obtained from wounded nucleons $(\varepsilon_{WN})$ in MC-Glauber model for $20\%-30\%$ collision centrality. As in ideal hydrodynamics, the results of $\mathrm{NSC}(m,n)$ for first order phase transition are again larger than that with crossover. Increasing $\eta/s$ increases the difference between them for $\mathrm{NSC}(3,4)$. The same is not true for $\mathrm{NSC}(2,3)$, which is too sensitive to the shear viscosity and responds to it in a rather non-trivial manner.

This indicates that we can utilize $\mathrm{NSC}(3,4)$ to probe the EoS of the system which implies that one can possibly use this observable to locate the QCD critical point. For example we can calculate $\mathrm{NSC}(3,4)$ from available experimental data for various $\sqrt{s_{NN}}$ and pinpoint the energies where $c(v_{n},v_{m})$ shows a sudden change in magnitude.

III Conclusion and outlook

In this paper, we studied correlation $\mathrm{NSC}(m,n)$ and $c(\epsilon_{m},v_{n})$ for $m,n=2-4$ in ultra-relativistic heavy-ion collisions using event-by-event viscous fluid dynamics model for two different EoSs. The comparison of numerical results from the hydro code to the corresponding analytic solution in 1+1D and in 2+1D shows very good agreement between them. Using optical Glauber model as an initial condition, we confirm that our code nicely explains experimentally measured charged pion invariant yield and elliptic flow in Au+Au collision at $\sqrt{s_{NN}}=$ 200 GeV and 62.4 GeV. First using ideal hydrodynamics with two different EoS’s we showed that the $\mathrm{NSC}(2,3)$, and $\mathrm{NSC}(3,4)$ clearly distinguishes a first order phase transition scenario and a cross over transition. Although the absolute values of the $\mathrm{NSC}(m,n)$ varies for different initial condition (energy density scales with binary collisions or wounded nucleons) the difference in them remains almost same for two different EoSs. We found in the mid central collisions $\mathrm{NSC}(m,n)$ is larger for the EoS with first order phase transition irrespective of the initial conditions used here. In the presence of finite viscosity, the result that the value of $\mathrm{NSC}(m,n)$ is larger for EoS with first order phase still remains valid. However by inreasing the viscosity we found the difference between the results from first order phase transition and that with crossover increases for $\mathrm{NSC}(3,4)$.

This indicates that we may utilize $\mathrm{NSC}(3,4)$ to probe the EoS of the system which implies that one can possibly use this observable to locate the QCD critical point. For example, we can calculate $\mathrm{NSC}(m,n)$ from available experimental data for various $\sqrt{s_{NN}}$ and pinpoint the energies where $\mathrm{NSC}(m,n)$ shows a sudden change in magnitude.

The coefficients $\mathrm{NSC}(m,n)$ can be directly calculated from the available experimental data, whereas calculation of $c(\epsilon_{n},v_{m})$ requires input from an initial condition model hence proves to be less attractive observable than $\mathrm{NSC}(m,n)$. However, we would like to point out that our study also shows $c(\epsilon_{2},v_{2})$ is a promising observable for probing the EoS as it shows $\sim 10\%-30\%$ difference in the two scenarios.

Acknowledgement

We thank J.Y. Ollitrault and R. Bhalerao for suggestions and helpful discussion. V.R. is supported by the DST Inspire faculty research grant (IFA-16-PH-167), India. AD acknowledges financial support from DAE, Government of India.